C 
      SUBROUTINE FDIN(LU,NX,NZ,AHL,BHL,CHL,AVL,BVL,CVL,
     *                         AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC,WORK)
C
      INTEGER LU,NX,NZ
      REAL AHL(NZ  ,NX-1),BHL(NZ  ,NX-1),CHL(NZ  ,NX-1)
      REAL AVL(NZ-1,NX  ),BVL(NZ-1,NX  ),CVL(NZ-1,NX  )
      REAL AHS(NZ-1,NX-1),BHS(NZ-1,NX-1),CHS(NZ-1,NX-1)
      REAL AVS(NZ-1,NX-1),BVS(NZ-1,NX-1),CVS(NZ-1,NX-1)
      REAL DEN(NZ  ,NX  ),QFC(NZ  ,NX  ),WORK(NZ ,NX  )
C
C Subroutine FDIN reads the effective elastic parameters for 2-D finite
C difference programs.
C
C Input:
C     LU...   Logical unit number to be used to read all files.
C     NX...   Number of grid points in the horizontal direction.
C     NZ...   Number of grid points in the vertical direction.
C     AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC,WORK...
C             Arrays of the dimensions as declared above or greater.
C             WORK is a temporary working array, its content will be
C             destroyed.
C
C Output:
C     AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC... Values
C             read from formatted or unformatted files given by SEP
C             parameters AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS,
C             CVS,DEN,QFC, respectively.
C     WORK... Undefined.
C
C Date: 1998, November 11
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      INTEGER NEWPAR
C     NEWPAR..NEWPAR=0: Average parameters along the legs are read.
C             NEWPAR=1: Average parameters over the grid facets are read
C               and the parameters corresponding to the legs are derived
C               from them.
C
C.......................................................................
C
      CALL RSEP3I('NEWPAR',NEWPAR,0)
      IF(NEWPAR.EQ.0) THEN
        CALL FDIN1(LU,NX-1,NZ-1,'A13',AHS,WORK,0,AHS)
        CALL FDIN1(LU,NX-1,NZ-1,'B13',BHS,WORK,0,BHS)
        CALL FDIN1(LU,NX-1,NZ-1,'C13',CHS,WORK,0,CHS)
        CALL FDIN1(LU,NX-1,NZ-1,'A31',AVS,WORK,0,AHS)
        CALL FDIN1(LU,NX-1,NZ-1,'B31',BVS,WORK,0,BHS)
        CALL FDIN1(LU,NX-1,NZ-1,'C31',CVS,WORK,0,CHS)
        CALL FDIN1(LU,NX-1,NZ  ,'A11',AHL,WORK,0,AHS)
        CALL FDIN1(LU,NX-1,NZ  ,'B11',BHL,WORK,0,BHS)
        CALL FDIN1(LU,NX-1,NZ  ,'C11',CHL,WORK,0,CHS)
        CALL FDIN1(LU,NX  ,NZ-1,'A33',AVL,WORK,0,AVS)
        CALL FDIN1(LU,NX  ,NZ-1,'B33',BVL,WORK,0,BVS)
        CALL FDIN1(LU,NX  ,NZ-1,'C33',CVL,WORK,0,CVS)
      ELSE
        CALL FDIN1(LU,NX-1,NZ-1,'AA13',AHS,WORK,0,AHS)
        CALL FDIN1(LU,NX-1,NZ-1,'BB13',BHS,WORK,0,BHS)
        CALL FDIN1(LU,NX-1,NZ-1,'CC13',CHS,WORK,0,CHS)
        CALL FDIN1(LU,NX-1,NZ-1,'    ',AVS,WORK,4,AHS)
        CALL FDIN1(LU,NX-1,NZ-1,'    ',BVS,WORK,4,BHS)
        CALL FDIN1(LU,NX-1,NZ-1,'    ',CVS,WORK,4,CHS)
        CALL FDIN1(LU,NX-1,NZ  ,'    ',AHL,WORK,2,AHS)
        CALL FDIN1(LU,NX-1,NZ  ,'    ',BHL,WORK,2,BHS)
        CALL FDIN1(LU,NX-1,NZ  ,'    ',CHL,WORK,2,CHS)
        CALL FDIN1(LU,NX  ,NZ-1,'    ',AVL,WORK,1,AVS)
        CALL FDIN1(LU,NX  ,NZ-1,'    ',BVL,WORK,1,BVS)
        CALL FDIN1(LU,NX  ,NZ-1,'    ',CVL,WORK,1,CVS)
      END IF
      CALL FDIN1(LU,NX  ,NZ  ,'DEN',DEN,WORK,0,DEN)
      CALL FDIN1(LU,NX  ,NZ  ,'QFC',QFC,WORK,0,QFC)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE FDIN1(LU,N1,N2,NAME,ARRAY,WORK,KDEF,DEF)
C
      INTEGER LU,N1,N2,KDEF
      CHARACTER*(*) NAME
      REAL ARRAY(N2*N1),WORK(N1*N2),DEF(*)
C
C Subroutine FDIN1 reads the values of a single effective elastic
C parameter for 2-D finite difference programs.
C
C Input:
C     LU...   Logical unit number to be used for input.
C     N1...   Number of values in the horizontal direction.
C     N2...   Number of values in the vertical direction.
C     NAME... String containing the name of the parameter specifying the
C             name of the file containing the grid values to be read.
C             Except for its case, it should match the parameter name in
C             a previously read input SEP parameter file.  If parameter
C             NAME is not defined in the input SEP parameter file, no
C             input is performed by this subroutine and the output
C             arrays are filled with zeros.
C     KDEF... Kind of default values:
C             KDEF=0: Zeros.  Array DEF is not used.
C             KDEF=1: Array DEF should have dimensions DEF(N2,N1-1).
C                     Then the default is ARRAY(I2,1)=DEF(I2,1),
C                     ARRAY(I2,N1)=DEF(I2,N1-1), and for 1.LT.I1.LT.N1
C                     ARRAY(I2,I1)=(DEF(I2,I1-1)+DEF(I2,I1))/2.
C             KDEF=2: Array DEF should have dimensions DEF(N2-1,N1).
C                     Then the default is ARRAY(1,I1)=0.5*DEF(1,I1),
C                     ARRAY(N2,I1)=DEF(N2-1,I1), and for 1.LT.I2.LT.N2
C                     ARRAY(I2,I1)=(DEF(I2-1,I1)+DEF(I2,I1))/2.
C             KDEF=4: Array DEF should have dimensions DEF(N2,N1).
C                     Then the default is ARRAY(I2,I1)=DEF(I2,I1).
C     DEF...  Array to specify default values.
C
C Output:
C     ARRAY...Array of N2*N1 values.  The values are read from the file
C             which name is given by parameter NAME in SEP input data
C             read previously to the invocation of this subroutine.
C             Since the inner loop in the input file is over the
C             horizontal grid lines (N1 values), whereas the inner loop
C             in 2-D FD programs is over the vertical grid lines
C             (N2 values), N1*N2 matrix WORK read from the file is
C             transposed to yield N2*N1 output array ARRAY.
C     WORK... Temporary storage array containing the grid values in the
C             same order as in the input file.
C
C Date: 1998, November 11
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      CHARACTER*80 FILE
      CHARACTER*11 FORM
      INTEGER I1,I2,I,J
C
C.......................................................................
C
C     Name and form of the input data file:
      IF(NAME.EQ.' ') THEN
        FILE=' '
      ELSE
        CALL RSEP3T(NAME,FILE,' ')
        CALL RSEP3T('FORM',FORM,'FORMATTED')
      END IF
C
C     Check whether the filename is specified in the input SEP data:
      IF(FILE.EQ.' ') THEN
C
C       No input file specified:
        I=1
        J=1
        IF(KDEF.EQ.1) THEN
C         Averaging along the X1 axis:
          DO 12 I2=1,N2
            ARRAY(I)=DEF(J)
            DO 11 I1=2,N1-1
              I=I+N2
              J=J+N2
              ARRAY(I)=0.5*(DEF(J-1)+DEF(J))
   11       CONTINUE
            I=I+N2
            ARRAY(I)=DEF(J)
            I=I-N2*(N1-1)+1
            J=J-N2*(N1-2)+1
   12     CONTINUE
        ELSE IF(KDEF.EQ.2) THEN
C         Averaging along the X2 axis:
          DO 22 I1=1,N1
            ARRAY(I)=0.5*DEF(J)
            DO 21 I2=2,N2-1
              I=I+1
              J=J+1
              ARRAY(I)=0.5*(DEF(J-1)+DEF(J))
   21       CONTINUE
            I=I+1
            ARRAY(I)=DEF(J)
            I=I+1
            J=J+1
   22     CONTINUE
        ELSE IF(KDEF.EQ.4) THEN
C         Copying:
          DO 40 I=1,N2*N1
            ARRAY(I)=DEF(I)
   40     CONTINUE
        ELSE
C         Zeros:
          DO 90 I=1,N1*N2
            ARRAY(I)=0.
   90     CONTINUE
        END IF
C
C       Transposing the array:
*       CALL GMTRA(ARRAY,WORK,N2,N1)
C
      ELSE
C
C       Reading the array (undefined values are replaced by zeros):
        CALL RARRAY(LU,FILE,FORM,.TRUE.,0.,N1*N2,WORK)
C
C       Transposing the array:
        CALL GMTRA(WORK,ARRAY,N1,N2)
C
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE FDOUT1(LU,N1,N2,NAME,ARRAY,WORK)
C
      INTEGER LU,N1,N2
      CHARACTER*(*) NAME
      REAL ARRAY(N2*N1),WORK(N1*N2)
C
C Subroutine FDOUT1 to write the discrete values of a single quantity
C for 2-D finite difference programs.
C
C Input:
C     LU...   Logical unit number to be used for output.
C     N1...   Number of values in the horizontal direction.
C     N2...   Number of values in the vertical direction.
C     NAME... String containing the name of the parameter specifying the
C             name of the output file to be written.
C             Except for its case, it should match the parameter name in
C             a previously read input SEP parameter file.  If parameter
C             NAME is not defined in the input SEP parameter file, no
C             output is performed.
C     ARRAY...Array of N2*N1 values.  The values will be written to the
C             file which name is given by parameter NAME in SEP input
C             data read previously to the invocation of this subroutine.
C             Since the inner loop in the input file is over the
C             horizontal grid lines (N1 values), whereas the inner loop
C             in 2-D FD programs is over the vertical grid lines
C             (N2 values), N2*N1 matrix ARRAY is transposed to yield
C             N1*N2 array WORK written to the output file.
C     WORK... Temporary working array, its content will be destroyed.
C
C Output:
C     WORK... Undefined.
C
C Date: 1998, May 31
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      CHARACTER*80 FILE
      CHARACTER*11 FORM
C
C.......................................................................
C
C     Name of the output data file:
      IF(NAME.EQ.' ') THEN
C       The output file is already open:
        FILE=' '
      ELSE
        CALL RSEP3T(NAME,FILE,' ')
C       Check whether the filename is specified in the input SEP data:
        IF(FILE.EQ.' ') THEN
          RETURN
        END IF
      END IF
C
C     Form of the output data file:
      CALL RSEP3T('FORM',FORM,'FORMATTED')
C
C     Transposing the array:
      CALL GMTRA(ARRAY,WORK,N2,N1)
C
C     Writing the array:
      CALL WARRAY(LU,FILE,FORM,.FALSE.,0.,.FALSE.,0.,N1*N2,WORK)
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE FDGRID(KMAX,LMAX,KLMAX,NX,NZ,DX)
C
      INTEGER KMAX,LMAX,KLMAX,NX,NZ
      REAL DX
C
C Subroutine to read receiver positions for 2-D finite differences.
C
C Input:
C     KMAX... Maximum number of vertical gridlines.
C     LMAX... Maximum number of horizontal gridlines.
C     KLMAX...Maximum number of gridpoints.
C
C Output:
C     NX...   Number of vertical gridlines.
C     NZ...   Number of horizontal gridlines.
C     DX...   Absolute value of the grid interval.  The absolute values
C             of the horizontal and vertical grid intervals must be
C             equal.
C
C Date: 1998, September 11
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER N1,N2,N3
      REAL D1,D2,D3,DZ
C
C.......................................................................
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of parameter in input data, Variable, Default):
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      IF (N1.LE.1) THEN
        NX=N2
        NZ=N3
        DX=ABS(D2)
        DZ=ABS(D3)
      ELSE IF (N2.LE.1) THEN
        NX=N1
        NZ=N3
        DX=ABS(D1)
        DZ=ABS(D3)
      ELSE IF (N3.LE.1) THEN
        NX=N1
        NZ=N2
        DX=ABS(D1)
        DZ=ABS(D2)
      ELSE
C       FDAUX-01
        CALL ERROR('FDAUX-01: 3-D grid specified for 2-D FD')
C       Exactly one of input parameters N1, N2 ,N3 must equal 1 for 2-D
C       finite differences.
      END IF
      IF (NX.LE.1.OR.NZ.LE.1) THEN
C       FDAUX-02
        CALL ERROR('FDAUX-02: 1-D grid specified for 2-D FD')
C       Exactly one of input parameters N1, N2 ,N3 must equal 1 for 2-D
C       finite differences.
      END IF
      IF (NX.GT.KMAX) THEN
C       FDAUX-03
        CALL ERROR('FDAUX-03: Small dimension KMAX')
C       Number NX of vertical gridlines exceeds maximum number KMAX of
C       vertical gridlines.
      END IF
      IF (NZ.GT.LMAX) THEN
C       FDAUX-04
        CALL ERROR('FDAUX-04: Small dimension LMAX')
C       Number NZ of horizontal gridlines exceeds maximum number LMAX of
C       horizontal gridlines.
      END IF
      IF (NX*NZ.GT.KLMAX) THEN
C       FDAUX-05
        CALL ERROR('FDAUX-05: Small dimension KLMAX')
C       Number N1*N2*N3 of gridpoints exceeds maximum number KLMAX of
C       gridpoints.
      END IF
      IF (DX.NE.DZ) THEN
C       FDAUX-06
        CALL ERROR('FDAUX-06: Different grid intervals')
C       Horizontal and vertical grid intervals must be equal for
C       2-D finite-difference programs.
      END IF
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE FDREC(LU,MAXREC,NREC,REC,KRE,LRE,KSOUR,LSOUR)
C
      INTEGER LU,MAXREC,NREC,KRE(MAXREC),LRE(MAXREC),KSOUR,LSOUR
      CHARACTER*(*) REC(MAXREC)
C
C Subroutine to read receiver and source positions for 2-D finite
C differences.
C
C Input:
C     LU...   Logical unit number to be used for input.
C     MAXREC..Maximum number of receivers.
C
C Output:
C     NREC... Number of receivers.
C     KREC,LREC... Horizontal and vertical indices of receiver
C             gridpoints.
C     KSOUR,LSOUR... Horizontal and vertical indices of the source
C             gridpoint.
C
C Date: 1998, October 7
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL UNDEF
      PARAMETER (UNDEF=-999999.)
C
C     Auxiliary storage locations:
      CHARACTER*80 FREC
      CHARACTER*6  TEXT
      INTEGER N1,N3
      REAL D1,D3,O1,O3,X1,X2,X3
C
C.......................................................................
C
      CALL RSEP3I('N1' ,N1  ,1  )
      CALL RSEP3I('N3' ,N3  ,1  )
      CALL RSEP3R('D1' ,D1  ,1. )
      CALL RSEP3R('D3' ,D3  ,1. )
      CALL RSEP3R('O1' ,O1  ,0. )
      CALL RSEP3R('O3' ,O3  ,0. )
C
      NREC=0
      CALL RSEP3T('REC',FREC,' ')
      IF(FREC.NE.' ') THEN
        OPEN(LU,FILE=FREC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
   11   CONTINUE
          X1=UNDEF
          X2=UNDEF
          X3=UNDEF
          READ(LU,*,END=12) TEXT,X1,X2,X3
          IF(X1.EQ.UNDEF.AND.X2.EQ.UNDEF.AND.X3.EQ.UNDEF) THEN
            GO TO 12
          END IF
          IF(NREC.GE.MAXREC) THEN
C           FDAUX-07
            CALL ERROR('FDAUX-07: Small dimension MAXREC')
C           Dimension MRAM must be at least
C             19*N1*N3+10*N1+12*N3+(2*NTFD+5)*MAX(1,NREC),
C           where N1, N3 and NTFD are parameters in the main input file
C           and NREC is the number of receivers specified in the file
C           which name is given by parameter REC.
          END IF
          NREC=NREC+1
          REC(NREC)=TEXT
          IF(X1.EQ.UNDEF) X1=0.
          IF(X2.EQ.UNDEF) X2=0.
          IF(X3.EQ.UNDEF) X3=0.
          KRE(NREC)=NINT((X1-O1)/D1)+1
          LRE(NREC)=NINT((X3-O3)/D3)+1
          IF(KRE(NREC).LT.1.OR.KRE(NREC).GT.N1.OR.
     *       LRE(NREC).LT.1.OR.LRE(NREC).GT.N3) THEN
C            FDAUX-08
           CALL ERROR('FDAUX-08: Receiver point outside the FD grid')
          END IF
        GO TO 11
   12   CONTINUE
        CLOSE(LU)
      END IF
C
      CALL RSEP3T('SRC',FREC,' ')
      IF(FREC.NE.' ') THEN
        OPEN(LU,FILE=FREC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        X1=0.
        X2=0.
        X3=0.
        READ(LU,*) TEXT,X1,X2,X3
        CLOSE(LU)
        KSOUR=NINT((X1-O1)/D1)+1
        LSOUR=NINT((X3-O3)/D3)+1
        IF(KSOUR.LT.1.OR.KSOUR.GT.N1.OR.
     *     LSOUR.LT.1.OR.LSOUR.GT.N3) THEN
C          FDAUX-09
         CALL ERROR('FDAUX-09: Source point outside the FD grid')
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C