C
C Program GRDCAL (GRiD CALculator) to perform vectorial calculations C with real-valued arrays stored in disk files. C C Version: 8.00 C Date: 2022, April 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 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 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 Whether a file contains real or integer grid values is C also 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 input/output files: C FORM='string'... Form of the input/output files: either C 'FORMATTED' or 'UNFORMATTED' (case insensitive). C Default: FORM='FORMATTED' C FORMW='string'... Form of the output file. If both FORM and FORMW C are specified, FORMW is used for output files. C Default: FORMW=FORM C FORMR='string'... Form of the input file. If both FORM and FORMR C are specified, FORMR is used for output files. C Default: FORMR=FORM 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 Value of undefined quantities: C UNDEF=real... The value to be used for undefined real quantities. C Default: UNDEF=undefined value used in 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 Above operations and functions are performed with reals. C However, the commands of forms like C $2=INT(A) C $3=NINT(A) C imply integer output values, which are written as integers to the C output files. C Analogously, the commands of forms like C A=REAL($1) C B=FLOAT($2) C imply integer input values, which are read as integers from the C input 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 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. END IF C UNDEF=UARRAY() 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,.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) END IF 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 END IF END IF 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 IF(LINE(KEQ+1:K-1).EQ.'real'.OR. * LINE(KEQ+1:K-1).EQ.'float') THEN KOM3(NKOM)=28 IF(KOM2(NKOM).NE.0) THEN C GRDCAL-33 CALL ERROR * ('GRDCAL-33: Redundant argument in REAL or FLOAT') END IF ELSE C GRDCAL-47 WRITE(*,'(2A)') ' ',LINE(KEQ:KEND) CALL ERROR ('GRDCAL-47: 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 58 IFILE=1,MFILE DO 56 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 C Input grid IF(KOM3(IKOM).EQ.28) THEN C Integer input values CALL RARRAI(LU(IFILE),FILE(IFILE),.TRUE., * IUNDEF,N1*N2*N3,IRAM(KGRID1(IFILE)+1)) DO 52 I=KGRID1(IFILE)+1,KGRID1(IFILE)+N1*N2*N3 IF(IRAM(I).EQ.IUNDEF) THEN GRID(I)=UNDEF ELSE GRID(I)=FLOAT(IRAM(I)) END IF 52 CONTINUE ELSE C Real input values CALL RARRAY(LU(IFILE),FILE(IFILE),.TRUE., * UNDEF,N1*N2*N3,GRID(KGRID1(IFILE)+1)) END IF ELSE C Input matrix CALL RMAT(LU(IFILE),FILE(IFILE), * M1,MATRIX*N2*N3,GRID(KGRID1(IFILE)+1)) END IF GO TO 57 END IF 56 CONTINUE 57 CONTINUE 58 CONTINUE END IF C 4-D space-time grids (@) DO 68 JFILE=MFILE+1,2*MFILE IFILE=MOD(JFILE-1,MFILE)+1 DO 66 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 61 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 61 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 C Input grid IF(KOM3(IKOM).EQ.28) THEN C Integer input values CALL RARRAI(LU(IFILE),' ',.TRUE.,IUNDEF, * N1*N2*N3,IRAM(KGRID1(IFILE)+1)) DO 62 I=KGRID1(IFILE)+1,KGRID1(IFILE)+N1*N2*N3 IF(IRAM(I).EQ.IUNDEF) THEN GRID(I)=UNDEF ELSE GRID(I)=FLOAT(IRAM(I)) END IF 62 CONTINUE ELSE C Real input values CALL RARRAY(LU(IFILE),' ',.TRUE.,UNDEF, * N1*N2*N3,GRID(KGRID1(IFILE)+1)) END IF ELSE C Input matrix CALL RMAT(LU(IFILE),' ', * M1,MATRIX*N2*N3,GRID(KGRID1(IFILE)+1)) END IF IF(ITIME.EQ.NTIME) THEN CLOSE(LU(IFILE)) END IF ELSE IF (MATRIX.LT.0) THEN C Input grid IF(KOM3(IKOM).EQ.28) THEN C Integer input values CALL RARAI(LU1,FILE(IFILE),.TRUE.,IUNDEF, * N1*N2*N3,N4,IRAM(KGRID1(IFILE)+1)) DO 63 I=KGRID1(IFILE)+1,KGRID1(IFILE)+N1*N2*N3*N4 IF(IRAM(I).EQ.IUNDEF) THEN GRID(I)=UNDEF ELSE GRID(I)=FLOAT(IRAM(I)) END IF 63 CONTINUE ELSE C Real input values CALL RARAY(LU1,FILE(IFILE),.TRUE.,UNDEF, * N1*N2*N3,N4,GRID(KGRID1(IFILE)+1)) END IF ELSE C Input matrix CALL OMAT(LU1,FILE(IFILE),1,FORMM) DO 65 I=0,N4-1 CALL RMAT(LU1,' ', * M1,MATRIX*N2*N3,GRID(KGRID1(IFILE)+N1*N2*N3*I+1)) 65 CONTINUE CLOSE(LU1) END IF END IF GO TO 67 END IF 66 CONTINUE 67 CONTINUE 68 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),' ',.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),.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),' ',.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),FILE(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),.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,FILE(IFILE), * 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