	SUBROUTINE PLFILE(ISP,JFIL,STATUS)
C
C	PLOT OUT X-RAY DIFFRACTION DATA FOR COMPARISON OF UP TO 8 FILES
C
C adapted for PGPLOT graphics package 	BHT 12/86-1/87
C
	include 'cmpr.common'
	include '[-.shared]flabel.common'

	INTEGER ISP,JFIL(*)
	LOGICAL STATUS

	LOGICAL flag,TTTOQ

	LOGICAL*1  QUESTF,ONEBOX

	LOGICAL QF,DF,SQRTF,LINF,LOGF
	LOGICAL clip,GRIDF,TITLEF,FILELBL
	LOGICAL sizef
	LOGICAL STKF,STKSET,FIRSTF,INDIVF,HARDF
	LOGICAL XLIMLO,XLIMHI,YLIMLO,YLIMHI
	INTEGER ixstack, iystack
	REAL XLOW,XHIGH,YLOW,YHIGH
	REAL XLEN,YLEN,XMAR,YMAR

	BYTE IDO,fchar,ians

	real XT(2),YT(2)
	character LBLX*50,LBLY*50
	INTEGER ln,fn

	INTEGER IFIL(10),nfil
C......................................................................
	CALL GETPLQUAL(clip,QF,DF,SQRTF,LINF,LOGF,
	1	sizef,TITLEF,FILELBL,GRIDF,STKF,STKSET,FIRSTF,INDIVF,
	2	HARDF,XLIMLO,XLIMHI,YLIMLO,YLIMHI,ixstack,iystack,
	3	XLOW,XHIGH,YLOW,YHIGH,XLEN,YLEN,XMAR,YMAR,wavel)
	QUESTF = .false.
C......................................................................

	NFIL = ISP
	DO i=1,NFIL
	  IFIL(I) = JFIL(I)
	ENDDO
C check that the conversions are valid and get the minimum and maximum
C x values
	FLAG = .TRUE.
	DO 30 I=1,NFIL
	  fn= Ifil(I)
	  npout1 = 0
	  npout2 = 0
	  IF (PLTYP(fn) .GT. 0 .and. NPOINT(fn) .gt. 0) THEN
C try to convert one point -- see if we get an error; display message
	    IF (QF) THEN
	      call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_Q,npout,0.0,.true.,*21)
	    ELSEIF (DF) THEN
	      call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_D,npout,0.0,.true.,*21)
	    ELSEIF (WAVEL .GT. 0) THEN
	      call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_TT,npout,wavel,.true.,*21)
	    ENDIF
C conversion was OK
C   find first valid point
	    DO J=1,NPOINT(fn)
	      IF (QF) THEN
	       call XYCONV(BX(J,fn),B(J,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_Q,npout1,0.0,.false.,*23)
	      ELSEIF (DF) THEN
	        call XYCONV(BX(J,fn),B(J,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_D,npout1,0.0,.false.,*23)
	      ELSEIF (WAVEL .GT. 0) THEN
	        call XYCONV(BX(J,fn),B(J,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_TT,npout1,wavel,.false.,*23)
	      ELSE
	        npout1 = 1
	        xout = BX(J,fn)
	      ENDIF
	      IF (npout1 .gt. 0) THEN
	        IF (FLAG) THEN
	          FLAG = .FALSE.
	          XMIN = XOUT
	          XMAX = XOUT
		ELSE
		  XMIN = MIN(XMIN,XOUT)
		  XMAX = MAX(XMAX,XOUT)
	        ENDIF
	        GOTO 26
	      ENDIF
	    ENDDO
	  ENDIF

26	  IF (PPEAKS(fn) .GT. 0) THEN
	    DO J=1,NPKS(fn)
	      IF (QF) THEN
	        call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_Q,npout2,0.0,.false.,*23)
	      ELSEIF (DF) THEN
	        npout2 = 1
	        xout = DSPACE(J,fn)
	      ELSEIF (WAVEL .GT. 0) THEN
	        call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_TT,npout2,wavel,.false.,*23)
	      ELSE
	        npout2 = 1
	        xout = TT(J,fn)
	      ENDIF
	      IF (npout2 .gt. 0) THEN
	        IF (FLAG) THEN
	          FLAG = .FALSE.
	          XMIN = XOUT
	          XMAX = XOUT
		ELSE
		  XMIN = MIN(XMIN,XOUT)
		  XMAX = MAX(XMAX,XOUT)
	        ENDIF
	        goto 25
	      ENDIF
	    ENDDO
	  ENDIF

C no points in range
25	  IF (npout1+npout2 .le. 0) THEN
	    type *,'No valid points after conversion for file #',fn
	    Ifil(I) = 0
	    goto 30
	  ENDIF

C find last valid point
	  IF (PLTYP(fn) .GT. 0) THEN
	    DO J=NPOINT(fn),1,-1
	      IF (QF) THEN
	       call XYCONV(BX(J,fn),B(J,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_Q,npout,0.0,.false.,*23)
	      ELSEIF (DF) THEN
	        call XYCONV(BX(J,fn),B(J,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_D,npout,0.0,.false.,*23)
	      ELSEIF (WAVEL .GT. 0) THEN
	        call XYCONV(BX(J,fn),B(J,fn),FXTYP(fn),1,WAVE(fn),
	1		   XOUT,YOUT,XTYP_TT,npout,wavel,.false.,*23)
	      ELSE
	        npout = 1
	        xout = BX(J,fn)
	      ENDIF
	      IF (npout .gt. 0) THEN
	        XMIN = MIN(XMIN,XOUT)
	        XMAX = MAX(XMAX,XOUT)
	        goto 28
	      ENDIF
	    ENDDO
	  ENDIF
28	  IF (PPEAKS(fn) .GT. 0) THEN
	    DO J=NPKS(fn),1,-1
	      IF (QF) THEN
	        call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_Q,npout,0.0,.false.,*23)
	      ELSEIF (DF) THEN
	        npout = 1
	        xout = DSPACE(J,fn)
	      ELSEIF (WAVEL .GT. 0) THEN
	        call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_TT,npout,wavel,.false.,*23)
	      ELSE
	        npout = 1
	        xout = TT(J,fn)
	      ENDIF
	      IF (npout .gt. 0) THEN
	        XMIN = MIN(XMIN,XOUT)
	        XMAX = MAX(XMAX,XOUT)
	        goto 30
	      ENDIF
	    ENDDO
	  ENDIF
	  GOTO 30

21	  type *,'Unable to plot file #',fn
	  Ifil(I) = 0
	  goto 30
23	  type *,'Unexpected conversion error for file #',fn
	  Ifil(I) = 0
	  goto 30

30	CONTINUE

	NFIL = 0
	DO i=1,ISP
	  IF (IFIL(I) .GT. 0) THEN
	    NFIL = NFIL + 1
	    IFIL(NFIL) = JFIL(I)
	  ENDIF
	ENDDO
	IF (NFIL .le. 0) THEN
	  Write (*,*) 'ERROR -- There are no remaining files to plot'
	  STATUS = .false.
	  RETURN
	ENDIF

C......................................................................
C set plot limits 
1	IF (XLIMLO) THEN 
	  cxmin1 = XLOW
	ELSEIF (XSCLDF .eq. 0) THEN
	  cxmin1 = xmin
	ELSE
	  QUESTF = .true.
	  IF (QF) then 
	    WRITE (LO,10) 'lower','Q',xmin
	  ELSEIF (DF) then 
	    WRITE (LO,10) 'lower','d-space',xmin
	  ELSEIF (WAVEL .GT. 0) then 
	    WRITE (LO,10) 'lower','two-theta',xmin
	  ELSE
	    WRITE (LO,10) 'lower','x-axis',xmin
	  ENDIF
10	  FORMAT (/'$Enter ',a,' limit for ',a,' [Default=',G10.3']: ')
	  READ (LI,'(Q,f10.0)',end=999,err=1) I,cxmin1
	  IF (I.EQ.0) cxmin1=xmin
	ENDIF
	IF (XLIMHI) THEN 
	  cxmax1 = XHIGH
	ELSEIF (XSCLDF .eq. 0) THEN
	  cxmax1 = xmax
	ELSE
	  QUESTF = .true.
	  IF (QF) then 
	    WRITE (LO,10) 'upper','Q',Xmax
	  ELSEIF (DF) then 
	    WRITE (LO,10) 'upper','d-space',Xmax
	  ELSEIF (WAVEL .GT. 0) then 
	    WRITE (LO,10) 'upper','two-theta',Xmax
	  ELSE
	    WRITE (LO,10) 'upper','x-axis',Xmax
	  ENDIF
	  READ (LI,'(Q,f10.0)',end=999,err=1) I,cxmax1
	  IF (I.EQ.0) cxmax1=xmax
	ENDIF
	IF (cxmin1 .eq. cxmax1) then
	  write(lo,'(2(a,f15.5))') ' Illegal x-axis range ',
	1	cxmin1,' to',cxmax1
	  XLIMLO = .false.
	  XLIMHI = .false.
	  IF (prompt .and. XSCLDF .eq. 3) goto 1
	  STATUS = .false.
	  return
	endif
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C loop for each plot -- get y-axis plot limits
	IF (INDIVF) THEN
C don't need limits here
	  J = 0
	ELSEIF (FIRSTF) THEN
C need limits for 1st file only
	  J = 1
	ELSE
C need limits for all files
	  j = NFIL
	ENDIF
C Find max & min intensity point in range of data to be plotted (CYMX, CYMN) 
	FLAG = .true.
	DO I = 1,J
C convert the file
	  fn=IFIL(I) 
	  IF (PLTYP(fn) .eq. 0) THEN
	    NPTS = 0
	  ELSEIF (QF) THEN 
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_Q,npts,0.0,.false.,*123)
	  ELSEIF (DF) THEN 
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_D,npts,0.0,.false.,*123)
	  ELSEIF (WAVEL .GT. 0) THEN
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_TT,npts,wavel,.false.,*123)
	  ELSE
	    NPTS = NPOINT(FN)
	    DO j=1,NPOINT(FN)
	      A(j)=B(j,FN)
	      AX(j)=BX(j,FN)
	    ENDDO
	  ENDIF
C search through the file of data points
	  do j=1,NPTS
	    if (AX(j).ge.cxmin1 .and. AX(j).le.cxmax1 .and. FLAG) THEN
	      FLAG = .false.
	      CYMX = A(J)
	      CYMN = A(J)
	    elseif (AX(j) .ge. cxmin1 .and. AX(j) .le. cxmax1) THEN
	      CYMX=max(A(J),CYMX)
	      CYMN=min(A(J),CYMN)
	    endif
	  enddo
C If peaks will be plotted, use them as well
	  IF (PPEAKS(fn) .GT. 0) THEN
	    DO J=1,NPKS(fn)
	      IF (QF) THEN
	        call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_Q,npout,0.0,.false.,*23)
	      ELSEIF (DF) THEN
	        npout = 1
	        xout = DSPACE(J,fn)
	      ELSEIF (WAVEL .GT. 0) THEN
	        call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_TT,npout,wavel,.false.,*23)
	      ELSE
	        npout = 1
	        xout = TT(J,fn)
	      ENDIF
	      IF (npout.gt.0 .and. xout.ge.cxmin1 .and. xout.le.cxmax1) then
	        if (FLAG) THEN
	          FLAG = .false.
		  CYMX= pkht(j,fn)
		ELSE
		  CYMX=max(pkht(j,fn),CYMX)
		ENDIF
	      ENDIF
	    ENDDO
	  ENDIF
	ENDDO
C   If /YLOW has been specified, this overrides all other scaling
	IF (YLIMLO) THEN
	  CYMIN = YLOW
	ELSEIF (FIRSTF) THEN
C for default plot from 0.0
	  CYMIN = 0.0
	ELSEIF (YSCLDF .eq. 0 .and. LOGF) THEN
	  CYMIN = cymn
C take LOG of points
	ELSEIF (YSCLDF .eq. 0) THEN
C for default plot from 0
	  CYMIN = 0.0
	ELSE
	  QUESTF = .true.
	  WRITE (LO,10) 'lower','y-axis',CYMN
	  READ (LI,'(Q,f10.0)',end=999,err=1) I,CYMIN
	  IF (I.EQ.0) CYMIN = CYMN
	ENDIF
C   If /YHIGH has been specified, this overrides all other scaling
	IF (YLIMHI) then
	  CYMAX = YHIGH
	ELSEIF (YSCLDF .eq. 0) THEN
	  CYMAX = cymx
	ELSE
	  QUESTF = .true.
	  WRITE (LO,10) 'upper','y-axis',cymx
	  READ (LI,'(Q,f10.0)',end=999,err=1) I,CYMAX
	  IF (I.EQ.0) CYMAX = CYMX
	ENDIF
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C	DETERMINE THE NUMBER OF GRAPHS TO BE PLOTTED AND THE GENERAL
C	NATURE OF THE PLOTS.
C
78	IF (SIZEF .and. (xlen .le. 0 .or. ylen .le. 0)) THEN
	  IF (.not. PROMPT) THEN
	    SIZEF = .FALSE.
	    GOTO 78
	  ENDIF
	  DO WHILE (XLEN .LE. 0.0)
	    WRITE (LO,'(A)') '$Enter X-axis length (inches) [no default]: '
	    READ (LI,FMT='(f10.0)',end=999,err=78) XLEN 
	  ENDDO
	  DO WHILE (YLEN .LE. 0.0)
	    WRITE (LO,'(A)') '$Enter Y-axis length (inches) [no default]: '
	    READ (LI,FMT='(f10.0)',end=999,err=78) YLEN 
	  ENDDO
	  WRITE (LO,'(A)') '$Enter X margin length [default=1 inch]: '
	  READ (LI,FMT='(q,f10.0)',end=999,err=78) ln,XMAR
	  IF (LN.EQ.0) XMAR = 1.
	  WRITE (LO,'(A)') '$Enter Y margin length [default=1 inch]: '
	  READ (LI,FMT='(q,f10.0)',end=999,err=78) ln,YMAR
	  IF (ln.EQ.0) YMAR = 1.
	ENDIF
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20	IF (NFIL.EQ.1) THEN
	  onebox = .true.
	ELSEIF (STKSET .and. STKF) THEN
	  onebox = .false.
	ELSEIF (STKSET .and. .not. STKF) THEN
	  onebox = .true.
	ELSEIF (NFIL.NE.1) THEN
	  QUESTF = .true.
C	  IF (.not. sizef) THEN 
3	  WRITE(LO,'(/a,i2,a,i2,a)') ' There are',NFIL,
	1	' files to be plotted. The files can be plotted in 1 or ',
	1	NFIL,' boxes.',
	3	'$Plot the files in 1 box? (Y/[N]): '
	  READ (LI,'(A1)',end=999) IANS
	  ONEBOX = .FALSE.
	  IF (ians .eq. 'Y' .or. ians .eq. 'y') ONEBOX = .TRUE.
	  IF (ONEBOX) then
	    WRITE (LO,76) NFIL
76	    FORMAT(/' Plotting ',i2,' files on a single set of axes:')
	    WRITE (LO,77) 'y',YSTAKDF
77	    FORMAT(' Choose the amount to offset consecutive plots as a '
	3	'percentage of the ',a,'-axis',/,
	1	' full scale, e.g. 0 to superimpose plots;',
	4	' 100 for no overlap.',/,
	1	'$Enter the offset in % [Default=',i2,']: ')
	    READ (LI,'(q,I5)',end=999,err=3) I,IYSTACK
	    IF (I .eq. 0) iystack = YSTAKDF
	    WRITE (LO,77) 'x',XSTAKDF
	    READ (LI,'(q,I5)',end=999,err=3) I,IXSTACK
	    IF (I .eq. 0) IXSTACK = XSTAKDF
	  endif
	endif
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C set up physical plot size
	chars = 1.3
	charss = chars/40.
C reduce character size if more than one box
	IF (.not. onebox .and. nfil .gt. 1) chars = 1.3 / (1. + nfil/6.)
	call pgsetc(chars)
C set plot size
	IF (sizef) then
C   external scaling
	  ybot = ymar
	  ymarb = .05
	  ymart = .05
	  IF (ONEBOX) THEN
	    ysize = ylen
	  ELSE
	    ysize = ylen / nfil
	  ENDIF
	  XMAR1 = XMAR
	  XLEN1 = XLEN
	ELSE
C internal scaling
C   overall plot borders
	  YTOP = 0.
	  ybot = 3.*charss
C   individual plot borders
	  ymart = 1.*charss
	  ymarb = 0.5*charss
	  IF (ONEBOX) THEN
	    ysize = (1. - ytop - ybot)
	  ELSE
	    ysize = (1. - ytop - ybot) / nfil
	  ENDIF
	  Xleft = 4.5*charss
	  Xright = 0.5*charss
	  XMAR1 = XLEFT
	  XLEN1 = 1. - XRIGHT - XMAR1
	ENDIF
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C plot labels
        IF (QF) THEN 
          lblx = Pxlabl(xtyp_Q) 
        ELSEIF (DF) THEN 
          lblx = Pxlabl(xtyp_D) 
        ELSEIF (WAVEL .GT. 0) THEN 
          lblx = Pxlabl(xtyp_TT) 
	   ln = lench(lblx)
          WRITE(lblx(ln+1:),'(a,f7.4,a)') ' (\gl=',wavel,'\A)'
        else 
	  lblx = Pxlabl(mod(fxtyp(ifil(1)),64))
        endif
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Open graphics device for hardcopy, if requested.
121	IF (HARDF) then
	  CALL OPNHD(lo,li)
	ELSE
	  CALL BGRAF(ICRT,lo,li)
	ENDIF
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	IF (ONEBOX) GOTO 300
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C each plot will be in a separate box
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
200 	DO I=1,NFIL
C++  NFIL is the number of diffractograms to be plotted  
C++  I is the number of the current diffractogram of the NFIL to be plotted 
C++  INDIVF means each file will be plotted to full scale
	  fn=IFIL(I) 
C++ convert the data points to the correct units
	  IF (PLTYP(fn) .eq. 0) THEN
	    NPTS = 0
	  ELSEIF (QF) THEN 
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_Q,npts,0.0,.false.,*123)
	  ELSEIF (DF) THEN 
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_D,npts,0.0,.false.,*123)
	  ELSEIF (WAVEL .GT. 0) THEN
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_TT,npts,wavel,.false.,*123)
	  ELSE
	    NPTS = NPOINT(FN)
	    DO j=1,NPOINT(FN)
	      A(j)=B(j,FN)
	      AX(j)=BX(j,FN)
	    ENDDO
	  ENDIF

C++ Find max & min intensity point in range of data to be plotted (CYMX, CYMN) 
	  IF (INDIVF .and. .not. YLIMHI) THEN
	    FLAG = .true.
	    do j=1,NPTS
	      if (AX(j).ge.cxmin1 .and. AX(j).le.cxmax1 .and. FLAG) THEN
	        FLAG = .false.
	        CYMX = A(J)
	        CYMN = A(J)
	      elseif (AX(j) .ge. cxmin1 .and. AX(j) .le. cxmax1) THEN
	        CYMX=max(A(J),CYMX)
	        CYMN=min(A(J),CYMN)
	      endif
	    enddo
C++ If peaks will be plotted, convert them too
	    IF (PPEAKS(fn) .GT. 0) THEN
	      DO J=1,NPKS(fn)
	        IF (QF) THEN
	          call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_Q,npout,0.0,.false.,*23)
	        ELSEIF (DF) THEN
	          npout = 1
	          xout = DSPACE(J,fn)
	        ELSEIF (WAVEL .GT. 0) THEN
	          call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_TT,npout,wavel,.false.,*23)
	        ELSE
	          npout = 1
	          xout = TT(J,fn)
	        ENDIF
		IF (npout.gt.0 .and. xout.ge.cxmin1 .and. xout.le.cxmax1) THEN
		  IF (FLAG) THEN
	            FLAG = .false.
		    CYMX = pkht(j,fn)
		  ELSE
		    CYMX=max(pkht(j,fn),CYMX)
		  ENDIF
		ENDIF
	      ENDDO
	    ENDIF
	    CYMIN = 0.0
	    CYMAX = CYMX
	  ENDIF
C--------------------------------------------------------------------------
C++ set up dimensions of box for current plot
	  ylw = ybot + ymarb + (i-1.)*ysize
	  yhgh = ybot - ymart + i*ysize
	  IF (sizef) then 
	    CALL PGVSIZE(XMAR1, XMAR1+XLEN1, ylw, yhgh)
	  ELSE
	    CALL PGVPORT(xmar1, xmar1+xlen1, ylw, yhgh)
	  ENDIF
C++    Label with the name/title of just the current file
	  call pgsetc(chars*0.9)
	  call PGSCI(PCOLOR(fn))
	  IF (FILELBL) call pgmtext('T',-1.5,0.98,1.0,cfilnam(fn))
	  IF (TITLEF) call pgmtext('T',-1.5,0.02,0.0,ATITLE(fn)(:60))
	  call pgsetc(chars)
	  call PGSCI(1)

C adjust the scaling by the plot offset/scale
	  IF (PLMAG(fn) .ne. 0) THEN
	    CYMAX1 = CYMAX1 / PLMAG(fn)
	    CYMIN1 = CYMIN1 / PLMAG(fn)
	  ENDIF
	  CYMAX1 = CYMAX - PLOFF(fn)
	  CYMIN1 = CYMIN - PLOFF(fn)
C++ take SQRT, if appropriate
	  IF (SQRTF) then
	    CYMIN1 = SQRT(MAX(CYMIN1,0.))
	    CYMAX1 = SQRT(MAX(CYMAX1,0.))
	    do j=1,NPTS
	      A(J)=SQRT(MAX(A(J),0.))
	    enddo
C++ take LOG, if appropriate
	  ELSEIF (LOGF) THEN
	    IF (CYMAX1 .le. 0) THEN
	      write(lo,'(A)') ' ERROR -- Plot max. is negative'
	      STATUS = .false.
	      RETURN
	    ENDIF
	    IF (CYMIN1 .le. 0.) CYMIN1 = CYMAX1 / 1000.
	    do j=1,NPTS
	      A(J) = LOG10(MAX(A(J),CYMIN1))
	    enddo
	    CYMAX1 = LOG10(CYMAX1)
	    CYMIN1 = LOG10(CYMIN1)
	  ENDIF
C++     make sure that YMAX <> YMIN
	  IF (CYMAX1 .EQ. CYMIN1) goto 92
	  CALL PGWINDOW(cxmin1,cxmax1,cymin1,cymax1)
	  IF (SQRTF) then
	    CALL PGBOX('bsit',0.0,0,'bstinv',0.0,0)
	  ELSEIF (LOGF) THEN
	    CALL PGBOX('bsit',0.0,0,'bstlinv',0.0,0)
	  ELSE
	    CALL PGBOX('bsit',0.0,0,'bstinv',0.0,0)
	  ENDIF
	  CALL PGBOX('cist',0.0,0,'ct',0.0,0)
	  if (i .eq. 1) CALL PGBOX('n',0.0,0,' ',0.0,0)
	  if (i .eq. 1) CALL PGLABEL(lbLX,' ',' ')
	  if (GRIDF) CALL PGBOX('g',0.0,0,'g',0.0,0)
	  LBLY = pylabl(mod(fytyp(fn),64))
	  IF (SQRTF) LBLY = '\(2267)'//lbly
	  call pgmtext('L',4.0,0.5,0.5,lblY)
C--------------------------------------------------------------------------
C++ set the color for the file
	  call PGSCI(PCOLOR(fn))
C++   PLOT THE DATA   (IF THE line type IS GREATER THAN ZERO) 
	  IF (PLTYP(fn) .GT. 0) THEN
C++     scale Y data; clip to CYMAX1 if necessary
	    if (clip) then 
	      do j=1,NPTS
	        a(j) = min(cymax1,a(j))
	        a(j) = max(cymin1,a(j))
	      enddo
	    endif
C++     set the width and line type for the diffractogram
	    IF (PWIDTH(fn) .GT. 0) CALL PGSLW(PWIDTH(fn))
	    IF (PLTYP(fn) .GE. 6) THEN
	      call PGPOINT(npts,AX,a,-1)
	    ELSE
	      CALL PGSLS(PLTYP(fn))
	      call PGLINE(npts,AX,a)
	    ENDIF
	  ENDIF
C++  PLOT THE PEAKS (IF THE PEAK LINE WIDTH IS GREATER THAN ZERO)
 	  IF (MOD(PPEAKS(fn),8) .GT. 0) THEN
 	    CALL PGSLS(MOD(PPEAKS(fn),8))
	    DO II=1,NPKS(fn)
	      XT(1) = TT(II,fn)
	      XT(2) = TT(II,fn)
C++   if PPEAKS > 8, intensity scale is arbitrary
	      IF (PPEAKS(fn) .gt. 8) THEN
	        YT(1) = CYMIN1
	        YT(2) = CYMAX1 
	      ELSE
	        YT(1) = CYMIN1
	        YT(2) = YT(1) + (PKHT(II,fn) + BKG(II,fn))
	      ENDIF
	      IF (QF) THEN
	        call XYCONV(XT,YT,FXTYP(fn),2,WAVE(fn),
	1		    XT,YT,XTYP_Q,npout,0.0,.false.,*123)
	      ELSEIF (DF) THEN
	        call XYCONV(XT,YT,FXTYP(fn),2,WAVE(fn),
	1		    XT,YT,XTYP_D,npout,0.0,.false.,*123)
	      ELSEIF (WAVEL .GT. 0) THEN
	        call XYCONV(XT,YT,FXTYP(fn),2,WAVE(fn),
	1		    XT,YT,XTYP_TT,npout,wavel,.false.,*123)
	      ELSE
	        npout = 2
	      ENDIF
	      IF (NPOUT .eq. 2) THEN
	        call PGLINE(2,XT,YT)
C++   draw a tick mark at the background position
	        IF (BKG(II,fn) .gt. 0) THEN
	          YT(1) = YT(1) + BKG(II,fn)
	          YT(2) = YT(1)
	          XT(1) = XT(2) - (cxmax1 - cxmin1)/200.
	          XT(2) = XT(2) + (cxmax1 - cxmin1)/200.
	          call PGLINE(2,XT,YT)
	        ENDIF
	      ENDIF
	    ENDDO
	  ENDIF

	  CALL PGSLS(1)
	  call PGSCI(1)
	  CALL PGSLW(1)
	ENDDO
	GOTO 400
C======================================================================
C all plots plot will be in one box (ONEBOX = .TRUE.)
C======================================================================
C== set up dimensions for current plot
300	ymar1 = ybot + ymarb
	ylen1 = ysize - ymart
C== set box size for current plot
	ylen2 = ylen1 * 100. / (100. + (NFIL-1) * IYSTACK)
	Xlen2 = xlen1 * 100. / (100. + (NFIL-1) * IXSTACK)
 	DO I=1,NFIL
C==  NFIL is the number of diffractograms to be plotted  
C==  I is the number of the current diffractogram of the NFIL to be plotted 
C==  INDIVF means each file will be plotted to full scale, plot will be labeled
C==         from 0 to 100(+)%
	  fn = IFIL(I) 
C== convert the data points to the correct units
	  IF (PLTYP(fn) .eq. 0) THEN
	    NPTS = 0
	  ELSEIF (QF) THEN 
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_Q,npts,0.0,.false.,*123)
	  ELSEIF (DF) THEN 
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_D,npts,0.0,.false.,*123)
	  ELSEIF (WAVEL .GT. 0) THEN
	    call XYCONV(BX(1,fn),B(1,fn),FXTYP(fn),NPOINT(fn),WAVE(fn),
	1		AX,A,XTYP_TT,npts,wavel,.false.,*123)
	  ELSE
	    NPTS = NPOINT(FN)
	    DO j=1,NPOINT(FN)
	      A(j)=B(j,FN)
	      AX(j)=BX(j,FN)
	    ENDDO
	  ENDIF

C== Find max & min intensity point in range of data to be plotted (CYMX, CYMN) 
	  IF (INDIVF) THEN
	    FLAG = .true.
	    do j=1,NPTS
	      if (AX(j).ge.cxmin1 .and. AX(j).le.cxmax1 .and. FLAG) THEN
	        FLAG = .false.
	        CYMX = A(J)
	        CYMN = A(J)
	      elseif (AX(j) .ge. cxmin1 .and. AX(j) .le. cxmax1) THEN
	        CYMX=max(A(J),CYMX)
	        CYMN=min(A(J),CYMN)
	      endif
	    enddo
C== If peaks will be plotted, convert them too
	    IF (PPEAKS(fn) .GT. 0) THEN
	      DO J=1,NPKS(fn)
	        IF (QF) THEN
	          call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_Q,npout,0.0,.false.,*23)
	        ELSEIF (DF) THEN
	          npout = 1
	          xout = DSPACE(J,fn)
	        ELSEIF (WAVEL .GT. 0) THEN
	          call XYCONV(dspace(J,fn),0.0,xtyp_D,1,0.0,
	1		   XOUT,YOUT,XTYP_TT,npout,wavel,.false.,*23)
	        ELSE
	          npout = 1
	          xout = TT(J,fn)
	        ENDIF
		IF (npout.gt.0 .and. xout.ge.cxmin1 .and. xout.le.cxmax1) THEN
		  IF (FLAG) THEN
	            FLAG = .false.
		    CYMX = pkht(j,fn)
		  ELSE
		    CYMX=max(pkht(j,fn),CYMX)
		  ENDIF
		ENDIF
	      ENDDO
	    ENDIF
C CYMAX is the full scale intensity
	    CYMIN1 = 0.
	    CYMAX1 = CYMX
	  ELSE
	    CYMIN1 = CYMIN
	    CYMAX1 = CYMAX
	  ENDIF
C----------------------------------------------------------------------
C--------------------------------------------------------------------------
C adjust the scaling by the plot offset/scale
	  IF (PLMAG(fn) .ne. 0) THEN
	    CYMAX1 = CYMAX1 / PLMAG(fn)
	    CYMIN1 = CYMIN1 / PLMAG(fn)
	  ENDIF
	  CYMAX1 = CYMAX1 - PLOFF(fn)
	  CYMIN1 = CYMIN1 - PLOFF(fn)
C== take SQRT, if appropriate
	  IF (SQRTF) then
	    CYMIN1 = SQRT(MAX(CYMIN1,0.))
	    CYMAX1 = SQRT(MAX(CYMAX1,0.))
	    do j=1,NPTS
	      A(J) = SQRT(MAX(A(J),0.))
	    enddo
C== take LOG, if appropriate
	  ELSEIF (LOGF) THEN
	    IF (CYMAX1 .le. 0) THEN
	      write(lo,'(A)') ' ERROR -- Plot max. is negative'
	      STATUS = .false.
	      RETURN
	    ENDIF
	    IF (CYMIN1 .le. 0.) CYMIN1 = CYMAX1 / 1000.
	    do j=1,NPTS
	      A(J) = LOG10(MAX(A(J),CYMIN1))
	    enddo
	    CYMAX1 = LOG10(CYMAX1)
	    CYMIN1 = LOG10(CYMIN1)
	  ENDIF
C== make sure that CYMAX1 <> CYMIN1
	  IF (CYMAX1 .EQ. CYMIN1) goto 92

	  call pgsetc(chars)
	  call PGSCI(1)
C----------------------------------------------------------------------
C==    Is this the first plot?
	  if (i .eq. 1) then
C== set the box we will use for the Nth plot
	    IF (sizef) then 
	      CALL PGVSIZE(xmar1, xmar1+xlen1, ymar1, ymar1+ylen1)
	    ELSE
	      CALL PGVPORT(xmar1, xmar1+xlen1, ymar1, ymar1+ylen1)
	    ENDIF
C==    Label here with all names on the 1st pass
	    call pgsetc(chars*0.9)
	    do ii=1,NFIL
	      iif = ifil(ii)
	      call PGSCI(PCOLOR(IFIL(ii)))
	      IF (FILELBL) call pgmtext('T',1.5*(II-NFIL)-2.0,0.98,1.0,
	1		cfilnam(iif))
	      IF (TITLEF) call pgmtext('T',1.5*(II-NFIL)-2.0,0.02,0.0
	1		,ATITLE(iif)(:60))
	    enddo
	    call PGSCI(1)
	    call pgsetc(chars)
C==   now label the axes
C==     set the full scale y axis values
	    IF (INDIVF) THEN
	      cymin2 = 0.
	      cymax2 = 100.
	    ELSE
	      cymin2 = cymin1
	      cymax2 = cymax1
	    ENDIF
C== increase the axis length by the expansion factor for stacking
	    cxmin2 = cxmin1
	    cxmax2 = cxmax1
	    cxmax2 = cxmin2 + (cxmax2 - cxmin2) * (xlen1 / xlen2)
	    cymax2 = cymin2 + (cymax2 - cymin2) * (ylen1 / ylen2)
C== label the axies & box
	    CALL PGWINDOW(cxmin2, cxmax2, cymin2, cymax2)
C for stacked plots compute a x-tick mark increment
	    i1 = 0
	    i2 = 0
	    IF (IYSTACK .gt. 0 .and. IXSTACK .GT. 0) THEN
	      CALL calctick(cxmin2, cxmax2, xinc1, nsubx, i1, i2)
	    else
	      xinc1 = 0.0
	      nsubx = 0	     
	    endif
	    CALL PGBOX('cst', xinc1, nsubx, 'ct',0.0,0)
	    CALL PGBOX('n', xinc1, nsubx, ' ',0.0,0)
	    CALL PGLABEL(lbLX,' ',' ')
	    if (GRIDF) CALL PGBOX('g',0.0,0,'g',0.0,0)
	    LBLY = pylabl(mod(fytyp(fn),64))
	    IF (INDIVF) LBLY = 'Relative Intensity'
	    IF (SQRTF) LBLY = '\(2267)'//lbly
	    call pgmtext('L',4.,0.5,0.5,lblY)
	    IF (SQRTF) then
	      CALL PGBOX('bsit', xinc1, nsubx, 'bstinv',0.0,0)
	    ELSEIF (LOGF) THEN
	      CALL PGBOX('bsit', xinc1, nsubx, 'bstlinv',0.0,0)
	    ELSE
	      CALL PGBOX('bsit', xinc1, nsubx, 'bstinv',0.0,0)
	    ENDIF
C draw a line along the "Z" axis
	    IF (IYSTACK .gt. 0 .and. IXSTACK .GT. 0) THEN
	      CALL PGWINDOW(xmar1, xmar1+xlen1, ymar1, ymar1+ylen1)
C== set margin size for last plot
	      ymar2 = ymar1 + ylen2 * (NFIL-1) * IYSTACK/100.
	      xmar2 = xmar1 + xlen2 * (NFIL-1) * IXSTACK/100.
	      CALL PGMOVE(xmar1, ymar1)
	      CALL PGDRAW(xmar2, ymar2)
	      do i3=i1,i2
        	x1 = xlen1 * (I3*xinc1 - cxmin2)/(cxmax2-cxmin2)
	        CALL PGMOVE(xmar1 + x1, ymar1)
	        CALL PGDRAW(xmar2 + x1, ymar2)
	      enddo
	      do i3=1,nfil-1
	        ymar2 = ymar1 + ylen2 * I3 * IYSTACK/100.
	        xmar2 = xmar1 + xlen2 * I3 * IXSTACK/100.
	        CALL PGMOVE(xmar2, ymar2)
	        CALL PGDRAW(xmar2 + xlen1, ymar2)
	      enddo
	    ENDIF
	  endif
C----------------------------------------------------------------------
C== set margin size for current plot
	  ymar2 = ymar1 + ylen2 * (I-1) * IYSTACK/100.
	  Xmar2 = xmar1 + xlen2 * (I-1) * IXSTACK/100.
C== set the box we will use for the Nth plot
	  IF (sizef) then 
	    CALL PGVSIZE(xmar2, xmar2+xlen2, ymar2, ymar2+ylen2)
	  ELSE
	    CALL PGVPORT(xmar2, xmar2+xlen2, ymar2, ymar2+ylen2)
	  ENDIF
C== scale the box we will use for the Nth plot
	  CALL PGWINDOW(cxmin1, cxmax1, cymin1, cymax1)
C--------------------------------------------------------------------------
C set the color for the file
	  call PGSCI(PCOLOR(fn))
C   PLOT THE DATA   (IF THE line type IS GREATER THAN ZERO) 
	  IF (PLTYP(fn) .GT. 0) THEN
C   scale Y data; clip to CYMAX1 if necessary
	    if (clip) then 
	      do j=1,NPTS
	        a(j) = min(cymax1,a(j))
	        a(j) = max(cymin1,a(j))
	      enddo
	    endif
C   set the width and line type for the diffractogram
	    IF (PWIDTH(fn) .GT. 0) CALL PGSLW(PWIDTH(fn))
	    IF (PLTYP(fn) .GE. 6) THEN
	      call PGPOINT(npts,AX,a,-1)
	    ELSE
	      CALL PGSLS(PLTYP(fn))
	      call PGLINE(npts,AX,a)
	    ENDIF
	  ENDIF
C----------------------------------------------------------------------
C PLOT THE PEAKS (IF THE PEAK LINE WIDTH IS GREATER THAN ZERO)
 	  IF (MOD(PPEAKS(fn),8) .GT. 0) THEN
 	    CALL PGSLS(MOD(PPEAKS(fn),8))
	    flag = .true.
	    xt1 = CYMIN1
	    xt2 = CYMIN1
	    DO II=1,NPKS(fn)
	      XT(1) = TT(II,fn)
	      XT(2) = TT(II,fn)
C if PPEAKS > 8, intensity scale is arbitrary
	      IF (PPEAKS(fn) .gt. 8) THEN
	        YT(1) = CYMIN1
	        YT(2) = CYMAX1 
	      ELSE
	        YT(1) = CYMIN1
	        YT(2) = YT(1) + (PKHT(II,fn) + BKG(II,fn))
	      ENDIF
	      IF (QF) THEN
	        call XYCONV(XT,YT,FXTYP(fn),2,WAVE(fn),
	1		    XT,YT,XTYP_Q,npout,0.0,.false.,*123)
	      ELSEIF (DF) THEN
	        call XYCONV(XT,YT,FXTYP(fn),2,WAVE(fn),
	1		    XT,YT,XTYP_D,npout,0.0,.false.,*123)
	      ELSEIF (WAVEL .GT. 0) THEN
	        call XYCONV(XT,YT,FXTYP(fn),2,WAVE(fn),
	1		    XT,YT,XTYP_TT,npout,wavel,.false.,*123)
	      ELSE
	        npout = 2
	      ENDIF
	      IF (NPOUT .eq. 2) THEN
	        call PGLINE(2,XT,YT)
C save the 1st and last peak position
	        IF (flag) then
	          flag = .false.
		  xt1 = XT(1)
		  xt2 = XT(1)
	        ELSE
		  xt1 = min(XT(1),xt1)
		  xt2 = max(XT(1),xt2)
	        ENDIF
C draw a tick mark at the background position
	        IF (BKG(II,fn) .gt. 0) THEN
	          YT(1) = YT(1) + BKG(II,fn)
	          YT(2) = YT(1)
	          XT(1) = XT(2) - (cxmax1 - cxmin1)/200.
	          XT(2) = XT(2) + (cxmax1 - cxmin1)/200.
	          call PGLINE(2,XT,YT)
	        ENDIF
	      ENDIF
	    ENDDO
C draw a line connecting the peaks at the baseline
	    IF (XT1 .ne. XT2) THEN
	      YT(1) = 0.0
	      YT(2) = 0.0
	      XT(1) = XT1
	      XT(2) = XT2
	      call PGLINE(2,XT,YT)
	    ENDIF
	  ENDIF

	  CALL PGSLS(1)
	  call PGSCI(1)
	  CALL PGSLW(1)
	ENDDO
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C close hardcopy graphics device and reopen terminal
400	IF (HARDF) then
	  CALL CLHD(lo,li)
	  hardf = .false.
	ELSE
	  WRITE(LO,'(A)')
	1	'$DO YOU WANT TO REPEAT THIS PLOT ON HARDCOPY? (Y/[N]): '
	  IANS = 'n'
	  READ (LI,'(A1)',end=91) IANS
C close the graphics terminal 
91	  CALL EGRAF(ICRT,lo,li)
	  IF(IANS.EQ.'Y' .or. IANS.EQ.'y') THEN 
	    hardf = .true.
	    goto 121
	  ENDIF
	ENDIF

	IF (QUESTF) THEN 
	  WRITE(LO,'(/A)')
	1 '$Do you want to change your answers to the previous '//
	1	'questions? (Y/[N]): '
	  READ (LI,'(A1)',end=999) IANS
	  IF(IANS.EQ.'Y' .or. IANS.EQ.'y') goto 1
	ENDIF

	STATUS = .true.
999	RETURN
92	CALL EGRAF(ICRT,lo,li)
	write (lo,'(A)') ' *** Scaling error: Ymax and Ymin are equal'
	IF (QUESTF) goto 1
	STATUS = .false.
	RETURN
123	type *,'Unexpected conversion error for file #',fn
	STATUS = .false.
	RETURN
	END
C----------------------------------------------------------------------
	subroutine calctick(XBLC, XTRC, xint, nsubx, i1, i2)
      INTEGER  I,ILOG
      REAL     FRAC,NICE(3),PWR,XLOG,XX
      INTRINSIC ABS, LOG10, SIGN
      DATA     NICE/2.0,5.0,10.0/

C Choose X tick intervals. Major interval = XINT,
C minor interval = XINT/NSUBX.
C
!      XINT = MAX(0.05, MIN(7.0*XSP/XLEN, 0.20))*(XTRC - XBLC)
      XINT = .1 * (XTRC - XBLC)

C
	I1 = 0
	I2 = 0
      IF (XINT.EQ.0.0) THEN
          XINT = 0.0
          NSUBX = 2
          RETURN
      END IF
      XX   = ABS(XINT)
      XLOG = LOG10(XX)
      ILOG = XLOG
      IF (XLOG.LT.0) ILOG=ILOG-1
      PWR  = 10.0**ILOG
      FRAC = XX/PWR
      DO 10 I=1,3
          IF (FRAC.LE.NICE(I)) GOTO 20
   10 CONTINUE
   20 XINT = SIGN(PWR*NICE(I),XINT)
      NSUBX = 5
      IF (NICE(I).EQ.2.0) NSUBX = 2

      NP = INT(LOG10(ABS(XINT)))-4
      NV = NINT(XINT/10.**NP)

      XLO = MIN(XBLC/XINT, XTRC/XINT)
      XHI = MAX(XBLC/XINT, XTRC/XINT)
      I1 = NINT(XLO)
      IF (I1.LT.XLO) I1 = I1+1
      I2 = NINT(XHI)
      IF (I2.GT.XHI) I2 = I2-1

	RETURN
	END
