C 
C Program GRDPS to Display GRiD values in Post Script
C
C Version: 5.40
C Date: 2000, April 19
C
C Coded by Ludek Klimes
C             Department of Geophysics, Charles University Prague
C             Ke Karlovu 3,  121 16  Praha 2,  Czech Republic
C             E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of the input and output files:
C     GRD='string' ... Name of the input ASCII file with the grid values
C             to be plotted.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Default: GRD='grd.out'
C     PS='string' ... Name of the output PostScript file.
C             If non-blank PSHEX is given, output file PS contains just
C             the prolog part of the output PostScript file, see PSHEX.
C             If the value of parameter N4 is greater than 1, PS is
C             the output filename for the first spatial grid.
C             The filenames for the subsequent spatial grids are
C             generated in such a way that all digits contained within
C             the filename are considered a number which is increased
C             by 1 for each new time slice.
C             Default: PS='grd.ps'
C     PSHEX='string' ... Optional name of the output file containing
C             hexadecimal grid data and the PostScript trailer.
C             If PSHEX is blank (default), file PS contains the whole
C             PostScript file.
C             Otherwise, file PS contains just the prolog, setup, and
C             object definition parts, whereas the grid data are written
C             to file PSHEX.
C             The separation into two files is designed to facilitate
C             batch or manual editing of the PostScript definitions.
C             Resulting PostScript file is then obtained by merging
C             PS, possible new definitions and PSHEX.
C             If the value of parameter N4 is greater than 1, PSHEX is
C             the output filename for the first spatial grid.
C             The filenames for the subsequent spatial grids are
C             generated in such a way that all digits contained within
C             the filename are considered a number which is increased
C             by 1 for each new time slice.
C             Default: PSHEX=' '
C     GRDPS2='string', GRDPS3='string', ... , GRDPS9='string' ... Names
C             of the optional input ASCII files with the grid values to
C             be plotted together with data of GRD.
C             For example, data GRD may control hue (colour) and data
C             GRDPS2 may control brightness or saturation.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Defaults: GRDSP2=' ', ... , GRDPS9=' '
C Data specifying grid dimensions:
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C     N4=positive integer... Number of time slices (snapshots).
C             Individual time slices from the input file 'FILEIN' are
C             separated into different files.  The name(s) of the first
C             file(s) is(are) given by input parameters 'PS' and
C             'PSHEX' described above.  Next filenames are generated in
C             such a way that all digits contained within the filenames
C             are considered as numbers which are increased by 1 for
C             each new time slice.
C             Default: N4=1
C     D4=real... Time interval between two consecutive time slices.
C             Useful only if explicitly specifying time-dependent colour
C             scaling of snapshots (see input parameters VCIRC4, VSAT4
C             and VBRI4 below).
C             Default: D4=1.
C Additional data specific to this program:
C (a) Dimensions and layout:
C     NH=positive integer... Layout of 2-D sections of a 3-D grid.
C             Individual 2-D sections of dimensions N1*N2 are merged
C             into rows NH sections long.  Resulting (N3-1)/NH+1 rows
C             are merged together to form the final image.
C             Default: NH=INT(SQRT(N3))
C     UNIT='string'... All lengths controlling the size and position of
C             the plot are assumed to be expressed in the units given
C             by the string.  The units also influence the default
C             paper size, plot size and margins.  Allowed values:
C             UNIT='cm': centimetres (default),
C             UNIT='in': inches (1in=2.54cm),
C             UNIT='pt': big points (1in=72pt).
C             Points (pt) are useful to generate plots for conversion to
C             bitmap forms, e.g., using GhostScript.
C             Default: UNIT='cm'
C     XSIGN=real... Determines the sign of the default value of HSIZE.
C             Default: XSIGN=1.
C     HSIZE=real... Size (in UNITs) of the image, corresponding to the X1
C             axis (horizontal before a possible rotation).
C             If negative, the values will be displayed from the right
C             to the left.
C             Default: HSIZE=SIGN( 16.0,XSIGN) for UNIT='cm',
C                      HSIZE=SIGN(  6.5,XSIGN) for UNIT='in',
C                      HSIZE=SIGN(NH*N1,XSIGN) for UNIT='pt'.
C     YSIGN=real... Determines the sign of the default value of VSIZE.
C             Default: YSIGN=1.
C     VSIZE=real... Size (in UNITs) of the image, corresponding to the X2
C             axis (vertical before a possible rotation).
C             If negative, the values will be displayed from the top to
C             the bottom.
C             Default (proportional display):
C               VSIZE=SIGN(HSIZE*((N3-1)/NH+1)*N2/(NH*N1),YSIGN), i.e.,
C               VSIZE=SIGN(HSIZE*N2/N1,YSIGN) for N3=1.
C     HOFFSET=real... Distance (in UNITs) of the image from the leftmost
C             paper edge (before a possible rotation).  Controls the
C             horizontal position of the figure.
C             Default: HOFFSET=2.5 for UNIT='cm',
C                      HOFFSET=1.0 for UNIT='in',
C                      HOFFSET=0.0 for UNIT='pt'.
C     VOFFSET=real... Distance (in UNITs) of the image from the bottom
C             paper edge (before a possible rotation).  Controls the
C             vertical position of the figure.
C             Default:
C               if VSIZE.LE.HEIGHT-2*2.5:     VOFFSET=HEIGHT-2.5-VSIZE
C               otherwise if VSIZE.LE.HEIGHT: VOFFSET=(HEIGHT-VSIZE)/2.
C               otherwise:                    VOFFSET=2.5
C     HEIGHT=real... Height of the paper in a portrait position.
C             Default: HEIGHT=29.7 for UNIT='cm',
C                      HEIGHT=11.0 for UNIT='in',
C                      HEIGHT=((N3-1)/NH+1)*N2 for UNIT='pt'.
C     ROTATE=real... Enables to rotate the image by angle specified in
C             degrees (positive counterclockwise).  The image is rotated
C             around the centre of the square paper of size HEIGHT.
C             If applied, the user will probably wish to specify the
C             value of ROTATE=90.
C             Parameters HSIZE,VSIZE,HOFFSET,VOFFSET apply to the image
C             before rotation.
C             Attention: BoundingBox is incorrect if ROTATE is not
C               multiple of 90 degrees.
C             Default: ROTATE=0.
C (b) Range of displayed values:
C     VMIN=real, VMAX=real... Values less than or equal to VMIN,
C             or greater than or equal to VMAX will be deemed
C             undefined.  For example, if VMIN=0 is set when plotting
C             velocities, zero velocities are rendered as undefined
C             values, that is usually desirable.
C             Defaults: VMIN=-999999, VMAX=999999000000.
C     R=real, G=real, B=real... Colour of the undefined
C             values.
C             Defaults: R=0.80, G=0.80, B=0.80 (light grey)
C (c) Colour interpretation of values contained within file GRD:
C     VPLUS=real, VSIGN=real, VCIRC=real... VCIRC is the extent of
C             values corresponding to the whole colour circle RGB.
C             Negative VCIRC corresponds to the opposite direction
C             around the colour circle (BGR).
C             Undefined values are displayed in gray.
C             Defaults: VPLUS=0, VSIGN=1,
C               VCIRC=SIGN(GMAX1-GMIN1+VPLUS,VSIGN),
C               where GMINi and GMAXi are the minimum and maximum grid
C               values defined in file 'GRDi'.
C               If the above value of VCIRC is zero (e.g., if
C               GMAX1=GMIN1), the default value of VCIRC is changed to
C               VCIRC=1.
C             Hint: Specify VPLUS=1 when plotting integer values.
C             Hint: To plot velocities with blue minimum velocity,
C               increasing through green, yellow and red again to blue,
C               specify VMIN=0 and VSIGN=-1 with default VREF and CREF.
C     VREF=real, CREF=real... Value VREF corresponds to colour CREF,
C             where Red=0, Yellow=1/6, Green=2/6, Cyan=3/6, Blue=4/6,
C             Magenta=5/6, Red=1, Yellow=7/6, Green=8/6, etc.
C             Default:  VREF=GMIN1, CREF=4/6 (blue).
C               where GMIN1 is the minimum defined grid value.
C               Default VMIN, VMAX, VPLUS, VSIGN, VCIRC, VREF and CREF
C               render the central value (GMIN1+GMAX1)/2 in yellow.
C             Hint: To plot the values oscillating around 0, you may
C               wish to set VREF=0 and CREF=0.166667 to render 0 in
C               yellow.
C     VCIRC4=real... If VCIRC is specified, VCIRC4 enables to modify it
C             in dependence on progressing time levels I4=1,2,...,N4:
C               VCIRC(I4)=VCIRC*EXP(-VCIRC4*(I4-1)*D4)
C             Default:  VCIRC4=0.
C (d) Saturation interpretation of values contained within file GRDPS2:
C     VSAT=real... The extent of values corresponding to the whole
C             saturation range from 0 (white) to 1 (fully saturated
C             colours).
C             Default:  VSAT=GMAX2-GMIN2 (or VSAT=1 for GMIN2=GMAX2).
C     VWHITE=real... Value corresponding to white (saturation=0).
C             Default:  VWHITE=GMIN2 (strictly speaking, =GMAX2-VSAT).
C     VSAT4=real... If VSAT is specified, VSAT4 enables to modify it in
C             dependence on progressing time levels I4=1,2,...,N4:
C               VSAT(I4)=VSAT*EXP(-VSAT4*(I4-1)*D4)
C             Default:  VSAT4=0.
C (e) Brightness interpretation of values contained within file GRDPS3:
C     VBRI=real... The extent of values corresponding to the whole
C             brightness range from 0 (black) to 1 (bright colours).
C             Default:  VBRI=GMAX3-GMIN3 (or VBRI=1 for GMIN3=GMAX3).
C     VBLACK=real... Value corresponding to black (brightness=0).
C             Default:  VBLACK=GMIN3 (strictly speaking, =GMAX3-VBRI).
C     VBRI4=real... If VBRI is specified, VBRI4 enables to modify it in
C             dependence on progressing time levels I4=1,2,...,N4:
C               VBRI(I4)=VBRI*EXP(-VBRI4*(I4-1)*D4)
C             Default:  VBRI4=0.
C (f) Form of the output PostScript file:
C     SHOWPAGE=integer... PostScript command 'showpage' at the end of
C             file directs printer to print and delete the picture.
C             This is usually what we want.  However, sometimes we may
C             wish to overlay the picture with another one.  In this
C             case, we wish to remove the 'showpage'.
C             SHOWPAGE=0: Two last lines of the file containing strings
C               'showpage' and '%%EOF' are not written.
C             SHOWPAGE=1: The lines are written.
C             Default:  SHOWPAGE=1
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER MGRID
      PARAMETER (MGRID=MRAM)
      REAL GRID(MGRID)
      EQUIVALENCE (GRID,RAM)
      INTEGER IGRID(MGRID)
      EQUIVALENCE (GRID,IGRID)
C
C.......................................................................
C
      CHARACTER*2 UNIT
      REAL UNITPT,HEIGHT,OFFSET,WIDTH
C
C     UNIT... One of: 'cm', 'in', 'pt'.
C     UNITPT...Size of the length unit, in which input data controlling
C             the size and position of the plot are expressed, in big
C             points (pt).  E.g., UNITPT=72./2.54 corresponds to plotting
C             in cm.
C     HEIGHT..Anticipated height of the paper sheet.
C     OFFSET..Left margin, and top or bottom margin for low or high
C             plots, respectively.
C     WIDTH...Default width of the plot.
C
      INTEGER MFILE,LU
      PARAMETER (MFILE=9,LU=10)
      CHARACTER*80 FSEP,FPS,FPSHEX,FILE(MFILE)
      CHARACTER*1 HEX1(0:15)
      CHARACTER*2 HEX2(0:255)
      CHARACTER*6 FGRDPS
      INTEGER LUIN(MFILE)
      REAL UNDEF
      PARAMETER (UNDEF=-999999.)
      INTEGER NFILE,IFILE,N1,N2,N3,N31,N32,N4,J,J0,I,I0,I1,I2,I31,I32,I4
      INTEGER ISHOW
      REAL XSIGN,YSIGN
      REAL VMIN,VMAX,VPLUS,VSIGN,VCIRC,VREF,CREF,ROTATE,RED,GREEN,BLUE
      REAL GMIN(MFILE),GMAX(MFILE),G,G0,GSTEP
      REAL BBOX1,BBOX2,BBOX3,BBOX4,BB1,BB2,BB3,BB4
      REAL BB1P,BB2P,BB3P,BB4P,BB2DEF,BB4DEF,AUX,VAUX,C,S
      DATA LUIN /1,2,3,4,5,6,7,8,9/
      DATA HEX1
     * /'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'/
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GRDPS: Enter input filename: '
      FSEP=' '
      READ(*,*) FSEP
      WRITE(*,'(A)') '+GRDPS: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.NE.' ') THEN
        CALL RSEP1(LU,FSEP)
      ELSE
C       GRDPS-04
        CALL ERROR('GRDPS-04: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
C     Reading input parameters from the SEP file:
      CALL RSEP3T('GRD',FILE(1),'grd.out')
      CALL RSEP3T('PS',FPS,'grd.ps')
      CALL RSEP3T('PSHEX',FPSHEX,' ')
      FGRDPS='GRDPS0'
      DO 12 IFILE=1,MFILE
        FGRDPS(6:6)=CHAR(ICHAR('0')+IFILE)
        IF (IFILE.NE.1) CALL RSEP3T(FGRDPS,FILE(IFILE),' ')
        IF(FILE(IFILE).EQ.' ') THEN
          NFILE=IFILE-1
          GO TO 13
        END IF
   12 CONTINUE
   13 CONTINUE
      CALL RSEP3I('SHOWPAGE',ISHOW,1)
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('N4',N4,1)
C     Transforming 2-D arrays to N1*N2*1 and 1-D arrays to N1*1*1
      IF(N1.EQ.1) THEN
        N1=N2
        N2=N3
        N3=1
      END IF
      IF(N1.EQ.1) THEN
        N1=N2
        N2=1
      END IF
      IF(N2.EQ.1) THEN
        N2=N3
        N3=1
      END IF
C     Determining the layout of 2-D slices of 3-D arrays
      CALL RSEP3I('NH',N31,INT(SQRT(FLOAT(N3))+.001))
      N32=(N3-1)/N31+1
      IF(N1*N2*N31*N32*NFILE.GT.MGRID) THEN
C       GRDPS-01
        CALL ERROR('GRDPS-01: Too small array RAM(MRAM)')
C       Array RAM(MRAM) allocated in include file 'ram.inc' is too small
C       to contain the grid values to be plotted.  You may wish to
C       increse the dimension MRAM in file 'ram.inc'.
C       ram.inc
      END IF
C     N3 slices will be arranged into N31*N32 slices.
C
C     Recalling the plotting unit and setting default dimensions:
      CALL RSEP3T('UNIT',UNIT,'cm')
      CALL LOWER(UNIT)
      IF(UNIT.EQ.'cm') THEN
        UNITPT=72./2.54
        HEIGHT=29.7
        OFFSET=2.5
        WIDTH=16.0
      ELSE IF(UNIT.EQ.'in') THEN
        UNITPT=72.
        HEIGHT=11.0
        OFFSET=1.0
        WIDTH=6.5
      ELSE IF(UNIT.EQ.'pt') THEN
        UNITPT=1.
        HEIGHT=FLOAT(N32*N2)
        OFFSET=0.0
        WIDTH=FLOAT(N31*N1)
      ELSE
C       GRDPS-02
        CALL ERROR('GRDPS-02: Unrecognized plotting units')
C       Allocated plotting units are UNIT='cm', UNIT='in' or UNIT='pt'.
      END IF

C.......................................................................
C
C     Loop over time slices:
      DO 90 I4=1,N4
        IF(I4.EQ.1) THEN
          IF(N4.GT.1) THEN
            DO 15 IFILE=1,NFILE
              OPEN(LUIN(IFILE),FILE=FILE(IFILE),
     *                                    FORM='FORMATTED',STATUS='OLD')
              FILE(IFILE)=' '
   15       CONTINUE
          END IF
        ELSE
C         Generating new output filenames:
          CALL NEWNAM(FPS)
          CALL NEWNAM(FPSHEX)
        END IF
C
C       Reading input grid values:
        DO 16 IFILE=1,NFILE
          CALL RARRAY(LUIN(IFILE),FILE(IFILE),'FORMATTED',.TRUE.,UNDEF,
     *                         N1*N2*N3,GRID(N1*N2*N31*N32*(IFILE-1)+1))
   16   CONTINUE
C
C       Rearranging 3-D input N1*N2*N3 grid into (N1*N31)*(N2*N32) grid:
        DO 29 IFILE=0,N1*N2*N31*N32*(NFILE-1),N1*N2*N31*N32
C         Extending the grid to N1*N2*N31*N32 points
          DO 21 I=IFILE+N1*N2*N3+1,IFILE+N1*N2*N31*N32
            GRID(I)=UNDEF
   21     CONTINUE
C         Transposing the N1*N2*N31*N32 grid to N1*N31*N2*N32 grid
          DO 25 I32=0,N2*N31*(N32-1),N2*N31
            DO 24 I0=0,N2*N31-1
              J0=I0
   22         CONTINUE
                I2=J0/N31
                I31=J0-I2*N31
                J0=I31*N2+I2
              IF(J0.LT.I0) GO TO 22
              I=IFILE+(I32+I0)*N1
              J=IFILE+(I32+J0)*N1
              DO 23 I1=1,N1
                I=I+1
                J=J+1
                AUX=GRID(I)
                GRID(I)=GRID(J)
                GRID(J)=AUX
   23         CONTINUE
   24       CONTINUE
   25     CONTINUE
   29   CONTINUE
C
        N1=N1*N31
        N2=N2*N32
C
C.......................................................................
C
C       Recalling the data for the plotting area:
        CALL RSEP3R('XSIGN'  ,XSIGN,1.)
        CALL RSEP3R('YSIGN'  ,YSIGN,1.)
        AUX=HEIGHT
        CALL RSEP3R('HEIGHT' ,HEIGHT,AUX)
        CALL RSEP3R('HSIZE'  ,BB3,SIGN(WIDTH,XSIGN))
        CALL RSEP3R('HOFFSET',BB1,OFFSET)
C       Default height of the figure
        BB4DEF=ABS(BB3)*FLOAT(N2)/FLOAT(N1)
        CALL RSEP3R('VSIZE'  ,BB4,SIGN(BB4DEF,YSIGN))
C       Default vertical position of the figure
        IF(ABS(BB4).LE.HEIGHT-2.*OFFSET) THEN
          BB2DEF=HEIGHT-OFFSET-ABS(BB4)
        ELSE IF(ABS(BB4).LE.HEIGHT) THEN
          BB2DEF=(HEIGHT-ABS(BB4))/2.
        ELSE
          BB2DEF=OFFSET
        END IF
        CALL RSEP3R('VOFFSET',BB2,BB2DEF)
        IF(BB3.LT.0.) THEN
          BB1=BB1-BB3
        END IF
        IF(BB4.LT.0.) THEN
          BB2=BB2-BB4
        END IF
        CALL RSEP3R('ROTATE',ROTATE,0.)
C
C       Transformation from plotting units (e.g. centimetres) to points:
        BB1P=BB1*UNITPT
        BB2P=BB2*UNITPT
        BB3P=BB3*UNITPT
        BB4P=BB4*UNITPT
C       Bounding box:
        BBOX1=AMIN1(BB1P,BB1P+BB3P)
        BBOX2=AMIN1(BB2P,BB2P+BB4P)
        BBOX3=AMAX1(BB1P,BB1P+BB3P)
        BBOX4=AMAX1(BB2P,BB2P+BB4P)
        IF(ROTATE.NE.0.) THEN
          C=COS(ROTATE*3.14159/180.)
          S=SIN(ROTATE*3.14159/180.)
          BBOX1=BBOX1-HEIGHT*UNITPT/2.
          BBOX2=BBOX2-HEIGHT*UNITPT/2.
          BBOX3=BBOX3-HEIGHT*UNITPT/2.
          BBOX4=BBOX4-HEIGHT*UNITPT/2.
          AUX  =C*BBOX1-S*BBOX2
          BBOX2=S*BBOX1+C*BBOX2
          BBOX1=AUX
          AUX  =C*BBOX3-S*BBOX4
          BBOX4=S*BBOX3+C*BBOX4
          BBOX3=AUX
          BBOX1=BBOX1+HEIGHT*UNITPT/2.
          BBOX2=BBOX2+HEIGHT*UNITPT/2.
          BBOX3=BBOX3+HEIGHT*UNITPT/2.
          BBOX4=BBOX4+HEIGHT*UNITPT/2.
          AUX  =AMIN1(BBOX1,BBOX3)
          BBOX3=AMAX1(BBOX1,BBOX3)
          BBOX1=AUX
          AUX  =AMIN1(BBOX2,BBOX4)
          BBOX4=AMAX1(BBOX2,BBOX4)
          BBOX2=AUX
        END IF
C
C.......................................................................
C
C       Calculating minimum and maximum values GMIN(*) and GMAX(*):
        CALL RSEP3R('VMIN',VMIN,-999999.)
        CALL RSEP3R('VMAX',VMAX, 999999000000.)
        DO 33 IFILE=1,NFILE
          NNN=N1*N2*(IFILE-1)
          GMINI=VMAX
          GMAXI=VMIN
          DO 31 I=NNN+1,NNN+N1*N2
            G=GRID(I)
            IF(G.NE.UNDEF) THEN
              IF(G.LE.VMIN.OR.VMAX.LE.G) THEN
                GRID(I)=UNDEF
              ELSE
                GMINI=AMIN1(GMINI,G)
                GMAXI=AMAX1(GMAXI,G)
              END IF
            END IF
   31     CONTINUE
C
C         Rescaling range GMIN--GMAX to 0--254
c         (undefined values are 255):
          IF(GMINI.NE.GMAXI) THEN
            GSTEP=(GMAXI-GMINI)/254.
          ELSE
            GSTEP=1.
          END IF
          G0=GMINI-GSTEP/2.
          DO 32 I=NNN+1,NNN+N1*N2
            IF(GRID(I).NE.UNDEF) THEN
              G=(GRID(I)-G0)/GSTEP
              IGRID(I)=INT(G)
            ELSE
              IGRID(I)=255
            END IF
   32     CONTINUE
C
          GMIN(IFILE)=GMINI
          GMAX(IFILE)=GMAXI
   33   CONTINUE
C
C       Preparing hexadecimal coding table:
        I=0
        DO 42 I2=0,15
          DO 41 I1=0,15
            HEX2(I)(1:1)=HEX1(I2)
            HEX2(I)(2:2)=HEX1(I1)
            I=I+1
   41     CONTINUE
   42   CONTINUE
C
C       Parameters controlling colours:
        CALL RSEP3R('D4'    ,D4    ,1.)
        TIME=FLOAT(I4-1)*D4
        CALL RSEP3R('VPLUS' ,VPLUS ,0.)
        CALL RSEP3R('VSIGN' ,VSIGN ,0.)
        CALL RSEP3R('VCIRC4',VAUX  ,0.)
        VAUX=EXP(VAUX*TIME)
        AUX=SIGN(GMAX(1)-GMIN(1)+VPLUS,VSIGN)*VAUX
        IF(AUX.EQ.0.) AUX=1.
        CALL RSEP3R('VCIRC' ,VCIRC ,AUX)
        VCIRC=VCIRC/VAUX
        CALL RSEP3R('VREF'  ,VREF  ,GMIN(1))
        CALL RSEP3R('CREF'  ,CREF  ,2./3.)
        IF(NFILE.GE.2) THEN
          CALL RSEP3R('VSAT4' ,VAUX  ,0.)
          VAUX=EXP(VAUX*TIME)
          AUX=(GMAX(2)-GMIN(2))*VAUX
          IF(AUX.EQ.0.) AUX=1.
          CALL RSEP3R('VSAT'  ,VSAT  ,AUX)
          VSAT=VSAT/VAUX
          CALL RSEP3R('VWHITE',VWHITE,GMAX(2)-AUX)
        END IF
        IF(NFILE.GE.3) THEN
          CALL RSEP3R('VBRI4' ,VAUX  ,0.)
          VAUX=EXP(VAUX*TIME)
          AUX=(GMAX(3)-GMIN(3))*VAUX
          IF(AUX.EQ.0.) AUX=1.
          CALL RSEP3R('VBRI'  ,VBRI  ,AUX)
          VBRI=VBRI/VAUX
          CALL RSEP3R('VBLACK',VBLACK,GMAX(3)-AUX)
        END IF
C
C       Writing PostScript prolog:
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Writing: ',FPS(1:MIN0(LEN(FPS),70))
        OPEN(LU,FILE=FPS)
        WRITE(LU,'(A/A,4I6,/(A))')
     *  '%!PS-Adobe-3.0',
     *  '%%BoundingBox:',INT(BBOX1+.5),INT(BBOX2+.5),
     *                   INT(BBOX3+.5),INT(BBOX4+.5),
     *  '%%EndComments',
     *  '%%BeginProlog',
     *  '%%BeginProcSet: (grdps)',
     *  '%%Creator: grdps',
     *  '%-----------------------------------------------------------',
     *  '% Default UNDEFINED procedure:'
C520    WRITE(LU,'(A)')
C520 *  '/UNDEFINED {0 0 0.80 sethsbcolor} bind def'
        CALL RSEP3R('R',RED  ,0.80)
        CALL RSEP3R('G',GREEN,0.80)
        CALL RSEP3R('B',BLUE ,0.80)
        WRITE(LU,'(A,3(F4.2,1X),A)')
     *  '/UNDEFINED {',RED,GREEN,BLUE,'setrgbcolor} bind def'
        WRITE(LU,'(A)')
     *  '% Default VALUEtoCOLOR procedure:'
        WRITE(LU,'(3(A,G13.6)/A,G13.6,A/A,G13.6,A/A)')
     *  '/VCIRC  ',VCIRC ,' def % Range of grid values:',
     *                                            GMIN(1),' to',GMAX(1),
     *  '/VREF   ',VREF  ,' def',
     *  '/CREF   ',CREF  ,' def',
     *  '/VRED VREF CREF VCIRC mul sub def'
        IF(NFILE.GE.2) THEN
          WRITE(LU,'(3(A,G13.6)/A,G13.6,A)')
     *    '/VSAT   ',VSAT  ,' def % Range of grid values:',
     *                                            GMIN(2),' to',GMAX(2),
     *    '/VWHITE ',VWHITE,' def'
        END IF
        IF(NFILE.GE.3) THEN
          WRITE(LU,'(3(A,G13.6)/A,G13.6,A)')
     *    '/VBRI   ',VBRI  ,' def % Range of grid values:',
     *                                            GMIN(3),' to',GMAX(3),
     *    '/VBLACK ',VBLACK,' def'
        END IF
        WRITE(LU,'(A)')
     *    '/VALUEtoCOLOR{'
        IF(NFILE.LE.1) THEN
          WRITE(LU,'(A)')
     *      ' VRED sub VCIRC div dup truncate sub dup 0 lt {1 add} if'
        ELSE
          WRITE(LU,'(A,I1,9A))')
     *      ' ',NFILE,' -1 roll',
     *      ' VRED sub VCIRC div dup truncate sub dup 0 lt {1 add} if'
        END IF
        IF(NFILE.LE.1) THEN
          WRITE(LU,'(A)')
     *    ' 1'
        ELSE
          WRITE(LU,'(A,I1,9A))')
     *      ' ',NFILE,' -1 roll',
     *      ' VWHITE sub VSAT div'
        END IF
        IF(NFILE.LE.2) THEN
          WRITE(LU,'(A)')
     *      ' 1'
        ELSE
          WRITE(LU,'(A,I1,9A))')
     *      ' ',NFILE,' -1 roll',
     *      ' VBLACK sub VBRI div'
        END IF
        DO 51 IFILE=4,NFILE
          WRITE(LU,'(A)')
     *    ' pop'
   51   CONTINUE
        WRITE(LU,'(A)')
     *  ' sethsbcolor} bind def',
     *  '%-----------------------------------------------------------',
     *  '% User-defined VALUEtoCOLOR procedure may be inserted here:'
        IF(FPSHEX.NE.' ') THEN
          CLOSE(LU)
          WRITE(*,'(''+'',79('' ''))')
          WRITE(*,'(2A)') '+Writing: ',FPSHEX(1:MIN0(LEN(FPSHEX),70))
          OPEN(LU,FILE=FPSHEX)
        END IF
C
C       Writing the plotting procedure:
        WRITE(LU,'(A)')
     *  '%-----------------------------------------------------------',
     *  '%%EndProcSet',
     *  '%%EndProlog',
     *  '%-----------------------------------------------------------',
     *  '%%BeginSetup',
     *  '% Numerical values describing the image size and position:'
        WRITE(LU,'(A,I6,A)')          '/N1',N1,' def'
        WRITE(LU,'(A,I6,A)')          '/N2',N2,' def'
        WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB1',BB1P,' def %',BB1,'cm'
        WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB2',BB2P,' def %',BB2,'cm'
        WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB3',BB3P,' def %',BB3,'cm'
        WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB4',BB4P,' def %',BB4,'cm'
        WRITE(LU,'(A,F8.1,A)')        '/PAPERSIZE',HEIGHT*UNITPT,' def'
        WRITE(LU,'(A,F8.1,A)')        '/ROTATE',ROTATE,' def'
        WRITE(LU,'(A)')
     *  '% Scaling of grid values:'
        DO 52 IFILE=1,NFILE
          WRITE(LU,'(A,I1,A,G13.6,A)')
     *    '/GMIN',IFILE,' ',GMIN(IFILE),' def'
          WRITE(LU,'(A,I1,A,G13.6,A)')
     *    '/GMAX',IFILE,' ',GMAX(IFILE),' def'
   52   CONTINUE
        WRITE(LU,'(A)')
     *  '%%EndSetup',
     *  '%-----------------------------------------------------------',
     *  '%%BeginObject: (grdps)',
     *  'PAPERSIZE  2 div dup translate ROTATE rotate',
     *  'PAPERSIZE -2 div dup translate',
     *  '/SCALE1 N1 BB3 div def',
     *  '/SCALE2 N2 BB4 div def'
        DO 53 IFILE=1,NFILE
          WRITE(LU,'(9(A,I1))')
     *    '/VSTEP',IFILE,' GMAX',IFILE,' GMIN',IFILE,' sub 254 div def'
   53   CONTINUE
        WRITE(LU,'(A)')
     *  '/RGB 3 string def',
     *  'N1 N2 8 [ SCALE1 0 0 SCALE2 BB1 SCALE1 mul neg',
     *  '                            BB2 SCALE2 mul neg ]'
        WRITE(LU,'(A,I1,A)')
     *  '{currentfile ',NFILE,' string readhexstring pop'
        DO 54 IFILE=1,NFILE
          WRITE(LU,'(9(A,I1))')
     *    ' dup ',IFILE-1,
     *                  ' get VSTEP',IFILE,' mul GMIN',IFILE,' add exch'
   54   CONTINUE
        WRITE(LU,'(11A)')
     *  ' 0 get 255 eq {',('pop ',IFILE=1,NFILE),
     *                                'UNDEFINED} {VALUEtoCOLOR} ifelse'
        WRITE(LU,'(A)')
     *  ' currentrgbcolor',
     *  ' 255 mul .5 add cvi RGB exch 2 exch put',
     *  ' 255 mul .5 add cvi RGB exch 1 exch put',
     *  ' 255 mul .5 add cvi RGB exch 0 exch put RGB}',
     *  '                            bind false 3 colorimage'
C
C       Writing output hexadecimal values:
        N12=N1*N2
        NNN=N1*N2*NFILE
        IF(NFILE.EQ.1) THEN
          WRITE(LU,'(40A2)') (HEX2(IGRID(I)),I=1,N12)
        ELSE
          DO 61 I=N12+1,NNN
            IF(IGRID(I).GE.255) THEN
              IGRID(MOD(I,N12))=255
            END IF
   61     CONTINUE
          WRITE(LU,'(40A2)')
     *                    ((HEX2(IGRID(IFILE)),IFILE=I,NNN,N12),I=1,N12)
        END IF
C
C       Writing PostScript trailer:
        WRITE(LU,'(A)')
     *  'PAPERSIZE  2 div dup translate ROTATE neg rotate',
     *  'PAPERSIZE -2 div dup translate',
     *  '%%EndObject'
        IF(ISHOW.NE.0) THEN
          WRITE(LU,'(A)')
     *    'showpage',
     *    '%%EOF'
        END IF
        CLOSE(LU)
C
   90 CONTINUE
C     End of the loop over time slices.
C
      IF(N4.GT.1) THEN
        DO 91 IFILE=1,NFILE
          CLOSE(LUIN(IFILE))
   91   CONTINUE
      END IF
      WRITE(*,'(A)') '+GRDPS: Done.                 '
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE NEWNAM(NAME)
      CHARACTER*(*) NAME
C
C-----------------------------------------------------------------------
C
      INTEGER N,I
C
      IF(NAME.NE.' ') THEN
        N=LEN(NAME)
   10   CONTINUE
        DO 11 I=N,1,-1
          IF(LLE('0',NAME(I:I)).AND.LLE(NAME(I:I),'8')) THEN
            NAME(I:I)=CHAR(ICHAR(NAME(I:I))+1)
            GO TO 12
          ELSE IF(NAME(I:I).EQ.'9') THEN
            NAME(I:I)='0'
            N=I-1
            GO TO 10
          END IF
   11   CONTINUE
C         GRDPS-03
          CALL ERROR('GRDPS-03: Too many output files')
C         The digits in the template name of the output files do not
C         allow for the generation of all output filenames.
C         The number of digits should be increased.
   12   CONTINUE
      END IF
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C