C 
C Program GRDFD to calculate gradient of the grid values by means of the
C second-order finite differences
C
C Version: 5.40
C Date: 2000, January 22
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 the input and output files:
C     GRD='string'...Name of the input ASCII file with the grid values.
C             Default: GRD='grd.out'
C     GRD1='string', GRD2='string', GRD3='string'... The names
C             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             Default: GRD1=' ', GRD2=' ', GRD3=' '
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 name of SEP file with input data:
      WRITE(*,'(A)') '+GRDFD: Enter input filename: '
      FILE1=' '
      READ(*,*) FILE1
      WRITE(*,'(A)') '+GRDFD: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FILE1.NE.' ') THEN
        CALL RSEP1(LU1,FILE1)
      ELSE
C       GRDFD-07
        CALL ERROR('GRDFD-07: 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
C
C     Reading input parameters from the SEP file:
      CALL RSEP3T('GRD',FILE2,'grd.out')
      CALL RSEP3T('GRD1',FILE3,' ')
      CALL RSEP3T('GRD2',FILE4,' ')
      CALL RSEP3T('GRD3',FILE5,' ')
C
C     Reading grid dimensions:
C     Original grid:
      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