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.20
C Date: 1998, November 6
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 diference scheme.
C             Only NORDER=2 is now coded.
C             Default: NORDER=2
C     NHALF=integer... Kind of finite diference 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