C 
C Program GRDNEW to trilinearly interpolate grid values into a new grid
C of different dimensions or density
C
C Version: 5.50
C Date: 2001, January 9
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 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 real quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
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 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
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)
            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=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