C
C Program GRDFD to calculate gradient of the grid values by means of the
C second-order finite differences
C
C Version: 7.40
C Date: 2017, May 18
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 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     GRD11='string', GRD22='string', GRD33='string'... The names
C             of the output ASCII files containing the homogeneous
C             second partial derivatives in the respective directions.
C             Derivatives corresponding to a blank filename are not
C             evaluated.
C             Default: GRD11=' ', GRD22=' ', GRD33=' '
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             Only NHALF=0 is now coded.
C             Default: NHALF=0
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,FILE4,FILE5,FILE6,FILE7
      INTEGER LU1,LU2,LU3,LU4,LU5,LU6,LU7
      REAL UNDEF
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=11,LU6=12,LU7=13)
C     Input data:
      INTEGER N1,N2,N3,N4,NORDER,NHALF
      REAL D1,D2,D3
C     Other variables:
      INTEGER N12,N123,N1234,I1,I2,I3,I4,I,J,J1
      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-01
        CALL ERROR('GRDFD-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
C
      UNDEF=UARRAY()
C
C     Reading input parameters from the SEP file:
      CALL RSEP3T('GRD'  ,FILE1,' ')
      CALL RSEP3T('GRD1' ,FILE2,' ')
      CALL RSEP3T('GRD2' ,FILE3,' ')
      CALL RSEP3T('GRD3' ,FILE4,' ')
      CALL RSEP3T('GRD11',FILE5,' ')
      CALL RSEP3T('GRD22',FILE6,' ')
      CALL RSEP3T('GRD33',FILE7,' ')
      IF (FILE1.EQ.' ') THEN
C       GRDFD-11
        CALL ERROR('GRDFD-11: File GRD not specified')
C       Name of the input grid file, specified by SEP parameter GRD
C       is blank or not specified.  It must be specified.
C       There is no default filename.
C       See the description of input data SEP.
      ENDIF
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-02
        CALL ERROR('GRDFD-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
C
C     Reading parameters of FD scheme:
      CALL RSEP3I('NORDER',NORDER,2)
      IF(NORDER.EQ.2) THEN
        CONTINUE
      ELSE
C       GRDFD-03
        CALL ERROR('GRDFD-03: 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-04
        CALL ERROR('GRDFD-04: Incorrect value of NHALF')
C       Allowed values of input parameter NHALF: 0.
C       See the description of input data SEP.
      END IF
C
C     Opening input and output files:
      OPEN(LU1,FILE=FILE1,FORM='FORMATTED',STATUS='OLD')
      IF(FILE2.NE.' ') THEN
        OPEN(LU2,FILE=FILE2,FORM='FORMATTED')
      END IF
      IF(FILE3.NE.' ') THEN
        OPEN(LU3,FILE=FILE3,FORM='FORMATTED')
      END IF
      IF(FILE4.NE.' ') THEN
        OPEN(LU4,FILE=FILE4,FORM='FORMATTED')
      END IF
      IF(FILE5.NE.' ') THEN
        OPEN(LU5,FILE=FILE5,FORM='FORMATTED')
      END IF
      IF(FILE6.NE.' ') THEN
        OPEN(LU6,FILE=FILE6,FORM='FORMATTED')
      END IF
      IF(FILE7.NE.' ') THEN
        OPEN(LU7,FILE=FILE7,FORM='FORMATTED')
      END IF
C
C     Loop over time levels:
      DO 90 I4=0,N4-1
        CALL RARRAY(LU1,' ','FORMATTED',.TRUE.,UNDEF,N1234,RAM)
C
C       Calculating X1 derivatives:
        IF(FILE2.NE.' ') THEN
          IF(N1.LT.2) THEN
C           GRDFD-05
            CALL ERROR('GRDFD-05: 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 X11 derivatives:
        IF(FILE5.NE.' ') THEN
          IF(N1.LT.3) THEN
C           GRDFD-06
            CALL ERROR('GRDFD-06: N1 is less than 3')
          END IF
          DD=D1*D1
          I=0
          J=N1234
          DO 18 I3=1,N3
            DO 17 I2=1,N2
              I=I+1
              J=J+1
              F2=RAM(I)
              F3=RAM(I+1)
              J1=J
              DO 16 I1=2,N1-1
                I=I+1
                J=J+1
                F1=F2
                F2=F3
                F3=RAM(I+1)
                IF(F1.NE.UNDEF.AND.F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN
                  RAM(J)=(F1-F2-F2+F3)/DD
                ELSE
                  RAM(J)=UNDEF
                END IF
   16         CONTINUE
              I=I+1
              J=J+1
              RAM(J1)=RAM(J1+1)
              RAM(J)=RAM(J-1)
   17       CONTINUE
   18     CONTINUE
C         Writing output grid values:
          CALL WARRAY(LU5,' ','FORMATTED',.TRUE.,UNDEF,.FALSE.,0.,
     *                                               N1234,RAM(N1234+1))
        END IF
C
C       Calculating X2 derivatives:
        IF(FILE3.NE.' ') THEN
          IF(N2.LT.2) THEN
C           GRDFD-07
            CALL ERROR('GRDFD-07: 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 X22 derivatives:
        IF(FILE6.NE.' ') THEN
          IF(N2.LT.3) THEN
C           GRDFD-08
            CALL ERROR('GRDFD-08: N2 is less than 3')
          END IF
          DD=D2*D2
          DO 28 I3=0,N3-1
            DO 27 I1=1,N1
              I=I3*N12+I1
              J=N1234+I
              F2=RAM(I)
              F3=RAM(I+N1)
              J1=J
              DO 26 I2=1,N2-2
                I=I+N1
                J=J+N1
                F1=F2
                F2=F3
                F3=RAM(I+N1)
                IF(F1.NE.UNDEF.AND.F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN
                  RAM(J)=(F1-F2-F2+F3)/DD
                ELSE
                  RAM(J)=UNDEF
                END IF
   26         CONTINUE
              I=I+N1
              J=J+N1
              RAM(J1)=RAM(J1+N1)
              RAM(J)=RAM(J-N1)
   27       CONTINUE
   28     CONTINUE
C         Writing output grid values:
          CALL WARRAY(LU6,' ','FORMATTED',.TRUE.,UNDEF,.FALSE.,0.,
     *                                               N1234,RAM(N1234+1))
        END IF
C
C       Calculating X3 derivatives:
        IF(FILE4.NE.' ') THEN
          IF(N3.LT.2) THEN
C           GRDFD-09
            CALL ERROR('GRDFD-09: 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
C       Calculating X33 derivatives:
        IF(FILE7.NE.' ') THEN
          IF(N3.LT.3) THEN
C           GRDFD-10
            CALL ERROR('GRDFD-10: N3 is less than 3')
          END IF
          DD=D3*D3
          DO 38 I2=0,N2-1
            DO 37 I1=1,N1
              I=I2*N1+I1
              J=N1234+I
              F2=RAM(I)
              F3=RAM(I+N12)
              J1=J
              DO 36 I3=1,N3-2
                I=I+N12
                J=J+N12
                F1=F2
                F2=F3
                F3=RAM(I+N12)
                IF(F1.NE.UNDEF.AND.F2.NE.UNDEF.AND.F3.NE.UNDEF) THEN
                  RAM(J)=(F1-F2-F2+F3)/DD
                ELSE
                  RAM(J)=UNDEF
                END IF
   36         CONTINUE
              I=I+N12
              J=J+N12
              RAM(J1)=RAM(J1+N12)
              RAM(J)=RAM(J-N12)
   37       CONTINUE
   38     CONTINUE
C         Writing output grid values:
          CALL WARRAY(LU7,' ','FORMATTED',.TRUE.,UNDEF,.FALSE.,0.,
     *                                               N1234,RAM(N1234+1))
        END IF
C
   90 CONTINUE
      CLOSE(LU1)
      IF(FILE2.NE.' ') THEN
        CLOSE(LU2)
      END IF
      IF(FILE3.NE.' ') THEN
        CLOSE(LU3)
      END IF
      IF(FILE4.NE.' ') THEN
        CLOSE(LU4)
      END IF
      IF(FILE2.NE.' ') THEN
        CLOSE(LU5)
      END IF
      IF(FILE3.NE.' ') THEN
        CLOSE(LU6)
      END IF
      IF(FILE4.NE.' ') THEN
        CLOSE(LU7)
      END IF
      WRITE(*,'(A)')
     *  '+GRDFD: 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