C 
C Program 'GRDCAL' (GRiD CALculator) to perform vectorial calculations
C with real-valued arrays stored in disk files.
C
C Version: 5.20
C Date: 1998, September 24
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                                                    
C Description of the data files:
C
C The data are read in by the list directed input (free format).
C In the description of data files, each numbered paragraph indicates
C the beginning of a new input operation (new READ statement).
C If the symbolic name of the input variable is enclosed in apostrophes,
C the corresponding value in input data is of the type CHARACTER, i.e.
C it should be a character string enclosed in apostrophes.  If the first
C letter of the symbolic name is I-N, the corresponding value is of the
C type INTEGER.  Otherwise, the input parameter is of the type REAL and
C may or may not contain a decimal point.
C
C Input data read from the * external unit:
C     The interactive * external unit may also be redirected to the file
C     containing the relevant data.
C (1) 'SEP','CAL','GRD1','GRD2',...,'GRDn',/
C     'SEP'...String in apostrophes containing the name of the input
C             file with the data specifying grid dimensions.
C             Description of file SEP
C     'CAL'...String in apostrophes containing the name of the input
C             file containing the commands to be carried out at each
C             gridpoint.
C             Description of file CAL
C     'GRD1','GRD2',...,'GRDn'... Strings in apostrophes containing the
C             names of the input/output ASCII files with the grid
C             values.
C     /...    Input data line must be terminated by a slash.
C     Default: 'SEP'='grd.h', 'CAL'=' ', 'GRD1'=' ', ..., 'GRDn'=' '.
C
C                                                     
C Data file SEP has the form of the SEP (Stanford Exploration Project)
C parameter file:
C     All the data are specified in the form of PARAMETER=VALUE, e.g.
C     N1=50, with PARAMETER directly preceding = without intervening
C     spaces and with VALUE directly following = without intervening
C     spaces.  The PARAMETER=VALUE couple must be delimited by a space
C     or comma from both sides.
C     The PARAMETER string is not case-sensitive.
C     PARAMETER= followed by a space resets the default parameter value.
C     All other text in the input files is ignored.  The file thus may
C     contain unused data or comments without leading comment character.
C     Everything between comment character # and the end of the
C     respective line is ignored, too.
C     The PARAMETER=VALUE couples may be specified in any order.
C     The last appearance takes precedence.
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             Default: N4=1
C
C                                                     
C File CAL containing commands to be performed at each gridpoint:
C     Each line may contain at most a single command.
C     The lines are read character-by-character.  The commands thus
C     should not be enclosed in parentheses.  The commands have the
C     structure like:
C       $3=$1+$2
C       C=A-B
C       C=$1-$2
C       $3=ABS(C)
C     or
C       A=SQRT($2)
C       A=$1*A
C       @4=@3*A
C     etc.
C     Here $i or @i corresponds to the i-th input/output file GRDi,
C     FUN(.) represents function FUN of a single argument,
C     FUN(.,.) represents function FUN of two arguments,
C     other names represent temporary variables.
C     Letter case is not distinguished.
C     A single line may contain a single operation.
C
C     Input and output files with data cubes:
C       $1,$2,$3,...,$9 are space data cubes of N1*N2*N3 grid values.
C             If N4 time levels are processed, the input space data
C             cubes apply to all time levels and output space data cubes
C             correspond to the first time level.
C             Examples: gridded P or S wave velocities or density.
C       @1,@2,@3,...,@9 are space-time data cubes of N1*N2*N3*N4 grid
C             values, read, processed and written by time levels.
C             Examples: components of wavefield snapshots.
C       The same input file cannot be denoted both by $i and @i.
C       The same output file cannot be denoted both by $i and @i.
C       If N4=1, there is no difference between $i and @i.
C       If an input or output $ file is considered, the @ files are read
C       and written by time levels.
C       If all input or all output @ files do not fit into array RAM,
C       the @ files are read and written by time levels.
C       If N4.GT.1, the @ files are read and written by time levels and
C       file @i is used for input, neither file $i nor @i can be used
C       for output.
C     Temporary variables:
C       If the value of a variable is not defined within the command
C       file, it is taken from the input SEP file with the zero default.
C       If the value of a variable is defined within the command file,
C       it must be defined before it is used on the right-hand side of
C       any command.
C     Allowed operators:
C       A=B+C
C       A=B-C
C       A=B*C
C       A=B/C
C       A=B**C
C     Allowed functions (= sign means equivalent function names):
C       ABS(.)
C       AINT(.)=INT(.)
C       ANINT(.)=NINT(.)
C       AMOD(.)=MOD(.)
C       SIGN(.)
C       DIM(.)
C       AMAX1(.,.)=AMAX(.,.)=MAX(.,.)
C       AMIN1(.,.)=AMIN(.,.)=MIN(.,.)
C       SQRT(.)
C       EXP(.)
C       ALOG(.)=LOG(.)=LN(.)
C       ALOG10(.)=LOG10(.)
C       SIN(.)
C       COS(.)
C       TAN(.)
C       ASIN(.)
C       ACOS(.)
C       ATAN(.)
C       ATAN2(.,.)
C       SINH(.)
C       COSH(.)
C       TANH(.)
C     Note that commands of forms like
C       $2=INT(A)
C       $3=NINT(A)
C     imply integer output values, which are thus written as integers
C     to the output ASCII files.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working array:
      INTEGER IRAM(MRAM)
      REAL GRID(MRAM)
      EQUIVALENCE (IRAM,RAM),(GRID,RAM)
C
C.......................................................................
C
      INTEGER MFILE,MNAME,MKOM
C     
      PARAMETER (MFILE=9,MNAME=2*MFILE+20,MKOM=100,LU1=1)
      CHARACTER*80 FGRD,FKOM,FILE(MFILE)
      INTEGER KGRID0(MFILE),KGRID1(MFILE),LU(MFILE)
      CHARACTER*80 NAME(MNAME)
      REAL RNAME(MNAME)
      INTEGER KOM0(MKOM),KOM1(MKOM),KOM2(MKOM),KOM3(MKOM)
C     KGRID0..Offsets in RAM for storing the output grids.
C     KGRID1..Offsets in RAM for storing the input grids.
C     NNAME...Number of all operands and results in the command file.
C             The first MFILE positions are reserved to files.
C     NAME... Strings identifying the operands and results.
C     RNAME...Values of the operands and results.
C
      CHARACTER*255 LINE
      CHARACTER*7 FORMAT
      LOGICAL LUNDEF,LARRAY
      INTEGER NNAME,NKOM
C     LARRAY..Logical variable identifying whether grid values must be
C             split into individual time levels to fit in the RAM.
C     NKOM... Number of commands
C
      DATA LU/1,2,3,4,5,6,7,8,9/
C
C.......................................................................
C
C     Main input data:
C     Default:
      FGRD='grd.h'
      FKOM=' '
      DO 10 IFILE=1,MFILE
        FILE(IFILE)=' '
        NAME(IFILE)='$'
        NAME(IFILE)(2:2)=CHAR(ICHAR('0')+IFILE)
        NAME(MFILE+IFILE)='@'
        NAME(MFILE+IFILE)(2:2)=CHAR(ICHAR('0')+IFILE)
   10 CONTINUE
C     Reading main input data:
      WRITE(*,'(A)')
     *    ' Enter filenames of Grid header, Commands, In/Out grids, /: '
      WRITE(*,'(A)')
     *    '+GRDCAL: Working...                                         '
      READ (*,*) FGRD,FKOM,FILE
C     Default extension of FKOM is '.cal':
      IF(INDEX(FKOM,'.').EQ.0) THEN
        FKOM(LENGTH(FKOM)+1:LENGTH(FKOM)+4)='.cal'
      END IF
C
C     Reading all the data from file FGRD to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU1,FGRD)

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
C.......................................................................
C
C     Reading the command file FKOM:
C
      NKOM=0
      NNAME=2*MFILE
      OPEN(LU1,FILE=FKOM,STATUS='OLD')
C
C     Loop over input lines
   11 CONTINUE
        READ(LU1,'(A)',END=19) LINE
        KEQ=INDEX(LINE,'=')
        IF(KEQ.NE.0) THEN
C
C         The line contains a new command
          NKOM=NKOM+1
          IF(NKOM.GT.MKOM) THEN
C           GRDCAL-01
            CALL ERROR('GRDCAL-01: Insufficient memory for commands')
C           Maximum number MKOM of the commands read from the command
C           file should probably be increased.  MKOM is declared by the
C           PARAMETER statement.
          END IF
          CALL LOWER(LINE)
C
C         Name of the result must precede '=':
          DO 12 K=KEQ-1,1,-1
            IF(LINE(K:K).EQ.' ') THEN
              GO TO 13
            END IF
   12     CONTINUE
   13     CONTINUE
          IF(K.GE.KEQ-1) THEN
C           GRDCAL-02
            CALL ERROR('GRDCAL-02: Missing identifier of the result')
          END IF
C         Registration of the name
          CALL REGNAM(LINE(K+1:KEQ-1),NAME,MNAME,NNAME,KOM0(NKOM))
C
C         End of the command:
          KEND=INDEX(LINE(KEQ+1:),' ')
          IF(KEND.EQ.0) THEN
C           GRDCAL-03
            CALL ERROR('GRDCAL-03: Too long command line')
          END IF
          IF(KEND.EQ.1) THEN
C           GRDCAL-39
            CALL ERROR('GRDCAL-39: Command line has no right-hand side')
          END IF
          KEND=KEQ+KEND-1
C
C         Search for left parenthesis:
          K=INDEX(LINE(KEQ+1:KEND),'(')
          IF(K.EQ.0) THEN
C
C           No left parenthesis - search for binary operators:
            K=INDEX(LINE(KEQ+1:KEND),'**')
            IF(K.NE.0) THEN
C             Two-letter binary operator **:
              KOM3(NKOM)=5
C             Registration of the name of the second operand
              K=KEQ+K
              CALL REGNAM(LINE(K+2:KEND-1),NAME,MNAME,NNAME,KOM2(NKOM))
            ELSE
C             Search for a one-letter binary operator:
              K=INDEX(LINE(KEQ+2:KEND+1),'+')
              IF(K.NE.0) THEN
                K=K+1
                KOM3(NKOM)=1
              ELSE
                K=INDEX(LINE(KEQ+2:KEND+1),'-')
                IF(K.NE.0) THEN
                  K=K+1
                  KOM3(NKOM)=2
                ELSE
                  K=INDEX(LINE(KEQ+1:KEND),'*')
                  IF(K.NE.0) THEN
                    KOM3(NKOM)=3
                IF(K.EQ.1) THEN
                  K=0
                  KOM3(NKOM)=0
                END IF
                  ELSE
                    K=INDEX(LINE(KEQ+1:KEND),'/')
                    IF(K.NE.0) THEN
                      KOM3(NKOM)=4
                    ELSE
C                     No binary operator:
                      KOM3(NKOM)=0
                    END IF
                  END IF
                END IF
              END IF
              K=KEQ+K
              IF(KOM3(NKOM).NE.0) THEN
C               Registration of the name of the second operand
                CALL REGNAM
     *                    (LINE(K+1:KEND),NAME,MNAME,NNAME,KOM2(NKOM))
                IF(K+1.GT.KEND) THEN
C                 
C                 GRDCAL-04
                  WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
                  CALL ERROR('GRDCAL-04: Missing second operand')
                END IF
              END IF
            END IF
            IF(KOM3(NKOM).NE.0) THEN
C             Registration of the name of the first operand
              IF(KEQ+1.GT.K-1) THEN
C               GRDCAL-05
                WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
                CALL ERROR('GRDCAL-05: Missing first operand')
              END IF
              CALL REGNAM(LINE(KEQ+1:K-1),NAME,MNAME,NNAME,KOM1(NKOM))
            ELSE
C             Registration of the name of the single operand
              IF(KEQ+1.GT.KEND) THEN
C               GRDCAL-06
                WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
                CALL ERROR('GRDCAL-06: Missing operand')
              END IF
              CALL REGNAM
     *                  (LINE(KEQ+1:KEND),NAME,MNAME,NNAME,KOM1(NKOM))
              KOM2(NKOM)=0
            END IF
C
          ELSE
C
C           Operator has the form of Fortran 77 intrinsic function
            K=KEQ+K
            IF(LINE(KEND:KEND).NE.')') THEN
C             GRDCAL-07
              WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
              CALL ERROR('GRDCAL-07: Missing closing parenthesis')
            END IF
C           Search for comma delimiting the arguments
            I=INDEX(LINE(K+1:KEND-1),',')
C           Registration of the arguments
            IF(I.EQ.0) THEN
C             Single argument:
              IF(K+1.GT.KEND-1) THEN
C               GRDCAL-08
                WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
                CALL ERROR('GRDCAL-08: Missing argument')
              END IF
              CALL REGNAM(LINE(K+1:KEND-1),NAME,MNAME,NNAME,KOM1(NKOM))
              KOM2(NKOM)=0
            ELSE
C             Two arguments:
              I=K+I
              IF(K+1.GT.I-1) THEN
C               GRDCAL-09
                WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
                CALL ERROR('GRDCAL-09: Missing first argument')
              END IF
              CALL REGNAM(LINE(K+1:I-1),NAME,MNAME,NNAME,KOM1(NKOM))
              IF(I+1.GT.KEND-1) THEN
C               GRDCAL-10
                WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
                CALL ERROR('GRDCAL-10: Missing second argument')
              END IF
              CALL REGNAM(LINE(I+1:KEND-1),NAME,MNAME,NNAME,KOM2(NKOM))
            END IF
C           Registration of the function
            IF(LINE(KEQ+1:K-1).EQ.'abs') THEN
              KOM3(NKOM)= 6
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-11
                CALL ERROR('GRDCAL-11: Redundant argument in ABS')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'aint'
     *          .OR.LINE(KEQ+1:K-1).EQ.'int') THEN
              KOM3(NKOM)= 7
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-12
                CALL ERROR('GRDCAL-12: Redundant argument in AINT')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'anint'
     *          .OR.LINE(KEQ+1:K-1).EQ.'nint') THEN
              KOM3(NKOM)= 8
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-13
                CALL ERROR('GRDCAL-13: Redundant argument in ANINT')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'amod'
     *          .OR.LINE(KEQ+1:K-1).EQ.'mod') THEN
              KOM3(NKOM)= 9
              IF(KOM2(NKOM).EQ.0) THEN
C               GRDCAL-14
                CALL ERROR('GRDCAL-14: Missing second argument of AMOD')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'sign') THEN
              KOM3(NKOM)=10
              IF(KOM2(NKOM).EQ.0) THEN
C               GRDCAL-15
                CALL ERROR('GRDCAL-15: Missing second argument of SIGN')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'dim') THEN
              KOM3(NKOM)=11
              IF(KOM2(NKOM).EQ.0) THEN
C               GRDCAL-16
                CALL ERROR('GRDCAL-16: Missing second argument of DIM')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'amax1'
     *          .OR.LINE(KEQ+1:K-1).EQ.'amax'
     *          .OR.LINE(KEQ+1:K-1).EQ.'max') THEN
              KOM3(NKOM)=12
              IF(KOM2(NKOM).EQ.0) THEN
C               GRDCAL-17
                CALL ERROR
     *                   ('GRDCAL-17: Missing second argument of AMAX1')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'amin1'
     *          .OR.LINE(KEQ+1:K-1).EQ.'amin'
     *          .OR.LINE(KEQ+1:K-1).EQ.'min') THEN
              KOM3(NKOM)=13
              IF(KOM2(NKOM).EQ.0) THEN
C               GRDCAL-18
                CALL ERROR
     *                   ('GRDCAL-18: Missing second argument of AMIN1')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'sqrt') THEN
              KOM3(NKOM)=14
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-19
                CALL ERROR('GRDCAL-19: Redundant argument in SQRT')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'exp') THEN
              KOM3(NKOM)=15
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-20
                CALL ERROR('GRDCAL-20: Redundant argument in EXP')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'alog'
     *          .OR.LINE(KEQ+1:K-1).EQ.'log'
     *          .OR.LINE(KEQ+1:K-1).EQ.'ln') THEN
              KOM3(NKOM)=16
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-21
                CALL ERROR('GRDCAL-21: Redundant argument in ALOG')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'alog10'
     *          .OR.LINE(KEQ+1:K-1).EQ.'log10') THEN
              KOM3(NKOM)=17
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-22
                CALL ERROR('GRDCAL-22: Redundant argument in ALOG10')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'sin') THEN
              KOM3(NKOM)=18
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-23
                CALL ERROR('GRDCAL-23: Redundant argument in SIN')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'cos') THEN
              KOM3(NKOM)=19
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-24
                CALL ERROR('GRDCAL-24: Redundant argument in COS')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'tan') THEN
              KOM3(NKOM)=20
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-25
                CALL ERROR('GRDCAL-25: Redundant argument in TAN')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'asin') THEN
              KOM3(NKOM)=21
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-26
                CALL ERROR('GRDCAL-26: Redundant argument in ASIN')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'acos') THEN
              KOM3(NKOM)=22
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-27
                CALL ERROR('GRDCAL-27: Redundant argument in ACOS')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'atan') THEN
              KOM3(NKOM)=23
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-28
                CALL ERROR('GRDCAL-28: Redundant argument in ATAN')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'atan2') THEN
              KOM3(NKOM)=24
              IF(KOM2(NKOM).EQ.0) THEN
C               GRDCAL-29
                CALL ERROR
     *                   ('GRDCAL-29: Missing second argument of ATAN2')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'sinh') THEN
              KOM3(NKOM)=25
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-30
                CALL ERROR('GRDCAL-30: Redundant argument in SINH')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'cosh') THEN
              KOM3(NKOM)=26
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-31
                CALL ERROR('GRDCAL-31: Redundant argument in COSH')
              END IF
            ELSE IF(LINE(KEQ+1:K-1).EQ.'tanh') THEN
              KOM3(NKOM)=27
              IF(KOM2(NKOM).NE.0) THEN
C               GRDCAL-32
                CALL ERROR('GRDCAL-32: Redundant argument in TANH')
              END IF
            ELSE
C             GRDCAL-33
              WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
              CALL ERROR     ('GRDCAL-33: Unknown function')
            END IF
C
          END IF
        END IF
      GO TO 11
   19 CONTINUE
      CLOSE(LU1)
C
C     Interpreting the constants:
      FORMAT='(F00.0)'
      DO 20 I=1,NNAME
        IF(('0'.LE.NAME(I)(1:1).AND.NAME(I)(1:1).LE.'9').OR.
     *       NAME(I)(1:1).EQ.'+'.OR.NAME(I)(1:1).EQ.'-'.OR.
     *       NAME(I)(1:1).EQ.'.') THEN
          L=LENGTH(NAME(I))
          FORMAT(3:3)=CHAR(ICHAR('0')+L/10)
          FORMAT(4:4)=CHAR(ICHAR('0')+MOD(L,10))
          READ(NAME(I),FORMAT) RNAME(I)
        ELSE
          CALL RSEP3R(NAME(I),RNAME(I),0.)
        END IF
   20 CONTINUE
C
C.......................................................................
C
C     Logical variable identifying whether grid values must be split
C     into individual time levels to fit in the RAM:
      LARRAY=.FALSE.
C
      DO 21 IFILE=1,MFILE
        KGRID0(IFILE)=-1
        KGRID1(IFILE)=-1
   21 CONTINUE
C
C     Determining storage for input grid values:
      IGRID=0
      DO 33 JFILE=1,2*MFILE
        IFILE=MOD(JFILE-1,MFILE)+1
        DO 31 IKOM=1,NKOM
          IF(KOM1(IKOM).EQ.JFILE.OR.KOM2(IKOM).EQ.JFILE) THEN
C           File appears at the R.H.S. of the command:
            IF(FILE(IFILE).EQ.' ') THEN
C             GRDCAL-34
              CALL ERROR('GRDCAL-34: Blank filename of input grid')
            END IF
            IF(JFILE.LE.MFILE) THEN
C             Space grid ($) on input,
C             calculation must be split into individual time levels
              IF(N4.GT.1) THEN
                LARRAY=.TRUE.
              END IF
            ELSE
              IF(KGRID1(IFILE).GE.0) THEN
C               GRDCAL-41
                CALL ERROR('GRDCAL-41: $ and @ for the same input file')
C               Coinciding input space ($) and space-time (@) data cubes
              END IF
            END IF
            KGRID1(IFILE)=IGRID
            IGRID=IGRID+N1*N2*N3
            GO TO 32
          END IF
   31   CONTINUE
   32   CONTINUE
        IF(JFILE.EQ.MFILE) THEN
C         The part of RAM reserved for input spatial grids ($-files):
          JGRID=IGRID
        END IF
   33 CONTINUE
      IF(IGRID.GT.MRAM) THEN
C       GRDCAL-35
        CALL ERROR('GRDCAL-35: Insufficient memory for input grids')
C       Dimension MRAM of array RAM in include file
C       ram.inc should probably be increased
C       to accommodate all input grids.
      END IF
      IF(IGRID*N4.GT.MRAM) THEN
C       Grid values must be split into individual time levels
        LARRAY=.TRUE.
      END IF
C
C     Determining storage for output grid values:
      IF(N4.LE.1) THEN
        IGRID=0
      ELSE
C       Protecting the part of RAM with input spatial grids ($-files):
        IGRID=JGRID
      END IF
      DO 43 JFILE=1,2*MFILE
        IFILE=MOD(JFILE-1,MFILE)+1
        DO 41 IKOM=1,NKOM
          IF(KOM0(IKOM).EQ.JFILE) THEN
C           File appears at the L.H.S. of the command:
            IF(FILE(IFILE).EQ.' ') THEN
C             GRDCAL-36
              CALL ERROR('GRDCAL-36: Blank filename of output grid')
            END IF
            IF(JFILE.LE.MFILE) THEN
C             Space grid ($) on output,
C             calculation must be split into individual time levels
              IF(N4.GT.1) THEN
                LARRAY=.TRUE.
              END IF
            ELSE
              IF(KGRID0(IFILE).GE.0) THEN
C               GRDCAL-42
                CALL ERROR
     *                   ('GRDCAL-42: $ and @ for the same output file')
C               Coinciding output space ($) and space-time (@) data
C               cubes.
              END IF
            END IF
            KGRID0(IFILE)=IGRID
            IGRID=IGRID+N1*N2*N3
            GO TO 42
          END IF
   41   CONTINUE
   42   CONTINUE
   43 CONTINUE
      IF(IGRID.GT.MRAM) THEN
C       GRDCAL-37
        CALL ERROR('GRDCAL-37: Insufficient memory for output grids')
C       Dimension MRAM of array RAM in include file
C       ram.inc should probably be increased
C       to accommodate all output grids.
      END IF
      IF(IGRID*N4.GT.MRAM) THEN
C       Grid values must be split into individual time levels
        LARRAY=.TRUE.
      END IF
      IF(LARRAY) THEN
        NGRID=N1*N2*N3
        NTIME=N4
      ELSE
        NGRID=N1*N2*N3*N4
        NTIME=1
        DO 44 IFILE=1,MFILE
          KGRID0(IFILE)=KGRID0(IFILE)*N4
          KGRID1(IFILE)=KGRID1(IFILE)*N4
   44   CONTINUE
      END IF
C
C     Loop over time slices:
      DO 900 ITIME=1,NTIME
C
C     Reading input grid values:
C     3-D space grids ($)
      IF(ITIME.EQ.1) THEN
        DO 52 IFILE=1,MFILE
          DO 51 IKOM=1,NKOM
            IF(KOM1(IKOM).EQ.IFILE.OR.KOM2(IKOM).EQ.IFILE) THEN
C             File appears at the R.H.S. of the command:
              CALL RARRAY(LU(IFILE),FILE(IFILE),'FORMATTED',.TRUE.,
     *                       -999999.,N1*N2*N3,   GRID(KGRID1(IFILE)+1))
            END IF
   51     CONTINUE
   52   CONTINUE
      END IF
C     4-D space-time grids (@)
      DO 58 JFILE=MFILE+1,2*MFILE
        IFILE=MOD(JFILE-1,MFILE)+1
        DO 56 IKOM=1,NKOM
          IF(KOM1(IKOM).EQ.JFILE.OR.KOM2(IKOM).EQ.JFILE) THEN
C           File appears at the R.H.S. of the command:
            IF(LARRAY) THEN
              DO 55 I=1,NKOM
                IF(KOM0(I).EQ.IFILE.OR.KOM0(I).EQ.JFILE) THEN
C                 
C                 GRDCAL-40
                  CALL ERROR('GRDCAL-40: Same input and output files')
C                 If the grid is as huge as it must be stored in RAM
C                 by individual time slices, the output filenames cannot
C                 coincide with input filenames.
                END IF
   55         CONTINUE
              IF(ITIME.EQ.1) THEN
                OPEN(LU(IFILE),FILE=FILE(IFILE),FORM='FORMATTED')
              END IF
              CALL RARRAY(LU(IFILE),' ' ,'FORMATTED',.TRUE.,-999999.,
     *                    N1*N2*N3,   GRID(KGRID1(IFILE)+1))
              IF(ITIME.EQ.NTIME) THEN
                CLOSE(LU(IFILE))
              END IF
            ELSE
              CALL RARAY(LU1,FILE(IFILE),'FORMATTED',.TRUE.,-999999.,
     *                    N1*N2*N3,N4,GRID(KGRID1(IFILE)+1))
            END IF
            GO TO 57
          END IF
   56   CONTINUE
   57   CONTINUE
   58 CONTINUE
C
C.......................................................................
C
C     Performing grid calculations:
      IF(ITIME.EQ.1) THEN
        IF(LARRAY) THEN
          WRITE(*,'(A)')
     *    '+GRDCAL: Reading, calculating, writing...                   '
        ELSE
          WRITE(*,'(A)')
     *    '+GRDCAL: Calculating...                                     '
        END IF
      END IF
C
C     Loop for gridpoints:
      DO 202 IGRID=1,NGRID
C
C       Loop for individual commands:
        DO 201 IKOM=1,NKOM
          I0=KOM0(IKOM)
          I1=KOM1(IKOM)
          I2=KOM2(IKOM)
          LUNDEF=.FALSE.
          IF(I1.LE.2*MFILE) THEN
            IF(I1.LE.MFILE) THEN
              RNAME(I1)=GRID(KGRID1(I1)+IGRID)
            ELSE
              RNAME(I1)=GRID(KGRID1(I1-MFILE)+IGRID)
            END IF
          END IF
          IF(RNAME(I1).LT.-999998.) THEN
            LUNDEF=.TRUE.
          END IF
          IF(I2.GT.0) THEN
            IF(I2.LE.2*MFILE) THEN
              IF(I2.LE.MFILE) THEN
                RNAME(I2)=GRID(KGRID1(I2)+IGRID)
              ELSE
                RNAME(I2)=GRID(KGRID1(I2-MFILE)+IGRID)
              END IF
            END IF
            IF(RNAME(I2).LT.-999998.) THEN
              LUNDEF=.TRUE.
            END IF
          END IF
          IF(LUNDEF) THEN
            RNAME(I0)=-999999.
          ELSE
C
            GO TO (101,102,103,104,105,106,107,108,109,110,
     *             111,112,113,114,115,116,117,118,119,120,
     *             121,122,123,124,125,126,127) KOM3(IKOM)
              RNAME(I0)=RNAME(I1)
              GO TO 199
  101       CONTINUE
              RNAME(I0)=RNAME(I1)+RNAME(I2)
              GO TO 199
  102       CONTINUE
              RNAME(I0)=RNAME(I1)-RNAME(I2)
              GO TO 199
  103       CONTINUE
              RNAME(I0)=RNAME(I1)*RNAME(I2)
              GO TO 199
  104       CONTINUE
              IF(RNAME(I2).EQ.0.) THEN
                IF(RNAME(I1).EQ.0.) THEN
                  RNAME(I0)=0.
                ELSE
                  RNAME(I0)=-999999.
                END IF
              ELSE
                RNAME(I0)=RNAME(I1)/RNAME(I2)
              END IF
              GO TO 199
  105       CONTINUE
              IF(RNAME(I1).LT.0.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=RNAME(I1)**RNAME(I2)
              END IF
              GO TO 199
  106       CONTINUE
              RNAME(I0)=ABS(RNAME(I1))
              GO TO 199
  107       CONTINUE
              RNAME(I0)=AINT(RNAME(I1))
              GO TO 199
  108       CONTINUE
              RNAME(I0)=ANINT(RNAME(I1))
              GO TO 199
  109       CONTINUE
              IF(RNAME(I2).EQ.0.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=AMOD(RNAME(I1),RNAME(I2))
              END IF
              GO TO 199
  110       CONTINUE
              RNAME(I0)=SIGN(RNAME(I1),RNAME(I2))
              GO TO 199
  111       CONTINUE
              RNAME(I0)=DIM(RNAME(I1),RNAME(I2))
              GO TO 199
  112       CONTINUE
              RNAME(I0)=AMAX1(RNAME(I1),RNAME(I2))
              GO TO 199
  113       CONTINUE
              RNAME(I0)=AMIN1(RNAME(I1),RNAME(I2))
              GO TO 199
  114       CONTINUE
              IF(RNAME(I1).LT.0.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=SQRT(RNAME(I1))
              END IF
              GO TO 199
  115       CONTINUE
              RNAME(I0)=EXP(RNAME(I1))
              GO TO 199
  116       CONTINUE
              IF(RNAME(I1).LE.0.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=ALOG(RNAME(I1))
              END IF
              GO TO 199
  117       CONTINUE
              IF(RNAME(I1).LE.0.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=ALOG10(RNAME(I1))
              END IF
              GO TO 199
  118       CONTINUE
              RNAME(I0)=SIN(RNAME(I1))
              GO TO 199
  119       CONTINUE
              RNAME(I0)=COS(RNAME(I1))
              GO TO 199
  120       CONTINUE
              RNAME(I0)=TAN(RNAME(I1))
              GO TO 199
  121       CONTINUE
              IF(ABS(RNAME(I1)).GT.1.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=ASIN(RNAME(I1))
              END IF
              GO TO 199
  122       CONTINUE
              IF(ABS(RNAME(I1)).GT.1.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=ACOS(RNAME(I1))
              END IF
              GO TO 199
  123       CONTINUE
              IF(ABS(RNAME(I1)).GT.1.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=ATAN(RNAME(I1))
              END IF
              GO TO 199
  124       CONTINUE
              IF(RNAME(I1).EQ.0..AND.RNAME(I2).EQ.0.) THEN
                RNAME(I0)=-999999.
              ELSE
                RNAME(I0)=ATAN2(RNAME(I1),RNAME(I2))
              END IF
              GO TO 199
  125       CONTINUE
              RNAME(I0)=SINH(RNAME(I1))
              GO TO 199
  126       CONTINUE
              RNAME(I0)=COSH(RNAME(I1))
              GO TO 199
  127       CONTINUE
              RNAME(I0)=TANH(RNAME(I1))
              GO TO 199
  199       CONTINUE
          END IF
C
          IF(I0.LE.2*MFILE) THEN
            IF(I0.LE.MFILE) THEN
              GRID(KGRID0(I0)+IGRID)=RNAME(I0)
            ELSE
              GRID(KGRID0(I0-MFILE)+IGRID)=RNAME(I0)
            END IF
          END IF
  201   CONTINUE
C
  202 CONTINUE
C
C.......................................................................
C
C     Writing output grid values:
      DO 339 JFILE=1,2*MFILE
        IFILE=MOD(JFILE-1,MFILE)+1
        DO 332 IKOM=1,NKOM
          IF(KOM0(IKOM).EQ.JFILE) THEN
C           File appears at the L.H.S. of the command:
            IF(KOM3(IKOM).EQ.7.OR.KOM3(IKOM).EQ.8) THEN
C             Integer values (Results of function INT or NINT):
              DO 331 I=KGRID0(IFILE)+1,KGRID0(IFILE)+NGRID
                IRAM(I)=NINT(GRID(I))
  331         CONTINUE
              IF(LARRAY) THEN
                IF(ITIME.EQ.1.OR.JFILE.GT.MFILE) THEN
                  IF(ITIME.EQ.1) THEN
                    OPEN(LU(IFILE),FILE=FILE(IFILE),FORM='FORMATTED')
                  END IF
                  CALL WARRAI(LU(IFILE),' ','FORMATTED',.TRUE.,-999998 ,
     *                     .FALSE.,0 ,N1*N2*N3,   IRAM(KGRID0(IFILE)+1))
                  IF(ITIME.EQ.NTIME.OR.JFILE.LE.MFILE) THEN
                    CLOSE(LU(IFILE))
                  END IF
                END IF
              ELSE
                CALL WARAI(LU1,FILE(IFILE),'FORMATTED',.TRUE.,-999998 ,
     *                     .FALSE.,0 ,N1*N2*N3,N4,IRAM(KGRID0(IFILE)+1))
              END IF
            ELSE
C             Output values may be real:
              IF(LARRAY) THEN
                IF(ITIME.EQ.1.OR.JFILE.GT.MFILE) THEN
                  IF(ITIME.EQ.1) THEN
                    OPEN(LU(IFILE),FILE=FILE(IFILE),FORM='FORMATTED')
                  END IF
                  CALL WARRAY(LU(IFILE),' ','FORMATTED',.TRUE.,-999998.,
     *                     .FALSE.,0.,N1*N2*N3,   GRID(KGRID0(IFILE)+1))
                  IF(ITIME.EQ.NTIME.OR.JFILE.LE.MFILE) THEN
                    CLOSE(LU(IFILE))
                  END IF
                END IF
              ELSE
                CALL WARAY(LU1,FILE(IFILE),'FORMATTED',.TRUE.,-999998.,
     *                     .FALSE.,0.,N1*N2*N3,N4,GRID(KGRID0(IFILE)+1))
              END IF
            END IF
            GO TO 333
          END IF
  332   CONTINUE
  333   CONTINUE
  339 CONTINUE
C
  900 CONTINUE
      WRITE(*,'(A)')
     *    '+GRDCAL: Done.                                              '
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE REGNAM(NAME0,NAME,MNAME,NNAME,KOM)
C
      INTEGER MNAME,NNAME,KOM
      CHARACTER*(*) NAME0,NAME(MNAME)
C
C-----------------------------------------------------------------------
C
      INTEGER INAME
C
      DO 10 INAME=1,NNAME
        IF(NAME(INAME).EQ.NAME0) THEN
          KOM=INAME
          GO TO 20
        END IF
   10 CONTINUE
      NNAME=NNAME+1
      IF(NNAME.GT.MNAME) THEN
C       GRDCAL-38
        CALL ERROR('GRDCAL-38: Insufficient memory for variable names')
C       Maximum number MNAME of variables used in the command file
C       should probably be increased.  MNAME is declared by the
C       PARAMETER statement.
      END IF
      NAME(NNAME)=NAME0
      KOM=NNAME
C
   20 CONTINUE
      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