C
C Program GRDCAL (GRiD CALculator) to perform vectorial calculations
C with real-valued arrays stored in disk files.
C
C Version: 6.50
C Date: 2011, May 23
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
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     CAL='string'...Name of the input file containing the commands
C             to be carried out at each gridpoint.
C             Should always be specified.
C             Description of file CAL
C             Default: CAL='.cal'
C     GRD1='string', GRD2='string', ... ,GRD9='string'... Names
C             of the input/output ASCII files with the grid values.
C             Whether a file is input, output, input/output or not
C             used is dependent on the pseudo-code in file CAL.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Default: GRD1=' ', ... , GRD9=' '
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 Modification to enable operations with multivalued grids:
C     GRDCALNUM='string'... Optional name of the file containing
C             N1*N2*N3 integers specifying the numbers of values
C             at the gridpoints of the multivalued grid.
C             If GRDCALNUM is specified, files GRD1, GRD2, ..., GRD9
C             are assumed to have form MGRD (Multivalued GRiDded data).
C             Note that parameter GRDCALNUM is equivalent to parameter
C             NUM of some other programs, and is called here differently
C             in order to avoid unintended application of parameter NUM
C             to this program.
C             Default: GRDCALNUM=' ' (no modification)
C Modification to enable operations with matrix elements:
C     If N1=0 or N2=0, parameters M1 or M2 are used, respectively.
C     M1='string'... Name of the file containing a single integer number
C             specifying N1.
C             Default: M1=' ' means that N1=1.
C     M2='string'... Name of the file containing a single integer number
C             specifying N2.
C             Option for symmetric matrices:
C             If the string value of M2 equals the string value of M1,
C             then N1=m1*(m1+1)/2 and N2=1, where m1 is read from M1.
C             Default: M2=' ' means that N2=1.
C     If N1=0 or N2=0, you probably wish to specify N3=1 and N4=1.
C Form of the output files with matrices:
C     FORMM='string' ... Form of the output files with matrices. Allowed
C             values are FORMM='formatted' and FORMM='unformatted'.
C Optional parameters specifying the form of the quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C     NUMLIN=positive integer ... Number of reals or integers in one
C             line of the output file. See the description in file
C             forms.for.
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
      EXTERNAL LENGTH,UARRAY
      INTEGER LENGTH
      REAL UARRAY
      INTEGER MFILE,MNAME,MKOM,LU1
      INTEGER JFILE,IFILE,KEQ,KEND,L,IGRID,NGRID,IKOM,ITIME,NTIME
C     
      PARAMETER (MFILE=9,MNAME=2*MFILE+20,MKOM=100,LU1=1)
      INTEGER IUNDEF
      REAL UNDEF
      PARAMETER (IUNDEF=-999999999)
      CHARACTER*80 FGRD,FKOM,FILE(MFILE),FM1,FM2
      CHARACTER*4 GRDN
      INTEGER KGRID0(MFILE),KGRID1(MFILE),LU(MFILE)
      CHARACTER*80 NAME(MNAME)
      REAL RNAME(MNAME)
      INTEGER KOM0(MKOM),KOM1(MKOM),KOM2(MKOM),KOM3(MKOM)
      CHARACTER*13 FORMM
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,MATRIX,I,K
      INTEGER N1,N2,N3,N4,M1,I0,I1,I2
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
      UNDEF=UARRAY()
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GRDCAL: Enter input filename: '
      FGRD=' '
      READ(*,*) FGRD
      WRITE(*,'(A)') '+GRDCAL: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FGRD.NE.' ') THEN
        CALL RSEP1(LU1,FGRD)
      ELSE
C       GRDCAL-45
        CALL ERROR('GRDCAL-45: 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 name of the file CAL:
      CALL RSEP3T('CAL',FKOM,'.cal')
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 names of the files GRDn:
      GRDN='GRD0'
      DO 9 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)
        GRDN(4:4)=CHAR(ICHAR('0')+IFILE)
        CALL RSEP3T(GRDN,FILE(IFILE),' ')
    9 CONTINUE
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
C     Modification to enable operations with multivalued grids:
      CALL RSEP3T('GRDCALNUM',FM1,' ')
      IF(FM1.NE.' ') THEN
        IF(N1*N2*N3.GT.MRAM) THEN
C         GRDCAL-46
          CALL ERROR('GRDCAL-46: Insufficient memory for GRDCALNUM')
C         Number N1*N2*N3 of integers specifying the numbers of values
C         at the gridpoints of the multivalued grid exceeds dimension
C         MRAM of array RAM in include file ram.inc.
        END IF
        CALL RARRAI(LU1,FM1,'FORMATTED',.TRUE.,0,N1*N2*N3,IRAM)
        I1=0
        DO 10 I=1,N1*N2*N3
          I1=I1+IRAM(I)
   10   CONTINUE
        N1=I1
        N2=1
        N3=1
      END IF
C
C     Modification to enable operations with matrix elements
      MATRIX=-1
      FM1=' '
      IF(N1.LE.0) THEN
        MATRIX=1
        CALL RSEP3T('M1',FM1,' ')
        IF(FM1.EQ.' ') THEN
          N1=1
        ELSE
          OPEN(LU1,FILE=FM1,STATUS='OLD')
          READ(LU1,*) N1
          CLOSE(LU1)
        ENDIF
      END IF
      M1=N1
      IF(N2.LE.0) THEN
        MATRIX=1
        CALL RSEP3T('M2',FM2,' ')
        IF(FM2.EQ.' ') THEN
          N2=1
        ELSE IF(FM2.EQ.FM1) THEN
          MATRIX=0
          M1=N1
          N1=N1*(N1+1)/2
          N2=1
        ELSE
          OPEN(LU1,FILE=FM2,STATUS='OLD')
          READ(LU1,*) N2
          CLOSE(LU1)
          IF(N1.EQ.1) THEN
            M1=N2
            N1=N2
            N2=1
          ENDIF
        ENDIF
      END IF
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),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
                  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
                IF(K+1.GT.KEND) THEN
C                 
C                 GRDCAL-04
                  WRITE(*,'(2A)') '                 ',LINE(KEQ:KEND)
                  CALL ERROR('GRDCAL-04: Missing second operand')
                END IF
                CALL REGNAM(LINE(K+1:KEND),NAME,MNAME,NNAME,KOM2(NKOM))
              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=2*NFILE+1,NNAME
      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
C-530   IF(JFILE.EQ.MFILE) THEN
C         The part of RAM reserved for input spatial grids ($-files):
C-530     JGRID=IGRID
C-530   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
C       be increased 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:
*530  IF(N4.LE.1) THEN
*530    IGRID=0
*530  ELSE
C       Protecting the part of RAM with input spatial grids ($-files):
*530    IGRID=JGRID
*530  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
            IF(KGRID1(IFILE).GE.0) THEN
              KGRID0(IFILE)=KGRID1(IFILE)
            ELSE
              KGRID0(IFILE)=IGRID
              IGRID=IGRID+N1*N2*N3
            END IF
            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
C       be increased 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 53 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:
              IF (MATRIX.LT.0) THEN
                CALL RARRAY(LU(IFILE),FILE(IFILE),'FORMATTED',.TRUE.,
     *                          UNDEF,N1*N2*N3,   GRID(KGRID1(IFILE)+1))
              ELSE
                CALL RMAT(LU(IFILE),FILE(IFILE),
     *                  M1,MATRIX*N2*N3,GRID(KGRID1(IFILE)+1))
              ENDIF
              GO TO 52
            END IF
   51     CONTINUE
   52     CONTINUE
   53   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 54 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
   54         CONTINUE
              IF(ITIME.EQ.1) THEN
                IF (MATRIX.LT.0) THEN
                  OPEN(LU(IFILE),FILE=FILE(IFILE),FORM='FORMATTED')
                ELSE
                  CALL OMAT(LU(IFILE),FILE(IFILE),1,FORMM)
                END IF
              END IF
              IF (MATRIX.LT.0) THEN
                CALL RARRAY(LU(IFILE),' ' ,'FORMATTED',.TRUE.,UNDEF,
     *                      N1*N2*N3,  GRID(KGRID1(IFILE)+1))
              ELSE
                CALL RMAT(LU(IFILE),' ',
     *                    M1,MATRIX*N2*N3,GRID(KGRID1(IFILE)+1))
              ENDIF
              IF(ITIME.EQ.NTIME) THEN
                CLOSE(LU(IFILE))
              END IF
            ELSE
              IF (MATRIX.LT.0) THEN
                CALL RARAY(LU1,FILE(IFILE),'FORMATTED',.TRUE.,UNDEF,
     *                     N1*N2*N3,N4,GRID(KGRID1(IFILE)+1))
              ELSE
                CALL OMAT(LU1,FILE(IFILE),1,FORMM)
                DO 55 I=0,N4-1
                  CALL RMAT(LU1,' ',
     *                 M1,MATRIX*N2*N3,GRID(KGRID1(IFILE)+N1*N2*N3*I+1))
   55           CONTINUE
                CLOSE(LU1)
              ENDIF
            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).EQ.UNDEF) 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).EQ.UNDEF) THEN
              LUNDEF=.TRUE.
            END IF
          END IF
          IF(LUNDEF) THEN
            RNAME(I0)=UNDEF
          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)=UNDEF
                END IF
              ELSE
                RNAME(I0)=RNAME(I1)/RNAME(I2)
              END IF
              GO TO 199
  105       CONTINUE
              IF(RNAME(I1).LT.0.) THEN
                RNAME(I0)=UNDEF
              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)=UNDEF
              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)=UNDEF
              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)=UNDEF
              ELSE
                RNAME(I0)=ALOG(RNAME(I1))
              END IF
              GO TO 199
  117       CONTINUE
              IF(RNAME(I1).LE.0.) THEN
                RNAME(I0)=UNDEF
              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)=UNDEF
              ELSE
                RNAME(I0)=ASIN(RNAME(I1))
              END IF
              GO TO 199
  122       CONTINUE
              IF(ABS(RNAME(I1)).GT.1.) THEN
                RNAME(I0)=UNDEF
              ELSE
                RNAME(I0)=ACOS(RNAME(I1))
              END IF
              GO TO 199
  123       CONTINUE
              IF(ABS(RNAME(I1)).GT.1.) THEN
                RNAME(I0)=UNDEF
              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)=UNDEF
              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 337 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
                IF(GRID(I).EQ.UNDEF) THEN
                  IRAM(I)=IUNDEF
                ELSE
                  IRAM(I)=NINT(GRID(I))
                END IF
  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.,IUNDEF,
     *                     .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.,IUNDEF,
     *                     .FALSE.,0 ,N1*N2*N3,N4,IRAM(KGRID0(IFILE)+1))
              END IF
            ELSE
C             Output values may be real-valued:
              IF(LARRAY) THEN
                IF(ITIME.EQ.1.OR.JFILE.GT.MFILE) THEN
                  IF(ITIME.EQ.1) THEN
                    IF(MATRIX.LT.0) THEN
                      OPEN(LU(IFILE),FILE=FILE(IFILE),FORM='FORMATTED')
                    ELSE
                      CALL OMAT(LU(IFILE),FILE(IFILE),2,FORMM)
                    END IF
                  END IF
                  IF(MATRIX.LT.0) THEN
                    CALL WARRAY(LU(IFILE),' ','FORMATTED',.TRUE.,UNDEF,
     *                     .FALSE.,0.,N1*N2*N3,   GRID(KGRID0(IFILE)+1))
                  ELSE
                    DO 333 I=KGRID0(IFILE)+1,KGRID0(IFILE)+N1*N2*N3
                      IF(GRID(I).EQ.UNDEF) THEN
C                       
C                       GRDCAL-43
                       CALL ERROR('GRDCAL-43: Undefined matrix element')
C                       Undefined values are not allowed during matrix
C                       operations, unlike for the grid operations.
                      END IF
  333               CONTINUE
                    CALL WMAT(LU(IFILE),' ',
     *                            M1,MATRIX*N2*N3,GRID(KGRID0(IFILE)+1))
                  END IF
                  IF(ITIME.EQ.NTIME.OR.JFILE.LE.MFILE) THEN
                    CLOSE(LU(IFILE))
                  END IF
                END IF
              ELSE
                IF(MATRIX.LT.0) THEN
                  CALL WARAY(LU1,FILE(IFILE),'FORMATTED',.TRUE.,UNDEF,
     *                     .FALSE.,0.,N1*N2*N3,N4,GRID(KGRID0(IFILE)+1))
                ELSE
                  DO 334 I=KGRID0(IFILE)+1,KGRID0(IFILE)+N1*N2*N3*N4
                    IF(GRID(I).EQ.UNDEF) THEN
C                     
C                     GRDCAL-44
                      CALL ERROR('GRDCAL-44: Undefined matrix element')
C                     Undefined values are not allowed during matrix
C                     operations, unlike for the grid operations.
                    END IF
  334             CONTINUE
                  CALL OMAT(LU1,FILE(IFILE),2,FORMM)
                  DO 335 I=0,N4-1
                    CALL WMAT(LU1,' ',
     *                 M1,MATRIX*N2*N3,GRID(KGRID0(IFILE)+N1*N2*N3*I+1))
  335             CONTINUE
                  CLOSE(LU1)
                END IF
              END IF
            END IF
            GO TO 338
          END IF
  337   CONTINUE
  338   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
      INCLUDE 'mat.for'
C     mat.for
      INCLUDE 'indexi.for'
C     indexi.for.
C
C=======================================================================
C