C file psgraph1.for
C ----------------------------------------------------------------------!-------
C includes PSGRAPH, USER_PROMPT, ADV_FEATURES, GET_DEFAULTS, INIT_PS_FILE, 
C FRAME, DRAW_TICKS_PS, SIGNIF_DIGITS, LABEL_AXES_PS, TRIM_STRING, 
C LEN_TRIMMED, CNVRT_LOWER, PLOT_PTS_PS, RAYPLOT_PTS_PS, RAYTXT_OUT_PS, 
C SUBHEAD_PS 
C ----------------------------------------------------------------------!-------
! 15 Aug 2011 - Added spaces after "cir", "cirs", and "cirfs" as recommended by
!               Larry Baker
C 23 Jun 2009 - Changed the default character size for DDTEXT to 1.
C 23 Jun 2009 - Increased size of text for DDTEXT for isize=2 and 3
C 23 Jun 2009 - Decreased size of text for DDTEXT for isize=0
C 22 Jun 2009 - Read in intensity, italics, and size of text
C 22 Aug 2008 - The string can be blank filled without a Do loop 
C               (complements of Larry Baker)
C 22 Mar 2008 - Fix READ raytxt to pad with blanks
C	Use full character variable assignment instead of loop to fill
C	raytxt with blanks
C 20 Mar 2006 - Generalized text printing
C 2 Sept 2004 - Replaced more F6.3s with G10.3s
C 27 June 2004 - Replaced F6.3 by G10.3 in 5 places
C ------------------------ PROGRAM PSGRAPH -----------------------------!-------
C
C This program creates a PostScript output file and graphics screen
C output from a TXTOUT input file.  TXTOUT is the output from one of the
C raytracing programs, PROFILE,CONPLT,EIGEN,PC-HARPO.  TXTOUT contains plotting
C commands that are compatible with the DISSPLA graphics software package.
C In this program, PSGRAPH, those commands are converted to their PostScript
C equivalents (if they exist).
C
C The graphics screen output from this program is PC-dependent since it
C uses Microsoft Graphics library routines.
C
C Cannot use IMPLICIT none with FGRAPH.FI.
C
C NOTE:  There are various hooks in PSGRAPH to handle dashed lines.  However,
C        the program that provides input to PSGRAPH e.g., PROFILE, would have
C        to be modified to create command 37 in the TXTOUT file.  Therefore,
C        if none of the raytracing group of programs has been modified to
C        create command 37, none of the dashline code in PSGRAPH will be
C        executed.  (You will find that MicroSoft Graphics do a poor job of
C        dashed lines when drawing line segments of varying length, as is
C        the usual case for the raytracing suite of programs.)

c      INCLUDE 'FGRAPH.FI'

      INTEGER
     &	     fontsize,
     &       iunit,
     &	     ounit,
     &	     ist,
     &	     icode,
     &	     title_len,
     &	     horiz_label_len,
     &	     vert_label_len,
     &	     numpts,
     &	     numpts_ray,
     &       numpts_con,
     &	     nremain,
     &	     marker,
     &	     i,
     &	     numxticks,
     &	     numyticks,
     &       num_pt_plots,
     &       lennew

      INTEGER*2
     &	     xtmp(0:4000), ytmp(0:4000),
     &       numtextrows


      REAL
     &	     scale,
     &	     xorig,
     &	     yorig,
     &	     x_axis_len,
     &	     y_axis_len,
     &	     xmin, ymin,
     &	     xinc, yinc,
     &	     xmax, ymax,
     &	     x(1000), y(1000),
     &       rlinewidth,
     &       deflwidth,
     &       slomofact


      PARAMETER (deflwidth = .01)

      CHARACTER*80 title,
     &		   horiz_label,
     &		   vert_label,
     &		   newtitle,
     &		   new_horiz_label,
     &		   new_vert_label

      CHARACTER*80 fname, outfname,
     &		   raytxt,
     &             newraytxt,
     &             resetxt

      PARAMETER (iunit = 1, ounit = 2)

      LOGICAL
     &	     start,
     &	     existflag,
     &       newline,
     &       rayflag,
     &       conflag,
     &       eigflag,
     &       dashline,
     &       translated,
     &       subheaddun


C==>
c      print*,'Beginning of psgraph1'

      call USER_PROMPT (fname, scale, fontsize, rlinewidth, edtfactor,
     &                  outfname, slomofact)

      call GET_DEFAULTS (fname, scale, fontsize, rlinewidth, edtfactor,
     &                   outfname, slomofact)

      start = .true.
      newline = .false.
      dashline = .false.
      translated = .false.
      subheaddun = .false.
      num_pt_plots = 0

      INQUIRE (FILE=fname, EXIST=existflag)
      if (.not.existflag) then
	   PRINT*,'ERROR:  Input file does not exist: ',fname
	   STOP 'Check file name and re-start.'
      endif

      OPEN (UNIT=iunit,FILE=fname,FORM='FORMATTED')
      REWIND iunit

C ----------------------------------------------------------------------!-------

C.....Loop until an ICODE for end of plotting is reached or an error occurs.

10    READ (iunit,111,IOSTAT=ist) icode
c      print*,'icode = ',icode

111   FORMAT (I2)

      if (ist.gt.0) then
	 PRINT 8888,ist
8888	 FORMAT(' ERROR: ',I5,' while reading input file.')
	 STOP ' Check error code in MicroSoft FORTRAN Reference Manual.'
      else if (ist.lt.0) then
	 PRINT *,'Unexpected end of input file.'
	 STOP 'ERROR: Input file should end with code 31 or code -1.'
      endif

C NOTE:  Several of the DISSPLA command codes can be ignored.   Those
C        commands are:
C		     10,11,12,13,20,21,23,24,26,32,33,34,35,36

      if (icode .eq. 10 .or. icode .eq. 11 .or. icode .eq. 12 .or.
     &    icode .eq. 13 .or. icode .eq. 20 .or. icode .eq. 21 .or.
     &	  icode .eq. 23 .or. icode .eq. 24 .or. icode .eq. 26 .or.
     &	  (icode .ge. 32 .and. icode .le. 36)) goto 10


      if (start) then
           call INIT_PS_FILE (outfname, scale, fontsize, rlinewidth)
           call GET_VIDEO_MODE (numtextrows)
           start = .false.
      endif

C********** NOTE:  Some icodes apply only to profile plots and some
C		   apply only to rayplots:
C
C		   RAYPLOTS
C		   --------
C
C		      0      Frame start
C		      1      Moveto with linetos
C		      3      Text location + text
C		     -2      Frame advance
C		     -1      End of plotting
C
C		   PROFILES
C		   --------
C
C		     10      Font indicator (vestigial, from DISSPLA)
C		     11      Case indicator (vestigial, from DISSPLA)
C		     12      Case indicator (vestigial, from DISSPLA)
C		     13      Text Height    (vestigial, from DISSPLA)
C		     20      File Compression (vestigial, from DISSPLA)
C		     21      Clipping Indicator (ignored)
C		     22      Lower Left origin in inches from corner
C		     23      Page Limits (ignored)
C		     24      Scale Factor for plot symbols (ignored)
C		     25      Title Text and Axis Labels
C		     26      Draw a frame (ignored)
C		     29      Points to be plotted
C		     30      Frame advance
C		     31      End of plotting
C		     32      Number of X axis ticks (ignored)
C		     33      Number of Y axis ticks (ignored)
C		     34      Axis labeling  (vestigial, from DISSPLA)
C		     35      Axis labeling  (vestigial, from DISSPLA)
C		     36      Suppress border (vestigial, from DISSPLA)
C		     37      Dash lines
C		     38      Reset from dash to solid lines
C		     39      Max, Min and increment in real units
C                    40      Eigen Plot indicator (from EIGEN program,not
C                            from PC-HARPO)

C.....Icode = -1 : End of plotting for rayplots.

      if (icode .eq. -1) then
           WRITE (2,*) 'grestore'
	   call SETVIDEO_DFT

	   STOP 'NORMAL END'

C.....Icode = 0 is a heading for the file.  The heading is ignored by this
C     program, but this is an opportune time to put a translation of the
C     PostScript page into effect since this is the beginning of the processing.
C     Note that this is only done ONCE per plotting session, not once per plot.

      else if (icode .eq. 0) then
           READ (1,*)
           if (.not. translated) then
                WRITE (2,*) '0.4 in 2.0 in translate'
                translated = .true.
           endif

C.....ICODE = 1:  A "moveto" command for rayplots.  This is followed
C     by NUMPTS_RAY, the number of points in the rayplot curve.  This is
C     followed by a set of x,y data points which will comprise the curve to
C     be plotted.
C     There are 8 pairs of data points per line in the input file.  If the
C     number data points is not a multiple of 8, then the remaining pairs are
C     each on a single line.

      else if (icode .eq. 1) then
	   BACKSPACE iunit
	   READ (iunit,77) icode, xtmp(0), ytmp(0)
77	   FORMAT (I2,2I5)
	   READ (iunit,88) numpts_ray
88	   FORMAT (I5)

           nremain = MOD(numpts_ray,8)

C..........Set RAYFLAG here. 3/5/91. Formerly set if NUMPTS > 0.

           rayflag = .true.

C..........Read all the points from TXTOUT

	   if (numpts_ray .ge. 8) then
		do ii = 1, numpts_ray-7, 8
		     READ (iunit,99,ERR=100,IOSTAT=ist)
     &				     (xtmp(j),ytmp(j),j=ii,ii+7)
99		     FORMAT (8(2I5))
		enddo
	   endif
	   if (nremain.gt.0) then
		READ (iunit,66,ERR=100,IOSTAT=ist)
     &		     (xtmp(k),ytmp(k),k=numpts_ray-nremain+1,numpts_ray)
66		FORMAT (2I5)
	   endif

           if (numpts_ray .gt. 0) then

                marker = 1
                conflag = .false.

C...............Note that X and Y are just dummy arguments here.  They
C               are REAL arrays which are only needed for profile plots.

		call PLOT_PTS_VT (1655.,9785.,605.,7515.,xtmp,ytmp,
     &                            x, y,
     &                            numpts_ray, marker, rayflag, conflag,
     &                            numxticks, numyticks, num_pt_plots,
     &                            slomofact, eigflag)

C...............Set the PS linewidth to user-selected value only if
C               we are plotting data, that is, not the frame for the
C               plot.  The frame plotting is made up of several (4)
C               passes and usually takes far less than 50 points per pass.
C               Also, don't bother setting the line width if RLINEWIDTH
C               is .01 since it was set at initialization.

                if (numpts_ray .gt. 50 .and. rlinewidth .ne.
     &              deflwidth) then
                     WRITE (2,4444) rlinewidth
                else if (rlinewidth .ne. deflwidth) then
                     WRITE (2,*) '0.01 in setlinewidth'
                endif
                call RAYPLOT_PTS_PS (numpts_ray, xtmp, ytmp, edtfactor,
     &                               dashline, conflag)

	   endif

	   numpts_ray = 0

C.....ICODE = 3:  Text for rayplots, e.g. titles, axis labels, tick labels.
C                 Also, user-supplied text for profiles in DINP.

      elseif (icode .eq. 3) then

	   BACKSPACE iunit
	   READ (iunit,11) icode, xtxt, ytxt
11	   FORMAT (I2,2F4.0)

C           READ (iunit,22) iorientn, nchar
C22	   FORMAT (2I3)

           READ (iunit,22) iorientn, nchar, int, ital, isize      ! 22 Jun 2009
22	   FORMAT (5I3)                                           ! 22 Jun 2009

C ! 22 Mar 2008	   READ (iunit,33) (raytxt(jj:jj),jj=1,nchar)
C ! 22 Mar 2008 33	   FORMAT (80A1)
           READ(iunit,'(A)') raytxt                                ! 22 Mar 2008

c           if (rayflag) then           ! 20 Mar 2006
c           print*,'icode = ',icode
c           print*,'xtxt,ytxt = ',xtxt,ytxt

C...............Get rid of any # or ! symbols before output to VT or PS.

                call TRIM_STRING (raytxt, newraytxt, lennew)
                call RAYTXT_OUT_VT (1655.,9785.,605.,7515.,iorientn,
     &                              lennew, newraytxt, xtxt, ytxt)
C                call RAYTXT_OUT_PS (newraytxt, xtxt, ytxt, iorientn,
C     &                              lennew)
                call RAYTXT_OUT_PS (newraytxt, xtxt, ytxt, iorientn,
     &                       lennew,fontsize,int,ital,isize)      ! 22 Jun 2009

C ! 22 Mar 2008                 do kkk = 1,80
C ! 22 Mar 2008                    raytxt(kkk:kkk) = ' '
C ! 22 Mar 2008                 enddo
               raytxt=''                                           ! 22 Mar 2008 

C..........Handle the special case of subheading text from Profile
C          and Eigen plots.  This is a text string which was put in the DINP
C          file by the user.  Usually, it contains some additional
C          info about a particular plot e.g. a name for the data set used.
C          This is sent to PS file by LABEL_AXES_PS and to screen by
C          LABEL_AXES_VT.

c           else                                      ! 20 Mar 2006
c               lentxt = LEN_TRIMMED (raytxt)          ! 20 Mar 2006
c           endif                                     ! 20 Mar 2006

	   GOTO 10

C.....ICODE = 4:  Handle contour plot arrays.  Note that contours are
C     composed of a series of "moveto,lineto"s, that is, one moveto followed
C     by one lineto.  Therefore, the rayplot code would not work.  In the
C     contour program (CONPLT), NUMPTS_CON is written to TXTOUT as the number
C     of moveto,lineto pairs because in that program, there are 4 arrays :
C     two "moveto" arrays and two "lineto" arrays (one for x,one for y).  In
C     PSGRAPH, I have just re-used the 2 XTMP,YTMP arrays, each of which will
C     contain lineto, moveto etc.   Therefore, we run the array to
C     2*NUMPTS_CON.

      elseif (icode .eq. 4) then

            conflag = .true.
            READ (iunit,2222,ERR=100,IOSTAT=ist) numpts_con
2222        format (I4)

C..........Subtract 1 from NUMPTS_CON since the arrays, XTMP,YTMP are
C          dimensioned from 0 to 4000 and for contour plots we do not
C          have a moveto in the TXTOUT file at the zeroth subscript of
C          the arrays as for rayplots.

            xtmp(0) = 0.0
            ytmp(0) = 0.0
            do ii = 1,numpts_con-7,8
                  READ (iunit,99,ERR=100,IOSTAT=ist)
     &                  xtmp(ii),ytmp(ii),xtmp(ii+1),ytmp(ii+1),
     &                  xtmp(ii+2),ytmp(ii+2),xtmp(ii+3),ytmp(ii+3),
     &                  xtmp(ii+4),ytmp(ii+4),xtmp(ii+5),ytmp(ii+5),
     &                  xtmp(ii+6),ytmp(ii+6),xtmp(ii+7),ytmp(ii+7)
           enddo

           if (MOD(numpts_con,8) .ne. 0) then
                nremain = MOD(numpts_con,8)
                do jj = numpts_con-nremain+1, numpts_con-1, 2
                     READ (iunit,3333,ERR=100,IOSTAT=ist)
     &                     xtmp(jj),ytmp(jj),
     &                     xtmp(jj+1),ytmp(jj+1)
                enddo
           endif
3333                 format (4(I5))

           if (numpts_con .gt. 0) then
                marker = 1    ! no symbol plotting
	        call PLOT_PTS_VT (1655.,9785.,605.,7515.,xtmp,ytmp,
     &                            x, y, numpts_con, marker, rayflag,
     &                            conflag,
     &                            numxticks, numyticks, num_pt_plots,
     &                            slomofact, eigflag)

                edtfactor = 0
                dashline = .false.
                if (rlinewidth .ne. deflwidth) WRITE (2,4444) rlinewidth
                call RAYPLOT_PTS_PS (numpts_con, xtmp, ytmp,
     &                               edtfactor, dashline, conflag)
           endif


C.....ICODE = -2:  Handle rayplot and contour plot frame advances.
C
C     Do a "gsave grestore" so that the user-selected values are maintained.

      elseif (icode .eq. -2) then
          WRITE (2,*) 'gsave'
          WRITE (2,*) 'showpage'
          WRITE (2,*) 'grestore'

          call PAUSEIT
	  call CLR_SCREEN

C.....ICODE = 22 gives the position (in inches) of the lower left corner
C     of the frame from the lower left corner of the page.

      elseif (icode .eq. 22) then

          BACKSPACE iunit
	  READ (iunit,222,IOSTAT=ist,ERR=100) icode, xorig, yorig
222	  FORMAT(I2,2F5.2)

      elseif (icode .eq. 13) then

C.....ICODE = 25 gives the lengths (in inches) of the X and Y axes.
C     Also, the text for the title of the graph and the labels of the
C     horizontal and vertical axes are read.

      elseif (icode .eq. 25) then

         BACKSPACE iunit
	 READ (iunit,333,IOSTAT=ist,ERR=100) icode, x_axis_len,
     &					     y_axis_len
333      FORMAT(I2,2F6.2)

	 READ (iunit,444,IOSTAT=ist,ERR=100) title_len, title
	 READ (iunit,444,IOSTAT=ist,ERR=100) horiz_label_len,
     &					     horiz_label
	 READ (iunit,444,IOSTAT=ist,ERR=100) vert_label_len,
     &					     vert_label
444      FORMAT (I5,1X,A80)

C.....ICODE = 29 provides the number of points to be plotted as well as
C     a MARKER which indicates whether to points are to be connected or are
C     to be indicated individually by a symbol (e.g. a circle).
C     The arrays of X (sound speed) and Y (depth) values are also provided at
C     this point.

      else if (icode .eq. 29) then

	   BACKSPACE iunit
	   READ (iunit,666) icode, numpts, marker
666	   FORMAT (I2,2I5)

	   do i = 1, numpts-1, 2
		READ (iunit,777,IOSTAT=ist,ERR=100) x(i),y(i),
     &						    x(i+1),y(i+1)
777		FORMAT (4E14.6)
	   enddo

C......... Handle the case of an odd number of points.

	   if (MOD(numpts,2) .ne.0) then
		READ (iunit,1111,IOSTAT=ist,ERR=100)
     &					      x(numpts),y(numpts)
1111		FORMAT(2E14.6)
	   endif

C..........The lengths of the X and Y axes (in inches) are given by the data
C          that is associated with ICODE = 25.  Typically, ICODE 22 will precede
C          ICODE 25 in the TXTOUT file.

           if (.not. eigflag .or. (eigflag .and. .not.subheaddun)) then
                call SUBHEAD_PS (raytxt, lentxt, x_axis_len, xorig,
     &                           yorig, y_axis_len, fontsize, marker)

                call SUBHEAD_VT (raytxt, lentxt, marker)

                subheaddun = .true.
           endif

	   call PLOT_PTS_PS (numpts, marker, x, y, xorig, yorig,
     &			     x_axis_len, y_axis_len, xmin, xmax,
     &			     ymin, ymax, rlinewidth, dashline,
     &                       num_pt_plots, eigflag)

           rayflag = .false.
           conflag = .false.

C..........Note that XTMP and YTMP are passed in but are not used
C          for profiles.  Only the X and Y REAL arrays are used.

	   call PLOT_PTS_VT (xmin,xmax,ymin,ymax, xtmp, ytmp,
     &                       x, y, numpts,
     &                       marker, rayflag, conflag,
     &			     numxticks, numyticks, num_pt_plots,
     &                       slomofact, eigflag)

C.....ICODE = 30:  End-of-frame (frame advance) for profile and eigen plots.

      else if (icode .eq. 30) then
           WRITE (2,*) 'showpage'
           WRITE (2,*) 'initclip'

           num_pt_plots = 0

C.........Need to pause here for eigen plots, because we don't pause after
C         every plot as for profiles.

           if (eigflag) call PAUSEIT

	   call CLR_SCREEN

C.....ICODE = 31:  End of all plotting for profile plots.

      else if (icode .eq. 31) then
           WRITE (2,*) 'grestore'
	   WRITE (2,*) '% END'

	   call SETVIDEO_DFT

	   STOP 'NORMAL END'

C.....ICODE = 37 indicates dashed lines to begin at this time.
C     Set the flag which is used by plotting routines.

      elseif (icode .eq. 37) then

           dashline = .true.

C.....ICODE = 38 indicates that the line style will be re-set to solid lines.

      elseif (icode .eq. 38) then

           BACKSPACE iunit
           READ (iunit,999,IOSTAT=ist,ERR=100) icode, resetxt
999        FORMAT (I2,A)

           lentxt = LEN_TRIMMED (resetxt)
           if (resetxt(1:lentxt) .eq. 'DASH') dashline = .false.
           WRITE (2,*) '[] 0 setdash'

C.....ICODE = 39 provides the minimum, maximum and increment values for
C     the X and Y axes.  These values are in real units (e.g. km, km/sec).

      elseif (icode .eq. 39) then

           BACKSPACE iunit
	   READ (iunit,555,IOSTAT=ist,ERR=100) icode, xmin, xinc,
     &		      xmax, ymin, yinc, ymax
555	   FORMAT (I2,6e13.7)

C..........NUMXTICKS and NUMYTICKS are output from DRAW_TICKS_PS.
C	   They are used by LABEL_TICKS_VT via PLOT_PTS_VT.

	   call DRAW_TICKS_PS (xorig,xmin,xinc,xmax,yorig,ymin,yinc,
     &			       ymax, x_axis_len, y_axis_len,
     &			       numxticks, numyticks)

C..........NEW_TITLE etc are output from LABEL_AXES_PS and become input to
C	   LABEL_AXES_VT.

	   call LABEL_AXES_PS (title, horiz_label, vert_label,
     &			    x_axis_len, y_axis_len, xorig, yorig,
     &                      fontsize,
     &			    newtitle, new_horiz_label, new_vert_label)

           call FRAME (x_axis_len, y_axis_len, xorig, yorig)
           WRITE (2,*) 'clip'
           WRITE (2,*) 's'

C..........Note that the frame is "stroked" outside of the FRAME routine.  This
C          enables the "frame" procedure to be called to set up a clipping path.
C          ("Stroke" causes a "newpath" to begin.  Therefore, calling "clip" after
C          a "stroke" essentially clips to a null path, not good.)

C..........Set the line width back to user-selected value since FRAME has set
C          it to the default = .01.

           if (rlinewidth .ne. deflwidth) WRITE (2,4444) rlinewidth
4444       FORMAT (F7.4,' in setlinewidth')

  	   call LABEL_AXES_VT (newtitle, new_horiz_label,new_vert_label,
     &                         numtextrows)

C..........Set back for next eigen plot.
           subheaddun = .false.

C.....ICODE = 40 indicates that an EIGEN plot is being made.  EIGFLAG is
C     sent into PLOT_PTS_VT to prevent there being a pause between every
C     eigen ray.  Also, it provides a way to tell PLOT_PTS_PS to plot
C     a dot instead of a symbol for discrete data point plots.

      elseif (icode .eq. 40) then

           eigflag = .true.
      endif

      go to 10

100   PRINT 888, ist
888   FORMAT (' ERROR: While reading input file. ',I5)
      STOP
C ----------------------------------------------------------------------!-------
      END
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      Subroutine USER_PROMPT (fname, scale, fontsize, rlinewidth,
     &                        edtfactor, outfname, slomofact)
C ----------------------------------------------------------------------!-------
C Prompts the user with a series of questions. 
C ----------------------------------------------------------------------!-------

      CHARACTER*(*) fname, outfname
      REAL scale, edtfactor, rlinewidth, slomofact
      INTEGER fontsize
      CHARACTER*3 advfeat

      call CLR_SCREEN
      WRITE (*,88)
88    FORMAT (8x,'****** PostScript and Screen Graphics Program *****')
      WRITE (*,*)
     &'       for use with PC-HARPO, PROFILE, CONPLT and EIGEN output'
      WRITE (*,99)
99    FORMAT(8X,'****************************************************')
      WRITE (*,*)
      WRITE (*,111)
111   FORMAT (' Input File Name: [TXTOUT] '\)
      READ (*,222) fname
222   FORMAT (A)
      WRITE (*,*)
      WRITE (*,333)
333   FORMAT (' PostScript Output File Name: [PS.OUT] '\)
      READ (*,222) outfname
      WRITE (*,*)
      WRITE (*,444)
444   FORMAT (' Advanced Features ? [N=No] '\)
      READ (*,222) advfeat

      call CLR_SCREEN

      if (advfeat(1:1) .eq. 'y' .or. advfeat(1:1).eq.'Y')
     &     call ADV_FEATURES (fname, scale, fontsize, rlinewidth,
     &                       edtfactor, outfname, slomofact)


      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine USER_PROMPT 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine ADV_FEATURES (fname, scale, fontsize, rlinewidth,
     &                        edtfactor, outfname, slomofact)
C ----------------------------------------------------------------------!-------
C ADVANCED FEATURES 
C ----------------------------------------------------------------------!-------
      REAL scale, edtfactor, rlinewidth, slomofact
      INTEGER fontsize
      CHARACTER*(*) fname, outfname

C.....These two tmp names are used because we don't want to change the
C     file names to defaults (that is done by GET_DEFAULTS).

      CHARACTER*80 tmpfname, tmpoutfname

      if (fname(1:6).eq.'      ') then
           tmpfname = 'txtout'
      else
           tmpfname = fname
      endif
      if (outfname(1:6).eq.'      ') then
           tmpoutfname = 'ps.out'
      else
           tmpoutfname = outfname
      endif

      WRITE (*,*) 'Input File: ',tmpfname(1:LEN_TRIMMED(tmpfname))
      WRITE (*,*) 'PostScript Output File: ',
     &             tmpoutfname(1:LEN_TRIMMED(tmpoutfname))
      WRITE (*,*)
      WRITE (*,*)
      WRITE (*,111)
111   FORMAT(8X,
     &'---------------------------------------------------------------')
      WRITE (*,222)
222   FORMAT(8X,
     &'_____________________  ADVANCED FEATURES ______________________')
      WRITE (*,111)
      WRITE (*,*)
      WRITE (*,333)
333   FORMAT(' Scale Factor for PostScript output: [1.0] '\)
      READ (*,444) scale
444   FORMAT (F6.2)
      WRITE (*,*)
      WRITE (*,555)
555   FORMAT (' Font Size for PostScript axis labels. [12] '\)
      READ (*,666) fontsize
666   FORMAT (I2)
      WRITE (*,*)
      WRITE (*,777)
777   FORMAT (' Line Width for PostScript plots. [.01] inches '\)
      READ (*,888) rlinewidth
888   FORMAT (F7.4)
      WRITE (*,*)
      WRITE (*,999)
999   FORMAT (' Ray-editing Factor for PostScript rayplots:')

10    WRITE (*,1111)
1111  FORMAT (
     &'    -1:  No editing (Highest resolution),'/
     &'    30:  Maximum editing (Lowest resolution). [10] '\)
      READ (*,2222) edtfactor
2222  FORMAT (F3.0)
      if (edtfactor .gt. 30 .or. edtfactor .lt. -1) then
           WRITE(*,*)' INVALID ENTRY. RE-ENTER.'
           goto 10
      endif

      WRITE (*,*)
20    WRITE (*,3333)
3333  FORMAT (' Slow Motion for Ray Plots: ')
      WRITE (*,4444)
4444  FORMAT ('   1:  No Slow Motion'/
     &        '   2 (slow) to 10 (slowest). [1] '\)
      READ (*,2222) slomofact
      if (slomofact .gt. 10. .or.
     &    (slomofact .lt. 1. .and. slomofact .ne. 0.0)) then
           WRITE(*,*)' INVALID ENTRY. RE-ENTER.'
           go to 20
      endif

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine ADV_FEATURES 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine GET_DEFAULTS (fname, scale, fontsize, rlinewidth,
     &                         edtfactor, outfname, slomofact)
C ----------------------------------------------------------------------!-------
C Defaults for input file name, scale zize, font size, line width, 
C editing factor, slow motion factor, and output file name
C ----------------------------------------------------------------------!-------

       CHARACTER*(*) fname, outfname
       REAL scale, rlinewidth, edtfactor, slomofact
       INTEGER fontsize

C      NOTE:  The default for EDTFACTOR is 10.0.  This number is used by
C             RAYPLOT_PTS_PS in determining whether two numbers are more than
C             EDTFACTOR from each other.  Therefore, when this number is 0.0, it
C             means that any two numbers will pass the test of being at least
C             0.0 away from each other, i.e. no editing is done.
C             Thus, another value (other than 0.0) was needed to determine
C             that no editing would done since the default was to be 10.0
C             when nothing (= 0.0) was entered by the user.

       if (fname .eq. '      ') fname = 'TXTOUT'
       if (scale .eq. 0.0) scale = 1.0
       if (fontsize .eq. 0) fontsize = 12
       if (rlinewidth .eq. 0.0) rlinewidth = 0.01
       if (edtfactor .eq. 0.0) edtfactor = 10.0
       if (edtfactor .eq. -1.) edtfactor = 0.0
       if (slomofact .eq. 0.0) slomofact = 1.0
       if (outfname .eq. ' ') outfname = 'PS.OUT'

       RETURN
C ----------------------------------------------------------------------!-------
       END Subroutine GET_DEFAULTS 
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
       Subroutine INIT_PS_FILE (outfname, scale, fontsize, rlinewidth)
C ----------------------------------------------------------------------!-------
C Initialize the postscript file. 
C ----------------------------------------------------------------------!-------

       CHARACTER*(*) outfname
       REAL scale
       INTEGER fontsize

       OPEN (UNIT=2,FILE=outfname,STATUS='UNKNOWN')
       WRITE (2,*) '% PostScript file for raytracing programs.'
       WRITE (2,*) '%'
       WRITE (2,*) 'grestore'   ! in case the printer was in a weird state.
       WRITE (2,*) '/in {72 mul} def'
       WRITE (2,*) '/psunits {72 div} def'
       WRITE (2,*) '/l /lineto load def'
       WRITE (2,*) '/rl /rlineto load def'
       WRITE (2,*) '/m /moveto load def'
       WRITE (2,*) '/rm /rmoveto load def'
       WRITE (2,*) '/s /stroke load def'
       WRITE (2,*) '/np /newpath load def'

C.....Define the symbol drawing procedures for discrete point plots.
C     See PLOT_PTS_PS.

       WRITE (2,*) '/cir { .037 in 0 360 arc stroke} def'
       WRITE (2,*) '/cirs { .025 in 0 360 arc stroke} def'

C......Define the symbol for Eigen point plots.
       WRITE (2,*) '/cirfs { .015 in 0 360 arc fill stroke} def'

C..... NOTE:  The following definition, "cirf" is supplied but is not used.
C             "cirf" and "cirfs" could be used as two extra symbol types for
C             profile point plots.  See subroutine PLOT_PTS_PS.
C             The symbols are a filled circle and a small filled circle.

ccccc  WRITE (2,*) '/cirf { .037 in 0 360 arc fill stroke} def'

       WRITE (2,*) '/sbox { m .025 in 0 in rm '//
     &             ' 0 in .025 in rl '
       WRITE (2,*) ' -.05 in 0 in rl 0 in -.05 in rl'//
     &             ' .05 in 0 in rl'//
     &             ' closepath stroke } def'
       WRITE (2,*) '/lbox { m .035 in 0 in rm '//
     &             ' 0 in .035 in rl '//
     &             ' -.07 in 0 in rl 0 in -.07 in rl'//
     &             ' .07 in 0 in rl'//
     &             ' closepath stroke } def'

       WRITE (2,*) '/strleninches {stringwidth pop psunits} def'
       WRITE (2,*) '/ts {72 171 div mul} def'      ! text scaling: RAYPLOTS.
       WRITE (2,*) '/g {72 1710 div mul} def'      ! data scaling: RAYPLOTS.

C.....Set up the user-defined values.

       WRITE (2,111) rlinewidth
111    FORMAT (F7.4,' in setlinewidth')
       WRITE (2,222) scale, scale
222    FORMAT (2F5.2,' scale')
       WRITE (2,333) fontsize
333    FORMAT ('/Helvetica findfont ',I2,' scalefont setfont')

       WRITE (2,*) '%end INIT_PS_FILE'


       RETURN
C ----------------------------------------------------------------------!-------
       END Subroutine INIT_PS_FILE 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine FRAME (x_axis_len, y_axis_len, xorig, yorig)
C ----------------------------------------------------------------------!-------

C Draw the box that is to contain the plot.  We will later "clip" to this
C box.

      WRITE (2,*) '/frame {'
      WRITE (2,*) '0.01 in setlinewidth'
      WRITE (2,*) 'np'
      WRITE (2,111) xorig, yorig, x_axis_len, y_axis_len,
     &             -x_axis_len
111   FORMAT (F5.2,' in ',F5.2,' in m ',F5.2,' in 0 in rl 0 in ',F5.2,
     & ' in rl ',F5.2,' in 0 in rl closepath } def')

      WRITE (2,*) 'frame'

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine FRAME 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine DRAW_TICKS_PS (xorig,xmin,xinc,xmax,yorig,ymin,yinc,
     &				ymax, x_axis_len, y_axis_len,
     &				numxticks, numyticks)
C ----------------------------------------------------------------------!-------
C
C.....Profile and Eigen Plots Only!
C
C     Draw tick marks and label them.
C
C     NOTE:  Some values for the lengths of lines have been hard-coded
C	     in this routine, e.g. TICKLEN and the distance away from
C	     the axes where the labels for the tick marks are drawn.
C	     These values were chosen by experiment.
C     INPUT:
C            Everything except output
C            NOTE:  xorig, yorig and x_axis_len, y_axis_len are in INCHES.
C                   xmax, xinc, xmin etc are in REAL UNITS.
C     OUTPUT:
C	     numxticks, numyticks
C
C     NOTE:  For non-integer intervals, there are variables with "_f"
C            in their names.

      IMPLICIT none
      INTEGER
     &	     numxticks,
     &	     numyticks,
     &       i,
     &       numdigits,
     &       numdigits_min,
     &       numdigits_inc,
     &       numdigits_max

      REAL
     &       ticklen,
     &	     xorig,xmin,xinc,xmax,
     &	     yorig,ymin,yinc,ymax,
     &	     x_axis_len, y_axis_len,
     &	     in_incr,
     &       lastloc,
     &       fracinterval,
     &       prevyval,
     &       numyticks_f,
     &       prevxval,
     &       numxticks_f,
     &       currval

      LOGICAL
     &       intnumticks,
     &       trunc_flag_min,
     &       trunc_flag_inc,
     &       trunc_flag_max,
     &       trunc_flag

      CHARACTER*5
     &       blanks

      CHARACTER*6
     &       charfmt,
     &       charfmt_min,
     &       charfmt_inc,
     &       charfmt_max

      CHARACTER*10
     &       charnum              !  the tick values are written to
                                  !  this variable using the FORMAT,
                                  !  CHARFMT
                                  !  Then, CHARNUM is written to the PS file.
C==>
      blanks = '     '
      intnumticks = .true.
      numdigits = 0
      charfmt = '      '

      WRITE (2,*) '% begin DRAW_TICKS_PS'
      ticklen = 0.1

      WRITE (2,*) 'np'

C.....First, the X axis.
C
C     Figure out what kinda numbers we're dealing with, i.e.
C     their magnitude and the number of significant digits.
C     Of the max, min and increment, we'll use the number of
C     significant digits which is the largest.

      numdigits_min = 0
      trunc_flag_min = .false.
      numdigits_inc = 0
      trunc_flag_inc = .false.
      numdigits_max = 0
      trunc_flag_max = .false.

C.....If the first argument of SIGNIF_DIGITS is greater than 10**7 or
C     less than 10**-7, then SIGNIF_DIGITS returns with value = 9.
C     Format 444 will be used later.

      if (xmin .ne. 0.0)
     &     call SIGNIF_DIGITS (ABS(xmin), numdigits_min, charfmt_min,
     &                         trunc_flag_min)

      call SIGNIF_DIGITS (ABS(xinc), numdigits_inc, charfmt_inc,
     &                         trunc_flag_inc)

      if (xmax .ne. 0.0)
     &      call SIGNIF_DIGITS (ABS(xmax), numdigits_max, charfmt_max,
     &                         trunc_flag_max)

C.....Remember CHARFMT is of the form:  "(Fn.d)" or "(In)" so the 3rd character
C     is the number of digits.  Count the open parenthesis as character
C     number 1 !

      if (xmin.lt.0.0) then
           numdigits_min = numdigits_min + 1
           write (charfmt_min(3:3),'(I1)') numdigits_min
      endif
      if (xmax.lt.0.0) then
           numdigits_max = numdigits_max + 1
           write (charfmt_max(3:3),'(I1)') numdigits_max
      endif

C.....Determine the maximum number of digits and use that format.

      if (MAX(numdigits_min,numdigits_inc,numdigits_max) .eq.
     &    numdigits_min) then
          charfmt = charfmt_min
          numdigits = numdigits_min
          trunc_flag = trunc_flag_min
      elseif (MAX(numdigits_min,numdigits_inc,numdigits_max) .eq.
     &    numdigits_inc) then
          charfmt = charfmt_inc
          numdigits = numdigits_inc
          trunc_flag = trunc_flag_inc
      elseif (MAX(numdigits_min,numdigits_inc,numdigits_max) .eq.
     &    numdigits_max) then
          charfmt = charfmt_max
          numdigits = numdigits_max
          trunc_flag = trunc_flag_max
      endif

C.....If the min and max are on different sides of zero, add one more
C     digit just to be sure to handle possible problems.

      if (xmin .lt. 0.0 .and. xmax .ge. 0.0) then
          numdigits = numdigits + 1
          write (charfmt(3:3),'(I1)') numdigits
      endif

C.....Truncate the number of ticks (there can be a non-integer number
C     of tick marks e.g. XMIN= -3. XMAX = 0.1 XINC = 1.0
C        --> (XMAX-XMIN)/XINC = 3.1.)  Actual number of ticks desired
C     is 5.  Ticks at -3, -2, -1, 0, 0.1.  If the calculated value
C     had been 3.0, then the actual number of ticks would be 4, there
C     always being n+1 endpoints for n intervals.
C     Set a flag indicating that there are not an integer number of
C     of tick marks for the interval given.  Then, an extra tick will
C     be drawn and labeled.

      numxticks_f = (xmax-xmin)/xinc
      numxticks = AINT(numxticks_f + 1.99)

C.....Determine whether there are an integer number of tick marks
C     and if not, whether rounding by the computer has caused it.

      if (NINT(numxticks_f)-numxticks_f .ne. 0 .and.
     &     ABS(NINT(numxticks_f)-numxticks_f) .ge. .01*xinc) then
           intnumticks = .false.
           numxticks = numxticks - 1
      endif

      in_incr = x_axis_len/numxticks_f      ! inches per tick interval
                                            ! can be non-integer
      do i = 1, numxticks

C..........Move to the tick mark location and draw the line, but don't
C          bother drawing the tick mark for the endpoints.  If there is
C          a non-integer number of tick intervals, we need to draw the
C          NUMXTICKth tick mark.

           if ((i.ne.1 .and. i.ne.numxticks) .or. (i.eq.numxticks.and.
     &         .not.intnumticks))
     &         WRITE (2,111) xorig+(i-1)*in_incr,yorig,
     &	                     xorig+(i-1)*in_incr, yorig+ticklen

C..........Move to the tick mark label location.

	   WRITE (2,222) xorig+(i-1)*in_incr,yorig-0.2

C..........Save the last location (offset from origin) for the possible extra
C          tick mark.  See below.

           lastloc = xorig + (i-1)*in_incr
C
C..........Now label the tick mark.
C          The following WRITE statement simplys move the cursor to a point
C          up to 4 character widths to the left of the tick mark.  The length
C          of the field can be up to 9 characters.  Therefore, we begin drawing
C          the field (i.e. the label) up to 4 characters to the left, thus
C          centering it around the tick mark.

           WRITE (2,2222) blanks(1:numdigits/2+1)

C..........Use F or E format depending on magnitude of tick mark label.

           if (ABS(xmin+(i-1)*xinc) .ge.  1.e+7 .or.
     &         (ABS(xmin+(i-1)*xinc) .lt. .0000001 .and.
     &          ABS(xmin+(i-1)*xinc) .gt. xinc .and.
     &          xmin+(i-1)*xinc .ne. 0.0) ) then

                WRITE (2,444) xmin+(i-1)*xinc
           else
                if (trunc_flag) then
                     write (charnum(1:numdigits),charfmt)
     &                      NINT(xmin+(i-1)*xinc)
                else
                     write (charnum(1:numdigits),charfmt)
     &                      xmin+(i-1)*xinc
                endif
                write (2,555) charnum(1:numdigits)
           endif

           charnum = '          '
      enddo

111	   FORMAT (F5.2,' in ',F5.2,' in m ',F5.2,' in ',F5.2,' in l')
222	   FORMAT (F5.2,' in ',F5.2,' in m')
444        FORMAT('(',e8.2,') show')
555        format ('(',A,') show')
2222       FORMAT (' (',A,') strleninches neg in 0 in rm')


C.....Label the extra tick mark, but don't draw it.

      if (.not. intnumticks) then
           prevxval = xmin + (numxticks-1)*xinc  ! in real units
           fracinterval = (xmax-prevxval)/xinc   ! fraction of the interval

	   WRITE (2,222) lastloc+fracinterval*in_incr, yorig-0.2
           WRITE (2,*) ' (     ) strleninches neg in 0 in rm'
           if (ABS(xmax) .ge.  1.e+5 .or.
     &         (ABS(xmax) .lt. .01 .and. xmax .ne. 0.0) ) then

                WRITE (2,444) xmax
           else
                WRITE (2,333) xmax
           endif

           numxticks = numxticks + 1
      endif

333        FORMAT('(',F8.2,') show')
C ----------------------------------------------------------------------!-------

C Now, the Y axis

      numdigits_min = 0
      trunc_flag_min = .false.
      numdigits_inc = 0
      trunc_flag_inc = .false.
      numdigits_max = 0
      trunc_flag_max = .false.

      if (ymin .ne. 0)
     &     call SIGNIF_DIGITS (ABS(ymin), numdigits_min, charfmt_min,
     &                         trunc_flag_min)

      call SIGNIF_DIGITS (ABS(yinc), numdigits_inc, charfmt_inc,
     &                    trunc_flag_inc)

      if (ymax .ne. 0)
     &      call SIGNIF_DIGITS (ABS(ymax), numdigits_max, charfmt_max,
     &                         trunc_flag_max)

      if (MAX(numdigits_min,numdigits_inc,numdigits_max) .eq.
     &    numdigits_min) then
          charfmt = charfmt_min
          numdigits = numdigits_min
          trunc_flag = trunc_flag_min
      elseif (MAX(numdigits_min,numdigits_inc,numdigits_max) .eq.
     &    numdigits_inc) then
          charfmt = charfmt_inc
          numdigits = numdigits_inc
          trunc_flag = trunc_flag_inc
      elseif (MAX(numdigits_min,numdigits_inc,numdigits_max) .eq.
     &    numdigits_max) then
          charfmt = charfmt_max
          numdigits = numdigits_max
          trunc_flag = trunc_flag_max
      endif

      if (ymin.lt.0 .and. ymax.gt.0) then
           numdigits = numdigits + 1
           write (charfmt(3:3),'(I1)') numdigits
      endif

      intnumticks = .true.
      numyticks_f = (ymax-ymin)/yinc
      numyticks = AINT(numyticks_f + 1.99)

      if (NINT(numyticks_f)-numyticks_f .ne. 0 .and.
     &    ABS(NINT(numyticks_f)-numyticks_f) .ge. .01*yinc) then
           intnumticks = .false.
           numyticks = numyticks-1
      endif

      in_incr = y_axis_len/numyticks_f      ! inches per tick interval

      do i = 1, numyticks

C..........Moveto the tick mark location and draw tick, but don't
C          bother drawing the tick mark for the endpoints.  Draw it
C          if we know that there will be an additional non-integer tick.

           if ((i.ne.1 .and. i.ne.numyticks) .or. (i.eq.numyticks.and.
     &         .not.intnumticks))
     &          WRITE (2,111) xorig, yorig+(i-1)*in_incr,
     &		              xorig+ticklen, yorig+(i-1)*in_incr

           currval = ymin+(i-1)*yinc

C..........If the profile does not extend above the surface of the ocean,
C          change the sign of all the y (depth) values.

           if (ymax .le. 0.0) currval = ABS(currval)

C..........Now label the tick mark.  Note that 0.2 inches was found by
C          experiment.

           WRITE (2,222) xorig-0.3, yorig+(i-1)*in_incr
           WRITE (2,2222) blanks(1:numdigits/2+1)

           if (ABS(currval) .ge. 1.e+7 .or.
     &         ABS(currval) .lt. .0000001  .and.
     &         ABS(currval) .gt. yinc .and.
     &         currval .ne. 0.0) then

                WRITE (2,444) currval

           else if (trunc_flag) then
                write (charnum(1:numdigits),charfmt) NINT(currval)
                write (2,555) charnum(1:numdigits)
           else
                write (charnum(1:numdigits),charfmt) currval
                write (2,555) charnum(1:numdigits)
           endif

           charnum = '          '

C..........Save the last location for the (possible) extra tick mark calculations.

           lastloc = yorig + (i-1) *in_incr

      enddo

C.....Label the extra tick mark, but don't draw it.

      if (.not. intnumticks) then
           prevyval = ymin + (numyticks-1)*yinc
           fracinterval = (ymax-prevyval)/yinc   ! fraction of the interval

           WRITE (2,222) xorig-0.2,
     &                   lastloc + fracinterval * in_incr

           if (ymax .le. 0.0) then
               currval = ABS(ymax)
           else
               currval = ymax
           endif
           if (ABS(ymax) .ge.  1.e+7 .or.
     &         (ABS(ymax) .lt. .00001 .and. ymax .ne. 0.0) ) then
                write(2,444) currval

           else if (trunc_flag) then
                write (charnum(1:numdigits),charfmt) NINT(currval)
                write (2,555) charnum(1:numdigits)
           else
                write (charnum(1:numdigits),charfmt) currval
                write (2,555) charnum(1:numdigits)
           endif

           numyticks = numyticks + 1
      endif

C Now stroke

      WRITE (2,*) 's'
      WRITE (2,*) '% end DRAW_TICKS_PS'

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine DRAW_TICKS_PS 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine SIGNIF_DIGITS (rnum, numdigits, charfmt, trunc)
C ----------------------------------------------------------------------!-------
C
C     This routine determines how many characters are necessary to construct
C     the smallest FORTRAN format edit descriptor (e.g. F4.1) which will
C     write all the significant digits of a given input number, RNUM.
C     It works as follows:
C
C     if RNUM = the truncated value of RNUM  (e.g. RNUM is essentially integer)
C          then   TRUNC is TRUE and an I format is used.
C     else
C           use LOG base 10 to determine how many digits there are.
C           Example, RNUM = 10.5 which would need format F4.1
C
C                  LOG(RNUM) = 1.02119
C                  but we need 2 digits plus some digits after the decimal.
C                  so numdigits = truncated(LOG(RNUM)) + 1
C
C                  Then the fractional part, .5, is examined to determine
C                  how many digits after the decimal, numdecimal, are needed.
C
C                  Increment numdecimal since there is at least 1 decimal
C                  place.
C                  The fractional part is multiplied by 10 and then truncated.
C                  (= 5 in our example).  If the fractional part is now equal
C                  to its truncated value, then we are done.
C
C                  Answer is numdigits = 2, numdecimal = 1.
C                  WAIT!!
C                  Add one digit for the decimal (numdigits = 3 now) AND
C                  add numdecimal to numdigits = 3 + 1 = 4.
C
C                  Then write these values into the appropriate locations
C                  in CHARDIG which is then written into CHARFMT.
C
C           NOTE:  if RNUM is less than 1, it has to be handled differently.
C                  The algorithm above will not work since a floating point
C                  number less than 1 can have numdecimal = numdigits.
C                  See the code below.
C
      logical trunc
      character*4 chardig
      character*6 charfmt

      numdecimal = 0
      numdigits = 0

      trunc = .false.
      if (aint(rnum) .eq. rnum) then
         trunc = .true.
      endif

      if (ABS(rnum) .ge. 1.e+7 .or. ABS(rnum) .le. 1.e-7) then
           numdigits = 9
           RETURN
      endif

      if (.not. trunc) then
           if (rnum.ge.1.0) then
                numdigits = AINT(LOG10(rnum)) + 1
                fracpart = rnum - AINT(rnum)
           endif
           if (rnum.lt.1.0) then
                numdigits = ABS(AINT(LOG10(rnum)))
                numdecimal = numdigits

C.............. Following statement handles RNUM in interval: 0.1 < RNUM < 1.0.

                if (numdigits .eq. 0) numdigits = 1
                if (AINT(LOG10(rnum)) .eq. LOG10(rnum) .or.
     &              ABS(LOG10(rnum)-AINT(LOG10(rnum))) .lt.
     &              1.e-5 ) goto 100
                fracpart = rnum*10**numdigits
           endif

C..........Now, determine how many (up to 7) significant digits there are.
C          When the fractional part = the truncated fractional part, then
C          we're done.  e.g. 5.0 = 5.

           do i = 1,6
                numdigits = numdigits + 1
                numdecimal = numdecimal + 1
                fracpart = fracpart *10
                if (
     &              (AINT(fracpart) .eq. fracpart) .or.
     &              (ABS(fracpart-NINT(fracpart)) .lt. 1.e-5*10**(i-1)))
     &              goto 100
           enddo
      else
           numdigits = AINT(LOG10(rnum)) + 1
           numdecimal = 0
      endif

100   continue

C.....CHARFMT has the parentheses added for the format to be output.

      if (.not. trunc) then
           numdigits = numdigits + 1  ! for the decimal point
           chardig(1:1) = 'f'
           write (chardig(2:2),'(I1)') numdigits
           chardig(3:3) = '.'
           write (chardig(4:4),'(I1)') numdecimal
           charfmt = '('//chardig//')'
      else
           chardig(1:1) = 'i'
           write (chardig(2:2),'(I1)') numdigits
           chardig(3:4) = '  '
           charfmt = '('//chardig//')'
      endif

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine SIGNIF_DIGITS 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------

      Subroutine LABEL_AXES_PS (title, horiz_label, vert_label,
     &			     x_axis_len, y_axis_len, xorig, yorig,
     &                       fontsize,
     &			     newtitle, new_horiz_label,
     &			     new_vert_label)
C ----------------------------------------------------------------------!-------

C This Routine is For PostScript Files for Profile and Eigen Plots Only:
C This routine labels the x and y axes and places the title at the top
C of the graph.
C
C In addition, the labels are stripped of trailing blanks and a routine
C is called which parses the 3 input character strings to determine
C if there are any case changes in the string.	Note that these character
C strings were created (possibly) on a machine (Cyber) which did not have
C lowercase available.	The DISSPLA software then used special characters
C (e.g. # and !) to determine when case changes occurred.  These characters
C are used here to convert the strings to a combination of upper/lower case.
C
C The PostScript procedure strleninches was created to handle centering the
C strings along the axes.
C
C OUTPUT:
C
C	 The title and labels (in their converted form e.g. no special
C	 characters embedded):
C
C	 newTITLE
C	 new_HORIZ_LABEL
C	 new_VERT_LABEL
C
      CHARACTER*(*) title, horiz_label, vert_label
      CHARACTER*80  newtitle, new_horiz_label, new_vert_label

      INTEGER fontsize


      WRITE (2,*) '% begin LABEL_AXES_PS'

C Find out the length of the TITLE, center it and move to left end of it.
C Then output the string.

      call TRIM_STRING (title,newtitle,lennew)


C For TITLE only, set font to a larger size.

      WRITE (2,*)
     &      '/Helvetica findfont ',fontsize,' 1.5 mul scalefont setfont'

      WRITE (2,111) newtitle(1:lennew)
111   FORMAT ('(',A,') strleninches')

C m x_axis_len/2 - (title length/2), taking units and origin into account

      WRITE (2,222) x_axis_len/2 + xorig, yorig + y_axis_len + 1.0
222   FORMAT (' neg 2 div ',F6.2,' add in ',F6.2,' in m')

      WRITE (2,333) newtitle(1:lennew)
333   FORMAT ('(',A,') show')

C Set font back to original size.

      WRITE (2,*) '/Helvetica findfont ',fontsize,' scalefont setfont'

C As above for the HORIZ_LABEL.

      call TRIM_STRING (horiz_label,new_horiz_label,lennew)

      WRITE (2,111) new_horiz_label(1:lennew)
      WRITE (2,222) x_axis_len/2 + xorig, yorig - 0.8
      WRITE (2,333) new_horiz_label(1:lennew)

C Now the VERT_LABEL.

      WRITE (2,444) xorig - 0.8
444   FORMAT (F6.2,' in ')

      call TRIM_STRING (vert_label,new_vert_label,lennew)

      WRITE (2,111) new_vert_label(1:lennew)
      WRITE (2,555) y_axis_len/2 + yorig
555   FORMAT (' neg 2 div ',F6.2,' add in m')

      WRITE (2,*) 'gsave'
      WRITE (2,*) '90 rotate'
      WRITE (2,333) new_vert_label(1:lennew)
      WRITE (2,*) 'grestore'

      WRITE (2,*) '% end LABEL_AXES_PS'

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine LABEL_AXES_PS 
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------

      Subroutine TRIM_STRING (oldstring, newstring, lennew)
C ----------------------------------------------------------------------!-------

C Trim the trailing blanks from the OLDSTRING and put them in the NEWSTRING.
C Output the LENNEW so that the NEWSTRING will be written to the PS
C file sans trailing blanks.
C
C The $ character at the end of the string is another DISSPLA feature which
C indicates the end of the string.
C

      CHARACTER*(*) oldstring
      CHARACTER*80  newstring

      INTEGER lenold

      EXTERNAL len_trimmed

      lenold = LEN_TRIMMED(oldstring)

C.....Avoid subscript out of range (i.e. = 0) when oldstring is all blanks.

      if (lenold .gt. 0) then
           if (oldstring(lenold:lenold) .eq. '$') lenold = lenold-1
      endif

      call CNVRT_LOWER (oldstring,lenold,newstring,lennew)

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine TRIM_STRING 
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------

      INTEGER FUNCTION LEN_TRIMMED (oldstring)
C ----------------------------------------------------------------------!-------

C This function determines the length of OLDSTRING without any trailing
C blanks.
C
C Some FORTRAN compilers have built-in functions which accomplish the
C same thing.
C

      CHARACTER*(*) oldstring

      len_trimmed = LEN(oldstring)
      do i = LEN(oldstring),1,-1
	   if (oldstring(i:i) .eq. ' ') then
	       len_trimmed = len_trimmed-1
	   else
	       RETURN
	   endif
      enddo

      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION LEN_TRIMMED 
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------

      Subroutine CNVRT_LOWER (oldstr, lenold, newstr, lennew)
C ----------------------------------------------------------------------!-------

C
      INTEGER lennew,
     &	      pound_idx,
     &	      exclam_idx,
     &	      lenold,
     &	      idx_tmp


      CHARACTER*(*) oldstr
      CHARACTER*80 newstr
      LOGICAL found_low

C Handle case of not having any lower case letters in the string.  When
C that occurs, need to set LENNEW to LENOLD.

      found_low = .false.

C Loop until no more lower case symbols are found.

10    pound_idx = INDEX (oldstr,'#')
      exclam_idx = INDEX (oldstr,'!')


      if (pound_idx .gt. 0) then

	   lennew = 0
	   found_low = .true.

	   do i = 1,pound_idx-1
		newstr(i:i) = oldstr(i:i)
		lennew = lennew + 1
	   enddo

	   if (exclam_idx.gt.0) then
	     do i = pound_idx+1, exclam_idx-1
		lennew = lennew+1
		if (oldstr(i:i) .ge.'A' .and. oldstr(i:i) .le.'Z') then
		     idx_tmp = ICHAR(oldstr(i:i)) + 32
		     newstr(lennew:lennew) = CHAR(idx_tmp)
		else
		     newstr(lennew:lennew) = oldstr(i:i)
		endif
	     enddo

	     do i = exclam_idx+1, lenold
		lennew = lennew + 1
		newstr(lennew:lennew) = oldstr(i:i)
	     enddo
	  endif
      else
C..........Handle the case of no lower case symbols in the string.

	   if (.not. found_low) then
		lennew = lenold
		do i = 1, lennew
		     newstr(i:i) = oldstr(i:i)
		enddo
C...............Blank fill the rest of the new string.  If this is not done,
C               a later string (which is shorter) will have some of the charac-
C               ters from this current string at the end of it.

c 22 Aug 2008		do i = lennew+1, 80
c 22 Aug 2008		     newstr(i:i) = ' '
c 22 Aug 2008		enddo
c The string can be blank filled without a Do loop (complements of Larry Baker)
		If(lennew.lt.LEN(newstr)) Then                     ! 22 Aug 2008
		newstr(lennew+1:) = ''                             ! 22 Aug 2008
		End If                                             ! 22 Aug 2008
	   endif

C..........Normal return point when all lower case have been found.

	   RETURN
      endif

      do i = 1, lennew
	   oldstr(i:i) = newstr(i:i)
      enddo

C.....Blank fill the rest of the new string.  See comment above.

c 22 Aug 2008      do i = lennew+1, 80
c 22 Aug 2008	   newstr(i:i) = ' '
c 22 Aug 2008      enddo
c The string can be blank filled without a Do loop (complements of Larry Baker)
		If(lennew.lt.LEN(newstr)) Then                     ! 22 Aug 2008
		newstr(lennew+1:) = ''                             ! 22 Aug 2008
		End If                                             ! 22 Aug 2008

      lenold = lennew

      go to 10
      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine CNVRT_LOWER 
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      Subroutine PLOT_PTS_PS (numpts, marker, x, y, xorig, yorig,
     &			      x_axis_len, y_axis_len,
     &			      xmin, xmax, ymin, ymax,
     &                        rlinewidth, dashline, num_pt_plots,
     &                        eigflag)
C ----------------------------------------------------------------------!-------

C.....For Profile Plots and Eigen Range-Elev Angle or Range-Time Travel Plots
C     Only!
C.....This routine plots curves (actually points connected by a series
C     of "lineto"s) and discrete data points (MARKER .lt 0) to the PostScript
C     file.  When discrete data points are plotted, they are centered at
C     a symbol which will be one of the following:
C          1. Open circle    (small radius)
C          2. Open box       (small size)
C          3. Open circle    (large radius)
C          4. Open box       (large size)
C          5. Filled circle  (small radius)   NOT IMPLEMENTED (see INIT_PS_FILE)
C          6. Filled circle  (large radius)   NOT IMPLEMENTED (see INIT_PS_FILE)
C
C     The symbols will be used for each consecutive point-plot in the order
C     given above.

      INTEGER numpts,
     &        marker,
     &        num_pt_plots,
     &        fmtlabel

      REAL    x(numpts), y(numpts),
     &	      scalgx, scalgy,
     &        rlinewidth,
     &        deflwidth

      PARAMETER (deflwidth = .01)

      LOGICAL dashline,
     &        eigflag

C.....Handle case of all values being equal.

      if (xmax .eq. xmin) then
	   scalgx = 1
	   scalgy = 1
      else
	   scalgx = x_axis_len/(xmax-xmin)
	   scalgy = y_axis_len/(ymax-ymin)
      endif

      WRITE (2,*) 'np'

C.....Handle plotting discrete points with a symbol for Profile plots.

      if (marker .lt. 0 .and. .not. eigflag) then

C..........Temporarily set line width to a very small width

           WRITE (2,*) '.001 in setlinewidth'

           num_pt_plots = num_pt_plots + 1

C..........Determine which symbol (i.e. how many point plots have
C          already been made) to use.  The symbol type determines
C          which FORMAT statement will be used (and vice versa).

           if (num_pt_plots .eq. 1) ASSIGN 111 TO fmtlabel
           if (num_pt_plots .eq. 2) ASSIGN 222 TO fmtlabel
           if (num_pt_plots .eq. 3) ASSIGN 333 TO fmtlabel
           if (num_pt_plots .eq. 4) ASSIGN 444 TO fmtlabel

	   do i = 1, numpts-1, 2
		WRITE (2,fmtlabel) xorig+((x(i)-xmin)*scalgx),
     &			yorig+((y(i)-ymin)*scalgy),
     &			xorig+((x(i+1)-xmin)*scalgx),
     &			yorig+((y(i+1)-ymin)*scalgy)
111             FORMAT (G10.3,' in ',G10.3,' in cirs ',            ! 15 Aug 2011
     &                  G10.3,' in ',G10.3,' in cirs')
222             FORMAT ('np ',G10.3,' in ',G10.3,' in sbox ',
     &                  'np ',G10.3,' in ',G10.3,' in sbox ')
333             FORMAT (G10.3,' in ',G10.3,' in cir ',             ! 15 Aug 2011
     &                  G10.3,' in ',G10.3,' in cir')
444             FORMAT ('np ',G10.3,' in ',G10.3,' in lbox ',
     &                  'np ',G10.3,' in ',G10.3,' in lbox ')
C Replaced F6.3 with G10.3 in above lines     2 Sept 2004
	   enddo

           if (mod(numpts,2).ne.0)
     &         WRITE (2,fmtlabel) xorig+((x(numpts)-xmin)*scalgx),
     &			yorig+((y(numpts)-ymin)*scalgy)

C..........Reset the line width to user-selected value.

           WRITE (2,555) rlinewidth
555        FORMAT (F7.4,' in setlinewidth')

C.....Eigen range-elevation and range-time plots with discrete points

      else if (marker .lt. 0 .and. eigflag) then

	   do i = 1, numpts-1, 2
		WRITE (2,666) xorig+((x(i)-xmin)*scalgx),
     &			      yorig+((y(i)-ymin)*scalgy),
     &		              xorig+((x(i+1)-xmin)*scalgx),
     &			      yorig+((y(i+1)-ymin)*scalgy)
666             FORMAT (G10.3,' in ',G10.3,' in cirfs ',           ! 15 Aug 2011
     &                  G10.3,' in ',G10.3,' in cirfs')
C Replaced F6.3 with G10.3 in above lines      2 Sept 2004
           enddo

           if (mod(numpts,2).ne.0 .or. numpts .eq. 1)
     &         WRITE (2,666) xorig+((x(numpts)-xmin)*scalgx),
     &		             yorig+((y(numpts)-ymin)*scalgy)


C.....Handle plotting continuous curves.

      else

           if (rlinewidth .ne. deflwidth) WRITE (2,555) rlinewidth
           if (dashline) WRITE (2,*) '[4 4] 0 setdash'

	   WRITE (2,777) xorig+((x(1)-xmin)*scalgx),
     &			yorig+((y(1)-ymin)*scalgy)
777	   FORMAT (G10.3,' in ',F6.3,' in m') ! 27 June 2004

	   do i = 2, numpts-2, 3
	       WRITE (2,888) xorig+((x(i)-xmin)*scalgx),
     &			yorig+((y(i)-ymin)*scalgy),
     &			xorig+((x(i+1)-xmin)*scalgx),
     &			yorig+((y(i+1)-ymin)*scalgy),
     &			xorig+((x(i+2)-xmin)*scalgx),
     &			yorig+((y(i+2)-ymin)*scalgy)
888	  FORMAT (G10.3,' in ',F6.3,' in l ',G10.3,' in ',F6.3,' in l ',
     &		  G10.3,' in ',F6.3,' in l') ! 27 June 2004

	   enddo

C..........Handle the extra 1 or 2 points that were not written by the loop
C          above.

           if (MOD(numpts-1,3).ne.0) then
                nremain = MOD(numpts-1,3)
                do j = nremain,1,-1
                    WRITE (2,999) xorig+((x(numpts-j+1)-xmin)*scalgx),
     &                            yorig+((y(numpts-j+1)-ymin)*scalgy)
                enddo
           endif
999        FORMAT (G10.3,' in ',f6.3,' in l') ! 27 June 2004
	   WRITE (2,*) 's'
      endif

C.....Reset to solid lines so as not to affect the frame, tick marks etc.

      if (dashline) WRITE (2,*) '[] 0 setdash'

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine PLOT_PTS_PS 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine RAYPLOT_PTS_PS (numpts, x, y, edtfactor, dashline,
     &                           conflag)
C ----------------------------------------------------------------------!-------


C This routine handles plotting the line graphs for rayplots only.
C Each x,y pair is scaled for PostScript plots by the "g" procedure
C which is defined by INIT_PS_FILE.  "g" causes the input Harpo plotting
C units to be scaled such that the graph and the accompanying text fill
C about 2/3 of an 8.5x11" sheet of paper (when a PostScript "scale" factor
C of 1.0 is used).
C
      INTEGER numpts
      INTEGER*2
     &        x(0:numpts), y(0:numpts)
      REAL
     &        edtfactor

      INTEGER fmtlabel, fmtlabel2

      LOGICAL
     &        dashline,
     &        conflag

C.....Note that the values in X and Y run from 0 to 10240.

C.....OUTX and OUTY are local arrays which are used to reduce the number of
C     points which are written to the PostScript file.  When 2 points are
C     within EDTFACTOR of each other, one of them is not written.  This
C     can be done without loss of resolution, depending on the value of
C     EDTFACTOR.
C
      INTEGER*2 OUTX(0:4000),OUTY(0:4000)



      WRITE (2,*) 'np'

      if (dashline) WRITE (2,*) '[4 4] 0 setdash'



C.....Recall that contour plots have alternating "moveto,lineto" pairs,
C     thus use a different FORMAT statement for contour plots.

      if (conflag) then
           ASSIGN 333 to fmtlabel
           ASSIGN 555 to fmtlabel2
           numremain = MOD(numpts,4)
      else
           WRITE (2,111) x(0), y(0)
111        FORMAT (I5,' g ',I5,' g m')
           ASSIGN 444 to fmtlabel2
           ASSIGN 222 to fmtlabel
           numremain = MOD(numpts,4)
      endif


C.....Write the linetos to the PS file 4 pairs (8 values) of X,Y per row.


C.....Handle contour plots and small ray plots with no editing the same way.
C     In general, each contour line (which is the same as a plot) comprises
C     only a few hundred points or less and actually there are only half
C     the number of NUMPTS because NUMPTS includes one point for a moveto and
C     one point for a lineto.  Also, the editing causes very little change
C     to the number of points due to the arrangement of the points in contour
C     plots.

      if (conflag .or. (numpts .ge. 4 .and. numpts .le. 100)) then
           do i = 1, numpts-3, 4
                WRITE (2,fmtlabel) x(i), y(i),
     &                        x(i+1), y(i+1),
     &                        x(i+2), y(i+2),
     &                        x(i+3), y(i+3)
222             FORMAT (I5,' g ',I5,' g l ',I5,' g ',I5,' g l ',
     &                  I5,' g ',I5,' g l ',I5,' g ',I5,' g l')

333             FORMAT (I5,' g ',I5,' g m ',I5,' g ',I5,' g l ',
     &                  I5,' g ',I5,' g m ',I5,' g ',I5,' g l')

                if (MOD(i,500) .eq. 0) WRITE (2,*) 's'
           enddo

           if (numremain .gt. 0) then
                if (.not. conflag) then
                     do ii = 1, numremain
                          WRITE (2,fmtlabel2)
     &                        x(numpts-numremain+ii),
     &                        y(numpts-numremain+ii)
444                  FORMAT (I5,' g ',I5,' g l')
555                  FORMAT (I5,' g ',I5,' g m ',I5,' g ',I5,' g l')

                     enddo
                else
                     do ii = 1, numremain-1, 2
                          WRITE (2,fmtlabel2)
     &                        x(numpts-numremain+ii),
     &                        y(numpts-numremain+ii),
     &                        x(numpts-numremain+ii+1),
     &                        y(numpts-numremain+ii+1)
                     enddo
                endif
           endif
           WRITE (2,*) 's'

      else if (numpts .lt. 4) then   ! Less than 4 points
           do ii = 1, numpts
                WRITE (2,444) x(ii), y(ii)
           enddo
           WRITE (2,*) 's'

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C.....Handle large plots.  Apply the editing factor.
C
C     NOTE:  The 4 lines starting at ".or." below can be commented out.
C            The PostScript file size will be reduced but the handling
C            of the peaks of the upturns and downturns for rayplots will
C            suffer.  The problem is:  the improvement is minimal but the
C            file size is increased significantly maybe 5-10%.  We chose
C            to have improved plots rather than smaller files.

      else if (numpts .gt. 100 ) then
           outx(0) = x(0)
           outx(1) = x(1)
           outy(0) = y(0)
           outy(1) = y(1)
           kk = 1
           do i = 2, numpts-1
                if (  (ABS(x(i)-outx(kk)) .ge. edtfactor .and.
     &                 ABS(y(i)-outy(kk)) .ge. edtfactor)
     &             .or.
     &             (y(i-1) .lt. y(i) .and. y(i) .ge. y(i+1))  ! downturn
     &             .or.
     &             (y(i-1) .gt. y(i) .and. y(i) .le. y(i+1))  ! upturn
     &           ) then

                     kk = kk + 1
                     outx(kk) = x(i)
                     outy(kk) = y(i)
                endif
           enddo
           kk = kk+1
           outx(kk) = x(numpts)
           outy(kk) = y(numpts)
           numremain = MOD (KK,4)
           do i = 1, kk-3, 4
                WRITE (2,fmtlabel) outx(i), outy(i),
     &                        outx(i+1), outy(i+1),
     &                        outx(i+2), outy(i+2),
     &                        outx(i+3), outy(i+3)

C...............Handle the PostScript limit for the maximum number of
C               points in a path (= 1500).  The "i+3"rd value is always
C               a "moveto", therefore by invoking a "newpath" we are
C               breaking up a path unless we write the "moveto" again.
C               Use 1200 to give a lot of slack (since the 1500 limit
C               includes clipping path as well as the current path).

                if (MOD(i+3,1200) .eq. 0) then
                     WRITE (2,*) 's'
                     WRITE (2,666) outx(i+3), outy(i+3)
666                  FORMAT (I5,' g ',I5,' g m ')
                endif
           enddo

           if (numremain .gt. 0) then
                     do ii = 1, numremain
                          WRITE (2,444)
     &                        outx(kk-numremain+ii),
     &                        outy(kk-numremain+ii)
                     enddo

           endif

           WRITE (2,*) 's'

      endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      if (dashline) WRITE (2,*) '[] 0 setdash'

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine RAYPLOT_PTS_PS 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
C      Subroutine RAYTXT_OUT_PS (raytxt, xtxt, ytxt, iorientn, lentxt)
      Subroutine RAYTXT_OUT_PS (raytxt, xtxt, ytxt, iorientn, lentxt
     1 ,fontsize,int,ital,isize)                                  ! 22 Jun 2009
C ----------------------------------------------------------------------!-------

C.....Handles text output for rayplots and contour plots.
C     IORIENTN = 1 -->  vertical orientation
C     IORIENTN = 0 -->  horizontal orientation

      REAL xtxt, ytxt
      CHARACTER*(*) raytxt
      CHARACTER*20 blanks
      INTEGER iorientn, lentxt
      INTEGER fontsize                                            ! 22 Jun 2009

      data blanks /'                    '/

      if(isize .eq. 0) WRITE (2,*)
     &'/Helvetica findfont ',fontsize,' .75 mul scalefont setfont' ! 23 Jun 2009

      if(isize .eq. 2) WRITE (2,*)
     &'/Helvetica findfont ',fontsize,' 1.5 mul scalefont setfont' ! 23 Jun 2009

      if(isize .eq. 3) WRITE (2,*)
     &'/Helvetica findfont ',fontsize,' 2.0 mul scalefont setfont' ! 22 Jun 2009

C.....First, Vertical Orientation

      if (iorientn .eq. 1) then
           WRITE (2,222) NINT(xtxt), NINT(ytxt)
222        FORMAT ('gsave ',I5,' ts ',I5,' ts translate')
           WRITE (2,*) '0 0 moveto 90 rotate'

           WRITE (2,333) raytxt(1:lentxt)
333        FORMAT ('(',A,') show')
           WRITE (2,*) 'grestore'

C.....Horizontal Orientation
C.....Put the title and the MODELS up 2.0 inches higher on the page.

      elseif (NINT(xtxt) .eq. 0 .or. NINT(xtxt) .eq. 1050) then
           WRITE (2,*) 'gsave 0.0 in 2.0 in translate'
           WRITE (2,111) NINT(xtxt), NINT(ytxt)
           WRITE (2,333) raytxt(1:lentxt)
           WRITE (2,*) 'grestore'

      else
           WRITE (2,111) NINT(xtxt), NINT(ytxt)
111        FORMAT (I5,' ts ',I5,' ts m')

C..........Trim off leading blanks so the text is better located.

           do i=1,lentxt
                blank = INDEX(raytxt(i:lentxt),' ')

                if (blank .eq. 1) then   ! leading blanks will always result
                                     ! in BLANK = 1 since raytxt is i:lentxt
                                     ! Embedded blanks will ==> BLANK > 1.
                    lastblank = i
                else if (blank .gt. 1) then ! No leading blanks, only embedded blanks
                    lastblank = 0
                    goto 10
                endif
           enddo

10         continue

C..........For small text, e.g. tick labels, move the text to the right depend-
C          on how long the string is (this is done because CONPLT and HARPO use
C          the same location to start writing the tick labels regardless of the
C          length of the string, always uses 6 places).  Therefore, we do the
C          moving here.

C..........Use all the non-blank text i.e. skip the leading blanks and trailing
C          blanks.

C..........If the label is e.g. ".5", change it to "0.5".

           if (raytxt(lastblank+1:lastblank+1) .eq. '.') then
                raytxt(lastblank:lastblank) = '0'
                if (lentxt .le. 6) then
C.....................Have to add 1 since we just added a '0'
                      lenmove = 6 - (lentxt+1-lastblank)
                      WRITE (2,444) blanks(1:lenmove)
444                   FORMAT (' (',A,') strleninches in 0 in rm')
                endif
                WRITE (2,333) raytxt(lastblank:lentxt)


C.....Scale Subheading to the same size as the fontsize chosen by the
C     user, smaller than TITLE which is 1.5 * fontsize.           ! 22 Jun 2009

      WRITE (2,*)
     &     '/Helvetica findfont ',fontsize,' scalefont setfont'

                RETURN
           endif

C..........Now handle those that do not contain a decimal point as first char.
           if (lentxt .le. 6) then
                lenmove = 6 - (lentxt-lastblank+1)
                WRITE (2,444) blanks(1:lenmove)
                WRITE (2,333) raytxt(lastblank+1:lentxt)

C..........Finally, handle those without decimal points as first or last and
C          are also more than 6 characters.  E.g., RANGE AT SEA LEVEL or
C          0.0 Deg. E

           else
                WRITE (2,333) raytxt(1:lentxt)
           endif
      endif

C.....Scale back to the same size as the fontsize chosen by the user, 
C     larger than for ISIZE=0, smaller than for isize=2 or 3.     ! 23 Jun 2009

      IF(ISIZE.EQ.0.OR.ISIZE.EQ.2.OR.ISIZE.EQ.3) WRITE (2,*)
     &     '/Helvetica findfont ',fontsize,' scalefont setfont'

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine RAYTXT_OUT_PS 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine SUBHEAD_PS (subtext, lentxt, x_axis_len, xorig, yorig,
     &                       y_axis_len, fontsize, marker)
C ----------------------------------------------------------------------!-------

C.....Handles subheading for regular profiles as well as subheading
C     for point plots (i.e. marker not equal to 0).  The point plot
C     subheadings are placed lower on the page than the other subheading,
C     thus, both subheadings can appear on the same page.

      CHARACTER*(*) subtext

C.....LENTXT is the length of SUBTEXT minus the trailing blanks
      INTEGER lentxt, fontsize

C     Place the Subheading just below the TITLE.

      WRITE (2,*) '% Subheading'

C.....Have to reset the clipping region to device's initial value since
C     at this point, the clipping region has been set to include
C     only the plotting area for the points.

      WRITE (2,*) 'gsave'
      WRITE (2,*) 'initclip'

C.....Handle Subheading.  Note:  Trailing blanks have already been removed
C     from SUBTEXT prior to entry to this routine.

C.....Scale Subheading to the same size as the fontsize chosen by the
C     user, smaller than TITLE which is 1.5 * fontsize.

      WRITE (2,*)
     &     '/Helvetica findfont ',fontsize,' scalefont setfont'

C.....Get length of SUBTEXT in PS units.

      WRITE (2,111) subtext(1:lentxt)
111   FORMAT ('(',A,') strleninches')

C.....Center SUBTEXT.
C     moveto x_axis_len/2 - (title length/2), taking units and origin into
C     account.

      if (marker .ge. 0) then
           WRITE (2,222) x_axis_len/2 + xorig, yorig + y_axis_len + 0.5
222        FORMAT (' neg 2 div ',F6.2,' add in ',F6.2,' in m')
      else
           WRITE (2,222) x_axis_len/2 + xorig, yorig + y_axis_len + 0.3
      endif
      WRITE (2,333) subtext(1:lentxt)
333   FORMAT ('(',A,') show')

C.....Restore previously set clipping region.

      WRITE (2,*) 'grestore'
      WRITE (2,*) '% End Subheading'

      RETURN
C ----------------------------------------------------------------------!-------
      END Subroutine SUBHEAD_PS 
C ----------------------------------------------------------------------!-------
