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