C
C Program GRDNEW to trilinearly interpolate grid values into a new grid C of different dimensions or density C C Version: 7.40 C Date: 2017, May 25 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 input and output files: C GRD='string'... Name of the input ASCII file with the grid values. C Default: GRD='grd.out' C GRDNEW='string'... Name of the output ASCII file containing the C grid values interpolated into a new grid. C Default: GRDNEW='grdnew.out' C For general description of the files with gridded data refer C to file forms.htm. 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 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 C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C....................................................................... C EXTERNAL UARRAY REAL UARRAY C C Filenames and parameters: CHARACTER*80 FILE1,FILE2,FILE3 INTEGER LU REAL UNDEF PARAMETER (LU=1) 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 LOGICAL LUNDEF C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+GRDNEW: Enter input filename: ' FILE1=' ' READ(*,*) FILE1 WRITE(*,'(A)') '+GRDNEW: Working ... ' C C Reading all data from the SEP file into the memory: IF (FILE1.EQ.' ') THEN C GRDNEW-01 CALL ERROR('GRDNEW-01: 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 CALL RSEP1(LU,FILE1) C UNDEF=UARRAY() C C Reading input parameters from the SEP file: CALL RSEP3T('GRD',FILE2,'grd.out') CALL RSEP3T('GRDNEW',FILE3,'grdnew.out') C C Reading grid dimensions: C Original grid: 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-02 CALL ERROR('GRDNEW-02: 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 IF(D1.LE.0..OR.D2.LE.0..OR.D3.LE.0.) THEN C GRDNEW-03 CALL ERROR('GRDNEW-03: Negative grid interval') C In this version of program 'grdnew.for', grid intervals D1 and C D1NEW must have equal signs, D2 and D2NEW must have equal signs, C D3 and D3NEW must have equal signs. 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(NINT(X3-0.5),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(NINT(X2-0.5),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(NINT(X1-0.5),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) LUNDEF=.FALSE. IF (A000.EQ.UNDEF) LUNDEF=.TRUE. IF (A100.EQ.UNDEF) LUNDEF=.TRUE. IF (A010.EQ.UNDEF) LUNDEF=.TRUE. IF (A110.EQ.UNDEF) LUNDEF=.TRUE. IF (A001.EQ.UNDEF) LUNDEF=.TRUE. IF (A101.EQ.UNDEF) LUNDEF=.TRUE. IF (A011.EQ.UNDEF) LUNDEF=.TRUE. IF (A111.EQ.UNDEF) LUNDEF=.TRUE. 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 IF (LUNDEF) RAM(INEW)=UNDEF 11 CONTINUE 12 CONTINUE 13 CONTINUE I1MIN=MAX0(0,I1MAX+1) 21 CONTINUE I2MIN=MAX0(0,I2MAX+1) 22 CONTINUE I3MIN=MAX0(0,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