C
C Program GRDFD to calculate gradient of the grid values by means of the C the second-order finite differences C C Version: 5.30 C Date: 1999, March 26 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','GRD','GRD1','GRD2','GRD3',/ C 'SEP'...String in apostrophes containing the name of the input C file with the data specifying dimensions of the input C grid. C 'GRD'...String in apostrophes containing the name of the input C ASCII files with the grid values. C 'GRD1','GRD2','GRD3',... Strings in apostrophes containing the C names of the output ASCII files containing the first C partial derivatives in the respective directions. C Derivatives corresponding to a blank filename are not C evaluated. C /... Input data line must be terminated by a slash. C Default: 'SEP'='grd.h', 'GRD'='grd.out', 'GRD1'=' ', C 'GRD2'=' ', 'GRD3'=' '. 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. C Default: N4=1 C D1=real... Grid interval in the direction of the first coordinate C axis. C Default: D1=1. C D2=real... Grid interval in the direction of the second coordinate C axis. C Default: D2=1. C D3=real... Grid interval in the direction of the third coordinate C axis. C Default: D3=1. C Data specific to this program: C NORDER=even positive integer... Order of finite difference scheme. C Only NORDER=2 is now coded. C Default: NORDER=2 C NHALF=integer... Kind of finite difference scheme: C NHALF=0: Ordinary centred 1-D scheme. Output grids have C the same dimensions as the input grid. C Differences at boundaries and in the vicinity of C undefined values are approximated by single-sided C differences. C Error of derivative F1 for NORDER=2: C Err(F1)=4*F111*D1**2/24 C NHALF=1: Half-interval 1-D scheme. Output grids have C dimensions: 'GRD1': (N1-1)*N2*N3 C 'GRD2': N1*(N2-1)*N3 C 'GRD3': N1*N2*(N3-1) C Error of derivative F1 for NORDER=2: C Err(F1)= F111*D1**2/24 C NHALF=2: Averaged half-interval 1-D schemes. Output grids C have dimensions (N1-1)*(N2-1)*(N3-1) gridpoints C if N1, N2 and N3 are greater than 1. C Naturally, any dimension Ni=1 is not changed. C Error of derivative F1 for NORDER=2: C Err(F1)=(F111*D1**2+F122*D2**2+F122*D2**2)/24 C Default: NHALF=0 C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C....................................................................... C C Filenames and parameters: CHARACTER*80 FILE1,FILE2,FILE3,FILE4,FILE5 INTEGER LU1,LU2,LU3,LU4 REAL UNDEF PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,UNDEF=-999999.) C Input data: INTEGER N1,N2,N3,N4 REAL D1,D2,D3 C Other variables: INTEGER N12,N123,N1234,I1,I2,I3,I4,I REAL D,DD,F1,F2,F3 C C....................................................................... C C Reading main input data: FILE1='grd.h' FILE2='grd.out' FILE3=' ' FILE4=' ' FILE5=' ' WRITE(*,'(2A)') '+GRDFD: Enter 5 filenames ', * '[''grd.h'',''grd.out'','' '','' '','' '']: ' READ(*,*) FILE1,FILE2,FILE3,FILE4,FILE5 WRITE(*,'(A)') * '+GRDFD: Reading, calculating, writing... ' C C Reading grid dimensions: C Original grid: CALL RSEP1(LU1,FILE1) CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3I('N4',N4,1) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) N12 =N1*N2 N123 =N1*N2*N3 C4 N1234=N1*N2*N3*N4 N1234=N1*N2*N3 IF(2*N1234.GT.MRAM) THEN C GRDFD-01 CALL ERROR('GRDFD-01: Too small array RAM(MRAM)') C Too small array RAM(MRAM) to allocate both input and output C grid values. If possible, increase dimension MRAM in include C file ram.inc. END IF C C Reading parameters of FD scheme: CALL RSEP3I('NORDER',NORDER,2) IF(NORDER.EQ.2) THEN CONTINUE ELSE C GRDFD-02 CALL ERROR('GRDFD-02: Incorrect value of NORDER') C Allowed values of input parameter NORDER: 2. C See the description of input data SEP. END IF CALL RSEP3I('NHALF',NHALF,0) IF(NHALF.EQ.0) THEN CONTINUE ccc ELSE IF(NHALF.EQ.1) THEN ccc CONTINUE ccc ELSE IF(NHALF.EQ.2) THEN ccc CONTINUE ELSE C GRDFD-03 CALL ERROR('GRDFD-03: Incorrect value of NHALF') C Allowed values of input parameter NHALF: 0. C See the description of input data SEP. END IF C C Reading input grid values: OPEN(LU1,FILE=FILE2,FORM='FORMATTED',STATUS='OLD') IF(FILE3.NE.' ') THEN OPEN(LU2,FILE=FILE3,FORM='FORMATTED') END IF IF(FILE4.NE.' ') THEN OPEN(LU3,FILE=FILE4,FORM='FORMATTED') END IF IF(FILE5.NE.' ') THEN OPEN(LU4,FILE=FILE5,FORM='FORMATTED') END IF C C Loop over time levels: DO 34 I4=0,N4-1 CALL RARRAY(LU1,' ','FORMATTED',.TRUE.,UNDEF,N1234,RAM) C C Calculating X1 derivatives: IF(FILE3.NE.' ') THEN IF(N1.LE.1) THEN C GRDFD-04 CALL ERROR('GRDFD-04: N1 is less than 2') END IF D=D1 DD=2.*D I=0 J=N1234 DO 13 I3=1,N3 DO 12 I2=1,N2 I=I+1 J=J+1 F2=RAM(I) F3=RAM(I+1) IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF DO 11 I1=2,N1-1 I=I+1 J=J+1 F1=F2 F2=F3 F3=RAM(I+1) IF(F1.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F1)/DD ELSE IF(F1.NE.UNDEF.AND.F2.NE.UNDEF) THEN RAM(J)=(F2-F1)/D ELSE IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF 11 CONTINUE I=I+1 J=J+1 IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF 12 CONTINUE 13 CONTINUE C Writing output grid values: CALL WARRAY(LU2,' ','FORMATTED',.TRUE.,UNDEF,.FALSE.,0., * N1234,RAM(N1234+1)) END IF C C Calculating X2 derivatives: IF(FILE4.NE.' ') THEN IF(N2.LE.1) THEN C GRDFD-05 CALL ERROR('GRDFD-05: N2 is less than 2') END IF D=D2 DD=2.*D DO 23 I3=0,N3-1 DO 22 I1=1,N1 I=I3*N12+I1 J=N1234+I F2=RAM(I) F3=RAM(I+N1) IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF DO 21 I2=1,N2-2 I=I+N1 J=J+N1 F1=F2 F2=F3 F3=RAM(I+N1) IF(F1.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F1)/DD ELSE IF(F1.NE.UNDEF.AND.F2.NE.UNDEF) THEN RAM(J)=(F2-F1)/D ELSE IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF 21 CONTINUE I=I+N1 J=J+N1 IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF 22 CONTINUE 23 CONTINUE C Writing output grid values: CALL WARRAY(LU3,' ','FORMATTED',.TRUE.,UNDEF,.FALSE.,0., * N1234,RAM(N1234+1)) END IF C C Calculating X3 derivatives: IF(FILE5.NE.' ') THEN IF(N3.LE.1) THEN C GRDFD-06 CALL ERROR('GRDFD-06: N3 is less than 2') END IF D=D3 DD=2.*D DO 33 I2=0,N2-1 DO 32 I1=1,N1 I=I2*N1+I1 J=N1234+I F2=RAM(I) F3=RAM(I+N12) IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF DO 31 I3=1,N3-2 I=I+N12 J=J+N12 F1=F2 F2=F3 F3=RAM(I+N12) IF(F1.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F1)/DD ELSE IF(F1.NE.UNDEF.AND.F2.NE.UNDEF) THEN RAM(J)=(F2-F1)/D ELSE IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF 31 CONTINUE I=I+N12 J=J+N12 IF(F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN RAM(J)=(F3-F2)/D ELSE RAM(J)=UNDEF END IF 32 CONTINUE 33 CONTINUE C Writing output grid values: CALL WARRAY(LU4,' ','FORMATTED',.TRUE.,UNDEF,.FALSE.,0., * N1234,RAM(N1234+1)) END IF C 34 CONTINUE CLOSE(LU1) IF(FILE3.NE.' ') THEN CLOSE(LU2) END IF IF(FILE4.NE.' ') THEN CLOSE(LU3) END IF IF(FILE5.NE.' ') THEN CLOSE(LU4) END IF WRITE(*,'(A)') * '+GRDFD: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= C