C
C Program GRDNEW to trilinearly interpolate grid values into a new grid C of different dimensions or density C C Version: 5.30 C Date: 1999, June 11 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','GRDNEW',/ C 'SEP'...String in apostrophes containing the name of the input C SEP parameter file with the data specifying the dimensions C of the input and output grids. C Description of file SEP C 'GRD'...String in apostrophes containing the name of the input C ASCII files with the grid values. C 'GRDNEW'... String in apostrophes containing the name of the C output ASCII file containing the grid values interpolated C into a new grid. C /... Input data line must be terminated by a slash. C Default: 'SEP'='grd.h', 'GRD'='grd.out', 'GRDNEW'='grdnew.out'. 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 dimensions of the input grid: 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 O1=real... First coordinate of the grid origin (first point of the C grid). C Default: O1=0. C O2=real... Second coordinate of the grid origin. C Default: O2=0. C O3=real... Third coordinate of the grid origin. C Default: O3=0. 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 specifying dimensions of the output grid: C N1NEW=positive integer C N2NEW=positive integer C N3NEW=positive integer C O1NEW=real C O2NEW=real C O3NEW=real C D1NEW=real C D2NEW=real C D3NEW=real... Analogous to N1, N2, N3, O1, O2, O3, D1, D2 and D3, C respectively, but for the output grid. C Defaults: N1NEW=N1, N2NEW=N2, N3NEW=N3, C O1NEW=O1, O2NEW=O2, O3NEW=O3, C D1NEW=D1, D2NEW=D2, D3NEW=D3. 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 INTEGER LU REAL UNDEF PARAMETER (LU=1,UNDEF=-999999.) C Input data: INTEGER N1,N2,N3,N1NEW,N2NEW,N3NEW REAL O1,O2,O3,D1,D2,D3,O1NEW,O2NEW,O3NEW,D1NEW,D2NEW,D3NEW C Other variables: INTEGER N1N2,N1N2N3,I1,I2,I3 INTEGER INEW,I1NEW,I2NEW,I3NEW,I1MIN,I2MIN,I3MIN,I1MAX,I2MAX,I3MAX INTEGER I000,I100,I010,I110,I001,I101,I011,I111 REAL A000,A100,A010,A110,A001,A101,A011,A111,A00,A10,A01,A11,A0,A1 REAL B0,B1,X1,X2,X3 C C....................................................................... C C Reading main input data: FILE1='grd.h' FILE2='grd.out' FILE3='grdnew.out' WRITE(*,'(A)') * '+GRDNEW: Enter 2 filenames [grd.h, grd.out, grdnew.out]: ' READ(*,*) FILE1,FILE2,FILE3 WRITE(*,'(A)') * '+GRDNEW: Working... ' C C Reading grid dimensions: C Original grid: CALL RSEP1(LU,FILE1) CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) C New grid: CALL RSEP3I('N1NEW',N1NEW,N1) CALL RSEP3I('N2NEW',N2NEW,N2) CALL RSEP3I('N3NEW',N3NEW,N3) CALL RSEP3R('O1NEW',O1NEW,O1) CALL RSEP3R('O2NEW',O2NEW,O2) CALL RSEP3R('O3NEW',O3NEW,O3) CALL RSEP3R('D1NEW',D1NEW,D1) CALL RSEP3R('D2NEW',D2NEW,D2) CALL RSEP3R('D3NEW',D3NEW,D3) C C Dimensions of the old grid related to the new one: O1=(O1-O1NEW)/D1NEW O2=(O2-O2NEW)/D2NEW O3=(O3-O3NEW)/D3NEW D1=D1/D1NEW D2=D2/D2NEW D3=D3/D3NEW N1N2 =N1*N2 N1N2N3=N1*N2*N3 C IF(N1N2N3+N1NEW*N2NEW*N3NEW.GT.MRAM) THEN C GRDNEW-01 CALL ERROR('GRDNEW-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 input grid values: CALL RARRAY(LU,FILE2,'FORMATTED',.TRUE.,UNDEF,N1N2N3,RAM) C C Trilinearly interpolating inside the grid: I3MIN=0 DO 23 I3=0,N3 IF(I3.LT.N3) THEN X3=O3+FLOAT(I3)*D3 I3MAX=MIN0(INT(X3),N3NEW-1) ELSE I3MAX=N3NEW-1 END IF I2MIN=0 DO 22 I2=0,N2 IF(I2.LT.N2) THEN X2=O2+FLOAT(I2)*D2 I2MAX=MIN0(INT(X2),N2NEW-1) ELSE I2MAX=N2NEW-1 END IF I1MIN=0 DO 21 I1=0,N1 IF(I1.LT.N1) THEN X1=O1+FLOAT(I1)*D1 I1MAX=MIN0(INT(X1),N1NEW-1) ELSE I1MAX=N1NEW-1 END IF I111=1+I1+N1*(I2+N2*I3) IF(I3.GT.0) THEN I110=I111-N1N2 IF(I3.GE.N3) THEN I111=I110 END IF ELSE I110=I111 END IF IF(I2.GT.0) THEN I100=I110-N1 I101=I111-N1 IF(I2.GE.N2) THEN I110=I100 I111=I101 END IF ELSE I100=I110 I101=I111 END IF IF(I1.GT.0) THEN I000=I100-1 I010=I110-1 I001=I101-1 I011=I111-1 IF(I1.GE.N1) THEN I100=I000 I110=I010 I101=I001 I111=I011 END IF ELSE I000=I100 I010=I110 I001=I101 I011=I111 END IF A000=RAM(I000) A100=RAM(I100) A010=RAM(I010) A110=RAM(I110) A001=RAM(I001) A101=RAM(I101) A011=RAM(I011) A111=RAM(I111) DO 13 I3NEW=I3MIN,I3MAX B0=(X3-FLOAT(I3NEW))/D3 B1=1.-B0 A00=A000*B0+A001*B1 A10=A100*B0+A101*B1 A01=A010*B0+A011*B1 A11=A110*B0+A111*B1 DO 12 I2NEW=I2MIN,I2MAX B0=(X2-FLOAT(I2NEW))/D2 B1=1.-B0 A0=A00*B0+A01*B1 A1=A10*B0+A11*B1 INEW=N1N2N3+I1MIN+N1NEW*(I2NEW+N2NEW*I3NEW) DO 11 I1NEW=I1MIN,I1MAX B0=(X1-FLOAT(I1NEW))/D1 B1=1.-B0 INEW=INEW+1 RAM(INEW)=A0*B0+A1*B1 11 CONTINUE 12 CONTINUE 13 CONTINUE I1MIN=I1MAX+1 21 CONTINUE I2MIN=I2MAX+1 22 CONTINUE I3MIN=I3MAX+1 23 CONTINUE C C Writing output grid values: CALL WARRAY(LU,FILE3,'FORMATTED',.TRUE.,UNDEF,.FALSE.,0., * N1NEW*N2NEW*N3NEW,RAM(N1N2N3+1)) WRITE(*,'(A)') * '+GRDNEW: Done. ' STOP 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