C 
C Program FD2D for 2-D P-SV elastic second-order finite differences
C
C Version: 5.20
C Date: 1998, November 11
C
C Program written by Jiri Zahradnik
C             Department of Geophysics, Charles University Prague
C             V Holesovickach 2, 180 00 Praha 8, Czech Republic
C             Tel.: +420-2-21912546, Fax: +420-2-21912555
C             E-mail: jz@karel.troja.mff.cuni.cz
C
C.......................................................................
C
C Contents:
C     Description of data files
C     General description
C     List of external procedures
C
C For the description of effective elastic parameters refer to file
C     modfd.for.
C
C.......................................................................
C                                                    
C Description of data files:
C
C Main input data read from the * input device:
C     The data are read in by the list directed input (free format), by
C     a single READ statement.  The data consist of character strings
C     which have to be enclosed in apostrophes.  Note that the
C     interactive * external unit may be piped or redirected to a file.
C     The following items are read:
C (1) 'FDSEP' /
C     'FDSEP'... String in apostrophes containing the name of the main
C             input data file with FD parameters in the SEP format.
C             The file specifies grid dimensions, FD time step, the form
C             and names of the files with effective elastic parameters
C             and snapshots, the type and position of the source,
C             the names of the files with receiver positions, snapshot
C             levels, the names of the files with synthetic seismograms
C             and snapshots, etc.
C     /...    It is recommended to terminate these data by a slash.
C     Default: 'FDSEP'='fd.h'.
C
C Data file 'FDSEP' has the form of the SEP (Stanford Exploration
C Project) 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 This program understands the following parameters common to nearly all
C programs dealing with SEP-like gridded data formats:
C     N1=positive integer... Number of gridpoints of the basic grid
C             along the X1 axis.  Note that 20 leftmost and 20 rightmost
C             gridlines are reserved for the dumpers and that an
C             explosive source requires 3 by 3 gridpoints.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints of the basic grid
C             along the X2 axis.  In 2-D, leave the default value.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints of the basic grid
C             along the X3 axis.  Note that the first (top) horizontal
C             gridline is assumed to be a free surface and that 20
C             bottom gridlines are reserved for the dumpers.
C             Default: N3=1
C     O1=real... X1 coordinate of the origin of the grid.
C             Default: O1=0.
C     O3=real... X3 coordinate of the origin of the grid.
C             Note that FD program 'fd2d.for' assumes the free Earth
C             surface at vertical coordinate X3=O3.
C             Default: O3=0.
C     D1=real... Grid spacing along the X1 axis.
C             The maximum relative error E in the wavefield is roughly
C               E = (2*pi*F)**3 * (D1/VSmin)**2 * Tmax / 24 ,
C             where VSmin is the minimum S-wave velocity in the model,
C             Tmax=(NTFD-1)*DTFD is the maximum time and F is the mean
C             frequency of the elastic waves, often roughly equal to the
C             reference frequency SIGF of the source time function.
C             Equation
C               D1 = VSmin * SQRT( 3 * E / Tmax / (pi*F)**3 )
C             may thus be used to estimate spatial steps D1 and D3.
C             Note that old-fashioned equation
C               D1 = VSmin / (10*Fmax) = VSmin / (30*SIGF)
C             corresponds to the maximum relative error
C               E = SIGF * Tmax * 1.15%
C             in the calculated wavefield.
C             Note that the grid intervals may also be negative.
C             Default: D1=1.
C     D3=real... Grid spacing along the X3 axis.
C             Absolute values of D1 and D3 must be equal.
C             Default: D3=1.
C Number and size of FD time steps:
C     DTFD=real... Time step for finite differences.
C             Choose DTFD as to satisfy the stability condition
C               DTFD < D1 / (1.65*VPmax) ,
C             close to the upper limit.
C             Here VPmax is the maximum P-wave velocity in the model.
C     NTFD=integer... Number of time levels.
C Other numerical parameters:
C     LA1=integer... It is recommended not to specify this parameter.
C             The medium is assumed homogeneous (in the 'bedrock') for
C             the grid lines L > LA1-2.  If there is no homogeneous
C             bedrock, choose the default of LA1=N3-1.
C             Default: LA1=N3-1
C Type and position of the source:
C     SRC='string'... Name of the input file with the position of the
C             point source.  If SRC is blank, the default coordinates
C             (0,0,0) of the source are used.
C             Description of file SRC
C             Default: SRC=' '.
C     KSIG=integer... Type of the source time function:
C             KSIG=1: Gabor signal.
C             KSIG=3: Kuepper (Mueller) signal.
C             Default: KSIG=0
C     SIGF=positive real... Reference frequency characterizing the
C             excitation wavelet.  The Kuepper (Mueller) wavelet has
C             time duration 1/SIGF.  For frequencies f > Fmax = 3*SIGF,
C             the spectral amplitude is negligibly small.
C             The excitation signal specifies the time dependence of the
C             density of the seismic force or the seismic moment along
C             a linear source.
C             To obtain the far-field displacement similar to that of
C             a point source, the excitation signal should specify the
C             time dependence of the halfth derivative of the point
C             seismic force or the point seismic moment.
C             Note that reference half-period TP=1/(2*SIGF) was used in
C             previous versions of this program instead of SIGF.
C             Default: SIGF=1.
C     SIGA=real... Reference amplitude of the excitation wavelet.
C             Maximum amplitude of the Kuepper wavelet for KSIG=3.
C             The amplitude describes the intensity of a linear seismic
C             force or moment.  The 3-D amplitude corresponding to
C             a point source is asymptotically 1/SQRT(2*pi*SIGMA)
C             greater.  Here SIGMA is the integral of the propagation
C             velocity along the raypath (in a homogeneous medium,
C             SIGMA=V*R, where V is the velocity and R the hypocentral
C             distance).  Thus, if we wish to obtain the amplitude at
C             hypocentral distance R roughly corresponding to point
C             source of intensity SIGA=A, we may choose the line source
C             of intensity SIGA=A/SQRT(2*pi*SIGMA), where SIGMA
C             corresponds to the wave type and hypocentral distance.
C             Default: SIGA=1.
C Files with receiver positions and synthetic seismograms:
C     REC='string'... Name of the input file with the positions of the
C             receivers.  If REC is blank, no receivers are considered.
C             Description of file REC
C             Default: REC=' '.
C     SS='string'... Name of the output GSE file with the synthetic
C             seismograms.  The file is created only if the receivers
C             are specified.
C             Default: SS='fd2d.gse'.
C     OUTU='string',OUTW='string'... Names of the output files with the
C             synthetic seismograms in the previously used format.  Just
C             for the backward compatibility with program 'trans.for'.
C             It is recommended not to specify them.
C               Files given by OUTU and OUTW, formerly named 'outu.dat'
C               and 'outw.dat', correspond to the horizontal (U) and
C               vertical (W) component, respectively. These output files
C               are  organized as 'time levels', not seismograms.  It
C               means that the computed displacement values for all
C               receivers at a single time level are recorded, and this
C               process is repeated subsequently for the all computed
C               time levels.  The output files can be processed by
C               program 'trans.for' which provides the output suitable
C               for plotting purposes (first column is time, second
C               column is the seismogram for the first receiver, etc.).
C               Attention: Only each 4th time level is on output.
C               The output time series have therefore the time step of
C               4*DTFD.  Files OUTU and OUTW may be tranformed by
C               'trans.for' into 'a.dat' and 'b.dat', respectively.
C               These files may serve for plotting by GRAPHER.
C             Defaults: OUTU=' ', OUTW=' '
C Output files with snapshots:
C     O1=real... Time of the first snapshot.
C             Default: O1=0.
C     D4=real... Time interval to store the snapshots.
C             Default: D4=DTFD
C     N4=integer... Number of the snapshots.
C             Default: N4=NTFD
C     SNAP1='string'... Name of the output file with the N1*N3*N4 grid
C             values of the horizontal component of the displacement.
C             The file is not created if SNAP1=' '.
C             Default: SNAP1=' '
C     SNAP3='string'... Name of the output file with the N1*N3*N4 grid
C             values of the vertical component of the displacement.
C             The file is not created if SNAP3=' '.
C             Default: SNAP3=' '
C Form of files with effective elastic parameters and snapshots:
C     FORM='string'... Form of the input and output files:
C             FORM='formatted': Formatted ASCII files,
C             FORM='unformatted': Unformatted files.
C Input files with effective elastic parameters:
C     A11='string', B11='string', C11='string'... Strings in apostrophes
C             containing the names of the input files with (N1-1)*N2*N3
C             values of the elastic parameters averaged in the direction
C             of the X1 axis, along the gridlines.
C             Defaults: A11=' ', B11=' ', C11=' '.
C     A33='string', B33='string', C33='string'... Strings in apostrophes
C             containing the names of the input files with N1*N2*(N3-1)
C             values of the elastic parameters averaged in the direction
C             of the X3 axis, along the gridlines.
C             Defaults: A33=' ', B33=' ', C33=' '.
C     A13='string', B13='string', C13='string'... Strings in apostrophes
C             containing the names of the input files with
C             (N1-1)*N2*(N3-1) values of the elastic parameters averaged
C             in the direction of the X1 axis, along the gridlines
C             shifted half a grid interval in the direction of the X3
C             axis.
C             Defaults: A13=' ', B13=' ', C13=' '.
C     A31='string', B31='string', C31='string'... Strings in apostrophes
C             containing the names of the input files with
C             (N1-1)*N2*(N3-1) values of the elastic parameters averaged
C             in the direction of the X3 axis, along the gridlines
C             shifted half a grid interval in the direction of the X1
C             axis.
C             Defaults: A31=' ', B31=' ', C31=' '.
C     Analogously A22, B22, C22, A12, B12, C12, A21, B21, C21, A23, B23,
C             C23, A32, B32 and C32.
C     DEN='string'... String in apostrophes containing the name of the
C             input file with N1*N2*N3 grid values of the density.
C             Default: DEN=' '
C Example of data set FDSEP
C
C                                                     
C Input file SRC containing the coordinates of the point source:
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) 'NAME',X1SRC,X2SRC,X3SRC,/
C     'NAME'... String reserved for the name of the initial point.
C             No meaning here, anything in apostrophes may be submitted.
C     X1SRC,X2SRC,X3SRC... Coordinates of a single initial point.
C     Default: X1SRC=0., X2SRC=0., X3SRC=0.
C (3) / (a slash) or end of file.
C Example of data set SRC
C
C                                                     
C Input file REC containing receiver coordinates:
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) Several times (2.1):
C (2.1) 'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR),/
C     'NAMER(IR)'... String reserved for the name of the receiver.
C             No meaning here, anything in apostrophes may be submitted.
C     X1REC(IR),X2REC(IR),X3REC(IR)... Coordinates of the receiver.
C     Default: X1REC(IR)=0, X2REC(IR)=0, X3REC(IR)=0.
C (3) / (a slash) or end of file.
C Example of data file REC
C
C.......................................................................
C                                                  
C General description:
C
C Calculation of P-SV waves in 2-D geological structures by the
C finite-difference method of the second order accuracy in both in space
C and time. No absorption included.
C
C The explicit PS2 scheme described in Zahradnik and Priolo (1995)
C is used. The uniform square grid: spatial step D1=D3.
C
C The computational region is a rectangle.
C The top side is the flat free Earth's surface.
C The left, right and bottom sides are nonreflecting after
C Emerman and Stephen. Simultaneously, artificial dumpers, 20*D1 wide,
C are applied at these three sides.
C
C The excitation is by a line source, described by Aboudi's body force
C or seismic moment corresponding to the explosive source.
C The excitation by a plane wave incident from below (well tested only
C for the vertical incidence), realized by the algorithm of Alterman and
C Karal, will be included later.
C
C Reference:
C     Zahradnik, J. and E. Priolo (1995): Heterogeneous
C     formulations of elastodynamic equations and finite-difference
C     schemes. Geophys. J. Int., 120, 663-676.
C
C If necessary, ask the author for more papers on the related subjects
C (absorption, filtration, convolution with earthquake
C ground motions, applications).
C
C Remarks:
C     Program FD2D should not be used to calculate elastic wavefields
C     in fluids!  Spurious waves, other than P waves, and instabilities
C     may occur in fluids.
C
C     The absorbing boundary conditions and dumpers do not work
C     perfectly!  The reflections from the boundaries may be found in
C     the seismograms by calculating the seismograms with two different
C     positions of the boundaries and comparing the results.
C
C     Mueller's time function
C     Program modification is possible for another time
C     function, see TIMFUN at the end of this code.
C
C     The line source is of the explosive (compressional) type. An easy
C     program modification is possible for another source type,
C     e.g., the vertical force, etc.
C
C     The nonreflecting boundaries after Emerman and Stephen are
C     applied.  Program modification is possible for changing the right
C     side to Reynolds.
C
C     The physical model is between the columns K=3 and K=NX-2.
C     When the right side is symmetrical, it is located at K=NX-2.
C     The columns K=1, 2, NX, NX-1 serve for the nonreflecting
C     (or symmetrical) sides.
C
C     The left, right and bottom sides = nonreflecting + dumpers.
C     Dumpers are 20 gridpoints wide.  Therefore, smaller models than
C     41(horiz)x21(vert) are not possible.
C     A program modification is possible for making the left or right
C     side symmetrical and/or to disable the dumpers.
C
C     It is useful to have some space, at least 5 gridlines,  between
C     the investigated structure and the nonreflecting sides.
C
C     The program has not been optimized. The main intention was to
C     permit many different models of the medium with no change in the
C     code.  Therefore, the largest part of the code form complicated
C     subroutines manipulating with the parameters of the medium.
C     The code has been originally developed as a 4th-order version.
C     Although the present version is 2nd-order only, many formal
C     complications from the planned 4th order still remain in the code.
C     Ignore numerous warnings 'formal argument ... never used'.
C
C.......................................................................
C                                                
C This file consists of the following external procedures:
C     FD2D... Main program allocating the memory for the FD code.
C             FD2D
C     FD...   Subroutine containing the main finite-difference code
C             (originally the main program).
C             FD
C     CDUMP...Function computes spatial variation of artificial dumpers.
C             CDUMP
C     AUXL1...Subroutine designed to transfer the displacements from
C             arrays UX,UY,etc., into auxiliary variables UOLD,U1,
C             U2,etc. (Selecting values from vicinity of a grid point
C             at line L=1.)
C             AUXL1
C     AUXL2...Subroutine designed to transfer the displacements from
C             arrays UX,UY,etc., into auxiliary variables UOLD,U1,
C             U2,etc. (Selecting values from vicinity of a grid point
C             at line L=2.)
C             AUXL2
C     AUXLA...Subroutine designed to transfer the displacements from
C             arrays UX,UY,etc., into auxiliary variables UOLD,U1,
C             U2,etc. (Selecting values from vicinity of a grid point
C             at line L>2.)
C             AUXLA
C     INDEX...Subroutine designed to compute 1-D index ICEN of a grid
C             point K,L, indices II(1)-II(25) of the grid-point
C             neighbours, indices isa,isb,...ish of the neighbouring
C             small legs (h/2), computing indices ila,ilb,...ilh of
C             the neighbouring large legs (h).
C             INDEX
C     PARL1...Subroutine designed to transfer the effective parameters
C             from arrays AHL,AHS,etc. into auxiliary variables.
C             (Selecting values from vicinity of a grid point at
C             line L=1.)
C             PARL1
C     PARL2...Subroutine designed to transfer the effective parameters
C             from arrays AHL,AHS,etc. into auxiliary variables.
C             (Selecting values from vicinity of a grid point at
C             line L=2.)
C             PARL2
C     PARL... Subroutine designed to transfer the effective parameters
C             from arrays AH,BH,etc. into auxiliary variables.
C             (Selecting values from vicinity of a grid point at
C             line L>2.)
C             PARL
C     PARHOM..Subroutine analogous to PARL but for a grid point with
C             a completely homogeneous vicinity.
C             PARHOM
C     DIFL1...Subroutine devoted to finite-difference scheme for a grid
C             point at the free surface at line L=1
C             DIFL1
C     DIFL2...Subroutine devoted to finite-difference scheme for a grid
C             point in a heterogeneous medium at line L=2.
C             DIFL2
C     DIFL... Subroutine devoted to finite-difference scheme for a grid
C             point in a heterogeneous medium at line L>2.
C             DIFL
C     DIFHOM..Subroutine devoted to finite-difference scheme for a grid
C             point in a homogeneous medium.
C             DIFHOM
C     FLERIB..Subroutine devoted to left and right nonreflecting
C             boundaries.
C             FLERIB
C     FBOB... Subroutine devoted to bottom nonreflecting boundaries.
C             FBOB
C     FCORNB..Subroutine devoted to left and right bottom corners of
C             the nonreflecting boundaries.
C             FCORNB
C     SAVER...Subroutine saving 'very old' values (from time level M-1)
C             at boundary lines for subsequent computations at
C             nonreflecting boundaries.
C             SAVER
C     TIMFUN..Function devoted to the excitation signal.
C             TIMFUN
C
C=======================================================================
C
C     
C
      PROGRAM FD2D
C
C Main program allocating the memory for the FD code.
C
C Subroutines and external functions required:
      EXTERNAL FD,RSEP1,RSEP3I
C     FD...   This file.
C     RSEP1,RSEP3I... File 'sep.for'.
C Note that the above subroutines reference many other external
C procedures from various subroutine files.  These indirectly
C referenced procedures are not named here, but are listed in the
C particular subroutine files.
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/ to allocate large array RAM:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C-----------------------------------------------------------------------
C
C     Input and output data filenames and logical unit numbers:
      INTEGER LU1,LU2,LU3,LU4
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4)
      CHARACTER*80 FSEP
C     LU1,LU2,LU3,LU4... Logical unit numbers to be used to read and
C             write data files.
C     FSEP... The name of the input data file with FD parameters in SEP
C             format.
C
      INTEGER MREC,N1,N3,N,NTFD,NR,NS,NGRID,I1,I2,I3,I4
      PARAMETER (MREC=1024)
      CHARACTER*6 REC(MREC)
C
C.......................................................................
C
C     Main input data:
C     Default:
      FSEP='fd.h'
C     Reading main input data:
      WRITE(*,'(A)') ' FD2D: Enter name of input SEP file [''fd.h'']: '
      READ (*,*) FSEP
      WRITE(*,'(A)') '+FD2D: Working...                               '
C
C     Reading all the data from input SEP parameter file:
      CALL RSEP1(LU1,FSEP)
C
C     Data specifying numbers of gridpoints and FD time steps:
      CALL RSEP3I('N1'  ,N1  ,1)
      CALL RSEP3I('N3'  ,N3  ,1)
      CALL RSEP3I('NTFD',NTFD,1)
C
C     Allocating the memory in array RAM:
      NGRID=19*N1*N3+10*N1+12*N3
      NS=2*NTFD+5
      IF(NGRID+NS.GT.MRAM) THEN
C       FD2D-01
        CALL ERROR('FD2D-01: Small dimension MRAM in ram.inc')
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 and
C       NREC is the number of receivers specified in the file which name
C       is given by parameter REC.
      END IF
      N =N1*N3
      NR=MIN0((MRAM-NGRID)/NS,MREC)
      NS=NTFD*NR
      I1= 1+19*N
      I2=I1+10*N1
      I3=I2+12*N3
      I4=I3+ 5*NR
C
      CALL FD(LU1,LU2,LU3,LU4,N1,N3,N,NR,NS,
     *        RAM(       1),RAM(   N + 1),RAM( 2*N + 1),RAM( 3*N + 1),
     *        RAM( 4*N + 1),RAM( 5*N + 1),RAM( 6*N + 1),RAM( 7*N + 1),
     *        RAM( 8*N + 1),RAM( 9*N + 1),RAM(10*N + 1),RAM(11*N + 1),
     *        RAM(12*N + 1),RAM(13*N + 1),RAM(14*N + 1),RAM(15*N + 1),
     *        RAM(16*N + 1),RAM(17*N + 1),RAM(18*N + 1),
     *        RAM(      I1),RAM(   N1+I1),RAM( 2*N1+I1),RAM( 3*N1+I1),
     *        RAM( 4*N1+I1),RAM( 5*N1+I1),RAM( 6*N1+I1),
     *        RAM(      I2),RAM(   N3+I2),RAM( 2*N3+I2),RAM( 3*N3+I2),
     *        RAM( 4*N3+I2),RAM( 5*N3+I2),RAM( 6*N3+I2),RAM( 7*N3+I2),
     *        RAM( 8*N3+I2),RAM( 9*N3+I2),RAM(10*N3+I2),RAM(11*N3+I2),
     *        IRAM(     I3),IRAM(  NR+I3),IRAM(2*NR+I3),
     *        RAM( 3*NR+I3),RAM( 4*NR+I3),
     *        RAM(      I4),RAM(   NS+I4),REC)
C
      WRITE (*,'(A)') '+FD2D: Done.                                    '
      STOP
      END
C
C
C=======================================================================
C
C     
C
      SUBROUTINE FD(LU1,LU2,LU3,LU4,KMAX,LMAX,KLMAX,MAXREC,MAXSS,
     *              UX,UY,WX,WY,AHS,BHS,CHS,AVS,BVS,CVS,
     *                          AHL,BHL,CHL,AVL,BVL,CVL,DEN,DUMP,WORK,
     *              UNZ1,UNZ2,UNZ,WNZ1,WNZ2,WNZ,EX,
     *              UNX1,UNX2,UNX,WNX1,WNX2,WNX,UK1,UK2,UK3,WK1,WK2,WK3,
     *              KRE,LRE,INREC,ULEV,WLEV,USS,WSS,REC)
      INTEGER LU1,LU2,LU3,LU4,KMAX,LMAX,KLMAX,MAXREC,MAXSS
      REAL UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX)
      REAL AHS(KLMAX),BHS(KLMAX),CHS(KLMAX)
      REAL AVS(KLMAX),BVS(KLMAX),CVS(KLMAX)
      REAL AHL(KLMAX),BHL(KLMAX),CHL(KLMAX)
      REAL AVL(KLMAX),BVL(KLMAX),CVL(KLMAX)
      REAL DEN(KLMAX),DUMP(KLMAX),WORK(KLMAX)
      REAL UNZ1(KMAX),UNZ2(KMAX),UNZ(KMAX)
      REAL WNZ1(KMAX),WNZ2(KMAX),WNZ(KMAX)
      REAL EX(4,KMAX)
      REAL UNX1(LMAX),UNX2(LMAX),UNX(LMAX)
      REAL WNX1(LMAX),WNX2(LMAX),WNX(LMAX)
      REAL UK1(LMAX),UK2(LMAX),UK3(LMAX)
      REAL WK1(LMAX),WK2(LMAX),WK3(LMAX)
      INTEGER KRE(MAXREC),LRE(MAXREC),INREC(MAXREC)
      REAL ULEV(MAXREC),WLEV(MAXREC)
      REAL USS(MAXSS),WSS(MAXSS)
      CHARACTER*(*) REC(MAXREC)
C
C Subroutine (originally main program) reading the input data and
C computing P-SV waves in 2-D geological structures by the
C finite-difference method.
C
C Input:
C     LU1,LU2,LU3,LU4... Logical unit numbers to be used to read and
C             write data files.
C     KMAX,LMAX,KLMAX,MAXREC,MAXSS... Array dimensions.
C     AHL,BHL,CHL,
C     AVL,BVL,CVL,
C     AHS,BHS,CHS,
C     AVS,BVS,CVS,
C     DEN,QFC...arrays with effective parameters
C     WORK...   temporary working array
C     REC...
C
C No output.
C
C Subroutines and external functions required:
      EXTERNAL INDEX
      EXTERNAL RSEP1,RSEP3I,RSEP3R
      EXTERNAL FDGRID,FDOUT1,FDREC,FDIN
C     INDEX...This file.
C     RSEP1,RSEP3I,RSEP3R...File 'sep.for'.
C     FDGRID,FDOUT1,FDREC,FDIN..File 'fdaux.for'.
C Note that the above subroutines reference many other external
C procedures from various subroutine files.  These indirectly
C referenced procedures are not named here, but are listed in the
C particular subroutine files.
C
C-----------------------------------------------------------------------
C
C     Storage locations:
C
C     Input and output data filenames:
      CHARACTER*80 FILNA1,FILNA2,OUTSS,FOUTU,FOUTW
      CHARACTER*11 FORM
C
      DIMENSION II(25),IJ(25)
      DIMENSION INU(6),INW(6),UP(6),WP(6)
C
      COMMON /DSPLC/ UOLD,WOLD,
     1           U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14,
     2           U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25,
     3           W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14,
     4           W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
      COMMON /DMP/ NZONE
C
C.......................................................................
C
C     Data specifying FD grid dimensions
      CALL FDGRID(KMAX,LMAX,KLMAX,NX,NZ,DX)
C     Data required to write synthetic seismograms
      CALL RSEP3R('O1',O1,0.0)
      CALL RSEP3R('O3',O3,0.0)
      CALL RSEP3R('D1',D1,1.0)
      CALL RSEP3R('D3',D3,1.0)
C
C     Number and size of FD time steps:
      CALL RSEP3I('NTFD',NT,1)
      CALL RSEP3R('DTFD',DT,1.)
C
C     Type and position of the source:
C     (arguments: Name of parameter in input data, Variable, Default):
C     The excitation (plane waves P=0, S=1; line source=2):
C     (for future extension, do not specify in input data, use ISOURC=2)
      CALL RSEP3I('ISOURC',ISOURC,2)
C     Grid indices of the source:
C     (for backward compatibility, do not specify in input data)
      CALL RSEP3I('KSOUR',KSOUR,1)
      CALL RSEP3I('LSOUR',LSOUR,1)
C     Semiduration TP of the excitation wavelet:
      CALL RSEP3R('SIGF',SIGF,1.)
      TP=0.5/SIGF
C     Medium must be homogeneous for L.GE.LA1-2. Choose LA1.GE.4.
C     If a plane wave is specified, wavefront t=0 is at L=LA1+3 and
C     L=LA1 to LA1+3 are prescribed.
      CALL RSEP3I('LA1',LA1,NZ-1)
C
C     Reading receiver (and source) positions:
      CALL FDREC(LU1,MAXREC,NREC,REC,KRE,LRE,KSOUR,LSOUR)
      IF(NT*NREC.GT.MAXSS) THEN
C       FD2D-02
        CALL ERROR('FD2D-02: Small dimension MAXSS')
      END IF
C
C     Reading data for snapshot levels:
      CALL RSEP3I('N4',N4,NT)
      CALL RSEP3R('D4',D4,DT)
      CALL RSEP3R('O4',O4,0.)
      IF(D4.LT.DT) THEN
C       FD2D-03
        CALL ERROR('FD2D-03: Small time step D4 for snapshot levels')
      END IF
      IF(O4+FLOAT(N4-1)*D4.GE.(FLOAT(NT)-0.5)*DT) THEN
C       FD2D-04
        CALL ERROR('FD2D-04: Large number N4 of snapshot levels')
      END IF
C
C     Reading the gridded effective elastic parameters:
      CALL FDIN(LU1,NX,NZ,AHL,BHL,CHL,AVL,BVL,CVL,
     *                    AHS,BHS,CHS,AVS,BVS,CVS,DEN,UX,UY)
C
C     The right, left and bottom sides of the comp. region are
C     nonreflecting when NSYM=0.  The right and bottom are
C     nonreflecting but the left is symmetrical when NSYM=1.
C     The latter is available for the incident S-wave only.
      CALL RSEP3I('NSYM'  ,NSYM  ,0)
C     The incidence angle, positive clockwise (degrees).
C     For vertical incidence the angle is 0.
      CALL RSEP3R('ANGINC',ANGINC,0.)
C     Azimuth of the profile with respect to the incident ray (degrees).
c     For the ray in (x,z) plane choose 0.
      CALL RSEP3R('AZIPRO',AZIPRO,0.)
C     Incident-wave velocity (m/s) in bedrock:
      CALL RSEP3R('ROCKVE',ROCKVE,999999.)
C
C     Opening output receiver files (old form):
      IF(NREC.GT.0) THEN
        CALL RSEP3T('OUTU',FOUTU,' ')
        IF(FOUTU.NE.' ') THEN
          OPEN(LU3,FILE=FOUTU)
        END IF
        CALL RSEP3T('OUTW',FOUTW,' ')
        IF(FOUTW.NE.' ') THEN
          OPEN(LU4,FILE=FOUTW)
        END IF
      END IF
C
C     Opening output snapshot files:
C     Form of the output files:
      CALL RSEP3T('FORM',FORM,'FORMATTED')
C     Names of the output snapshot files:
      CALL RSEP3T('SNAP1',FILNA1,' ')
      CALL RSEP3T('SNAP3',FILNA2,' ')
      IF(FILNA1.NE.' ') THEN
        OPEN(LU1,FILE=FILNA1,FORM=FORM)
      END IF
      IF(FILNA2.NE.' ') THEN
        OPEN(LU2,FILE=FILNA2,FORM=FORM)
      END IF
C
C     Auxiliary parameters
C     TSWI2=9999.
C
      DO 350 I=1,NREC
        INREC(I)=LRE(I)+(KRE(I)-1)*NZ
  350 CONTINUE
C
C     Preparatory calculations for artificial dumpers
C     fixed with NZONE=20.
C
C     goto 4001
      NZONE=20
C
      NXV=NX-NZONE
      NZV=NZ-NZONE
      NXP=NXV+1
      NZP=NZV+1
      DO 3200 K=1,NX
        DO 3100 L=1,NZ
          IDU=L+(K-1)*NZ
          DUMP(IDU)=1.
 3100 CONTINUE
 3200 CONTINUE
C
C     Change to K=1,NXV, and 'C' the loop 3700 if symmetrical left side.
C
      DO 3400 K=NZONE+1,NXV
        DO 3300 L=NZP,NZ
          IDU=L+(K-1)*NZ
          DUMP(IDU)=CDUMP(NZ+1-L)
 3300 CONTINUE
 3400 CONTINUE
      DO 3700 K=1,NZONE
        DO 3500 L=1,NZV
          IDU=L+(K-1)*NZ
          DUMP(IDU)=CDUMP(K)
 3500   CONTINUE
        DO 3600 L=NZP,NZ
          IDU=L+(K-1)*NZ
          DUMP(IDU)=CDUMP(K)*CDUMP(NZ+1-L)
 3600   CONTINUE
 3700 CONTINUE
C
      DO 4000 K=NXP,NX
        DO 3800 L=1,NZV
          IDU=L+(K-1)*NZ
          DUMP(IDU)=CDUMP(NX+1-K)
 3800   CONTINUE
        DO 3900 L=NZP,NZ
          IDU=L+(K-1)*NZ
          DUMP(IDU)=CDUMP(NX+1-K)*CDUMP(NZ+1-L)
 3900   CONTINUE
 4000 CONTINUE
C
 4001 CONTINUE
      NZ5=NZ-5
      NZ6=2*NZ+1
      CN=DT/DX
      CNST=CN*CN/12.
      NSIZE=NX*NZ
      LA2=LA1+1
      LA3=LA1+2
      LA4=LA1+3
C
      NSIZEA=NSIZE
      NXA=NX
      NZA=NZ
C
      IF(ISOURC.NE.2) THEN
        PI=3.141593
        ANG=PI*ANGINC/180.
        AZI=PI*AZIPRO/180.
        LZERO=LA3
        IF(ANG.GE.0.) KZERO=1
        IF(ANG.LT.0.) KZERO=NX
        SX=SIN(ANG)*COS(AZI)/ROCKVE
        SY=SIN(ANG)*SIN(AZI)/ROCKVE
        SZ=COS(ANG)/ROCKVE
        TIMEXC=SX*FLOAT(NX-1)*DX+SZ*3.*DX + 2.*TP
C
        IF(ISOURC.EQ.1) GOTO 353
C       The P-wave incidence
        UCOM=SIN(ANG)
        WCOM=-1.*COS(ANG)
        GOTO 354
C       The S-wave incidence
  353   UCOM=-1.*COS(ANG)
        WCOM=-1.*SIN(ANG)
  354   CONTINUE
      ELSE
C       Explosive point source
        UCOM=0.
        WCOM=0.
        TIMEXC=0.
        IF(KSOUR.LT.2.OR.KSOUR.GT.NX-1.OR.
     *     LSOUR.LT.2.OR.LSOUR.GT.NZ-1) THEN
C         FD2D-05
          CALL ERROR('FD2D-05: Explosive source exceeds the grid')
C         Explosive source is composed of 3*3 gridpoints.  It thus
C         cannot be centred at the boundary of the grid.  It is also not
C         recommended to place the source into a dumper.  The dumpers
C         are located along the bottom and sides of the grid and are
C         NZONE gridpoints thick (NZONE=20 in this version).
        END IF
        ICSR=LSOUR+(KSOUR-1)*NZ
        INU(1)=ICSR-NZ-1
        INU(2)=ICSR-NZ
        INU(3)=ICSR-NZ+1
        INU(4)=ICSR+NZ-1
        INU(5)=ICSR+NZ
        INU(6)=ICSR+NZ+1
        INW(1)=ICSR-NZ-1
        INW(2)=ICSR-1
        INW(3)=ICSR+NZ-1
        INW(4)=ICSR-NZ+1
        INW(5)=ICSR+1
        INW(6)=ICSR+NZ+1
C       Unit seismic moment
        CA=(DT/DX)**2/DEN(ICSR)/(4.*DX)
        CB=0.5*CA
        UP(1)=CB
        UP(2)=CA
        UP(3)=CB
        UP(4)=-CB
        UP(5)=-CA
        UP(6)=-CB
        WP(1)=CB
        WP(2)=CA
        WP(3)=CB
        WP(4)=-CB
        WP(5)=-CA
        WP(6)=-CB
      ENDIF
C
      DO 360 I=1,NSIZE
        UX(I)=0.
        UY(I)=0.
        WX(I)=0.
        WY(I)=0.
  360 CONTINUE
C
C.......................................................................
C
C     Starting loop over time levels
C
      DO 2000 M=1,NT
      WRITE(*,'(3(A,I6))') '+',M,' of',NT,' time levels'
      MPAR=1
      IF((FLOAT(M)/2.-FLOAT(M/2)).LT.0.1) MPAR=2
C
C      Excitation signal at lines L=LA1 to LA4
C
      T=FLOAT(M)*DT
      DO 650 I=1,4
      DO 640 K=1,NX
      EX(I,K)=0.
  640 CONTINUE
  650 CONTINUE
      IF(T.GT.TIMEXC) GOTO 900
      TRET=T-DT
C
C      We assume the incident wavefront striking the plane (x,z) at t=0
      TAUY=SY*0.
C
      DO 800 L=LA1,LA4
      I=(L-LA1)+1
      TAUZ=SZ*FLOAT(LZERO-L)*DX
      DO 700 K=1,NX
      TAUX=SX*FLOAT(K-KZERO)*DX
      TAU=TAUX+TAUY+TAUZ
      EX(I,K)=TIMFUN(TRET,TAU,TP)
  700 CONTINUE
  800 CONTINUE
  900 CONTINUE
C
C      Saving boundary-strip values from the 'very old' time level M-1
C
      IF(MPAR.EQ.1) THEN
       CALL SAVER(KMAX,LMAX,KLMAX,UY,WY,
     1            UK1,UK2,UK3,UNX2,UNX1,UNX,UNZ2,UNZ1,UNZ,
     2            WK1,WK2,WK3,WNX2,WNX1,WNX,WNZ2,WNZ1,WNZ)
      ELSE
       CALL SAVER(KMAX,LMAX,KLMAX,UX,WX,
     1            UK1,UK2,UK3,UNX2,UNX1,UNX,UNZ2,UNZ1,UNZ,
     2            WK1,WK2,WK3,WNX2,WNX1,WNX,WNZ2,WNZ1,WNZ)
      ENDIF
C
C      Computing 'new' time level M+1
C
C      Regular grid points (i.e. all except nonreflecting boundaries)
C
C      Line L=1
C
      L=1
      DO 1100 K=3,NX-2
      CALL INDEX(K,L)
      IF(MPAR.EQ.1) THEN
       CALL AUXL1(KLMAX,UX,UY,WX,WY)
      ELSE
       CALL AUXL1(KLMAX,UY,UX,WY,WX)
      ENDIF
      CALL PARL1(KLMAX,AHS,BVS,BVS,CHS,AHL,BVL,BVL,CHL,DEN)
      CALL DIFL1(UOLD,U3,U8,U13,U14,U15,U18,U23,
     1           W3,W5,W8,W9,W13,W14,W15,W18,W19,W23,W25,
     2           UNEW)
      CALL PARL1(KLMAX,BHS,CVS,AVS,BHS,BHL,CVL,AVL,BHL,DEN)
      CALL DIFL1(WOLD,W3,W8,W13,W14,W15,W18,W23,
     1           U3,U5,U8,U9,U13,U14,U15,U18,U19,U23,U25,
     2           WNEW)
C
      IF(MPAR.EQ.1) THEN
       UY(ICEN)=UNEW
       WY(ICEN)=WNEW
      ELSE
       UX(ICEN)=UNEW
       WX(ICEN)=WNEW
      ENDIF
 1100 CONTINUE
C
C      Line L=2
C
      L=2
      DO 1200 K=3,NX-2
      CALL INDEX(K,L)
      IF(MPAR.EQ.1) THEN
       CALL AUXL2(KLMAX,UX,UY,WX,WY)
      ELSE
       CALL AUXL2(KLMAX,UY,UX,WY,WX)
      ENDIF
      CALL PARL2(KLMAX,AHS,BVS,BVS,CHS,AHL,BVL,BVL,CHL,DEN)
      CALL DIFL2(UOLD,U3,U8,U12,U13,U14,U15,U18,
     1           U23,
     2           W3,W5,W7,W8,W9,W12,W13,W14,W15,W17,W18,W19,
     3           W23,W25,
     4           UNEW)
      CALL PARL2(KLMAX,BHS,CVS,AVS,BHS,BHL,CVL,AVL,BHL,DEN)
      CALL DIFL2(WOLD,W3,W8,W12,W13,W14,W15,W18,
     1           W23,
     2           U3,U5,U7,U8,U9,U12,U13,U14,U15,U17,U18,U19,
     3           U23,U25,
     4           WNEW)
C
      IF(MPAR.EQ.1) THEN
       UY(ICEN)=UNEW
       WY(ICEN)=WNEW
      ELSE
       UX(ICEN)=UNEW
       WX(ICEN)=WNEW
      ENDIF
 1200 CONTINUE
C
C      Lines L=3 to LA1-1 (heterogeneous, in general)
C
      DO 1400 L=3,LA1-1
      DO 1300 K=3,NX-2
c
      NOUT=0
c     if(l.eq.3.or.l.eq.4) then
c     NOUT=1
c     write(33,*) 'k,l',k,l
c     endif
c
      CALL INDEX(K,L)
      IF(MPAR.EQ.1) THEN
       CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UX,UY,WX,WY,T,
     1            TIMEXC,UCOM,WCOM)
      ELSE
       CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UY,UX,WY,WX,T,
     1            TIMEXC,UCOM,WCOM)
      ENDIF
      CALL PARL(KLMAX,AHS,BVS,BVS,CHS,AHL,BVL,BVL,CHL,DEN)
 1260 CALL DIFL(UOLD,U3,U8,U11,U12,U13,U14,U15,U18,
     1          U23,
     2          W1,W3,W5,W7,W8,W9,W11,W12,W13,W14,W15,W17,W18,
     3          W19,W21,W23,W25,
     4          UNEW,NOUT)
      CALL PARL(KLMAX,BHS,CVS,AVS,BHS,BHL,CVL,AVL,BHL,DEN)
      CALL DIFL(WOLD,W3,W8,W11,W12,W13,W14,W15,W18,
     1          W23,
     2          U1,U3,U5,U7,U8,U9,U11,U12,U13,U14,U15,U17,U18,
     3          U19,U21,U23,U25,
     4          WNEW,NOUT)
C
 1290 IF(MPAR.EQ.1) THEN
       UY(ICEN)=UNEW
       WY(ICEN)=WNEW
      ELSE
       UX(ICEN)=UNEW
       WX(ICEN)=WNEW
      ENDIF
 1300 CONTINUE
 1400 CONTINUE
C
C      Lines LA1 to NZ-2 (homogeneous), optional excitation at LA1-LA4
C
      DO 1600 L=LA1,NZ-2
      DO 1500 K=3,NX-2
      CALL INDEX(K,L)
      IF(MPAR.EQ.1) THEN
       CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UX,UY,WX,WY,T,
     1            TIMEXC,UCOM,WCOM)
      ELSE
       CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UY,UX,WY,WX,T,
     1            TIMEXC,UCOM,WCOM)
      ENDIF
      CALL PARHOM(KLMAX,AHL,BVL,BVL,CHL,DEN)
      CALL DIFHOM(UOLD,U3,U8,U11,U12,U13,U14,U15,U18,
     1            U23,
     2            W1,W5,W7,W9,W17,
     3            W19,W21,W25,
     4            UNEW)
      CALL PARHOM(KLMAX,BHL,CVL,AVL,BHL,DEN)
      CALL DIFHOM(WOLD,W3,W8,W11,W12,W13,W14,W15,W18,
     1            W23,
     2            U1,U5,U7,U9,U17,
     3            U19,U21,U25,
     4            WNEW)
C
C
      IF(MPAR.EQ.1) THEN
       UY(ICEN)=UNEW
       WY(ICEN)=WNEW
      ELSE
       UX(ICEN)=UNEW
       WX(ICEN)=WNEW
      ENDIF
 1500 CONTINUE
 1600 CONTINUE
C
C
C      Nonreflecting boundaries - inner lines
C
C      The right side E-S:
      CALL FLERIB(KLMAX,LMAX,1,NZA-3,NSIZEA-3*NZA,NSIZEA-2*NZA,
     1          NSIZEA-3*NZA,UNX2,UNX1,WNX2,WNX1,
     2          UX,UY,WX,WY,AHL,BHL,DEN)
C      The right side Reynolds
C     CALL FLERIB(KLMAX,LMAX,1,NZA-3,NSIZEA-3*NZA,NSIZEA-2*NZA,
C    1          NSIZEA-3*NZA,UNX,UNX2,WNX,WNX2,
C    2          UX,UY,WX,WY,AHL,BHL,DEN)
C      The left side E-S:
      CALL FLERIB(KLMAX,LMAX,1,NZA-3,2*NZA,NZA,
     1          NZA,UK3,UK2,WK3,WK2,
     2          UX,UY,WX,WY,AHL,BHL,DEN)
C      The bottom side E-S:
       CALL FBOB(KLMAX,KMAX,4,NXA-3,NZA-2,1,UNZ2,UNZ1,WNZ2,WNZ1,
     1           UX,UY,WX,WY,AVL,BVL,DEN)
C      The right bottom corner E-S:
       CALL FCORNB(KLMAX,NSIZEA-NZA-1,NZA,NZA,0,
     1             UX,UY,WX,WY,AHL,BHL,DEN)
C      The left bottom corner E-S:
       CALL FCORNB(KLMAX,2*NZA-1,0,-NZA,1,
     1             UX,UY,WX,WY,AHL,BHL,DEN)
C
C      Nonreflecting boundaries - outer lines
C
C      The right side E-S:
      CALL FLERIB(KLMAX,LMAX,1,NZA-2,NSIZEA-2*NZA,NSIZEA-NZA,
     1            NSIZEA-2*NZA,UNX1,UNX,WNX1,WNX,
     2            UX,UY,WX,WY,AHL,BHL,DEN)
C      The right side Reynolds:
C     CALL FLERIB(KLMAX,LMAX,1,NZA-2,NSIZEA-2*NZA,NSIZEA-NZA,
C    1            NSIZEA-2*NZA,UNX2,UNX1,WNX2,WNX1,
C    2            UX,UY,WX,WY,AHL,BHL,DEN)
C      The left side E-S:
      CALL FLERIB(KLMAX,LMAX,1,NZA-2,NZA,0,
     1            0,UK2,UK1,WK2,WK1,
     2            UX,UY,WX,WY,AHL,BHL,DEN)
C      The bottom side E-S:
       CALL FBOB(KLMAX,KMAX,3,NXA-2,NZA-1,0,UNZ1,UNZ,WNZ1,WNZ,
     1           UX,UY,WX,WY,AVL,BVL,DEN)
C      The right bottom corner E-S:
       CALL FCORNB(KLMAX,NSIZEA,NZA,NZA,0,
     1             UX,UY,WX,WY,AHL,BHL,DEN)
C      The left bottom corner E-S:
       CALL FCORNB(KLMAX,NZA,0,-NZA,1,
     1             UX,UY,WX,WY,AHL,BHL,DEN)
C
      IF(NSYM.EQ.0) goto 4300
C
C   Symmetrical bottom
c4200 DO 1624 L=NZ-1,NZ
c     DO 1623 K=3,NX-2
c     IND=L+(K-1)*NZ
c     IF(L.EQ.NZ-1) GOTO 1622
c     IF(MPAR.EQ.1) THEN
c     UY(IND)=UY(IND-4)
c     WY(IND)=-WY(IND-4)
c     ELSE
c     UX(IND)=UX(IND-4)
c     WX(IND)=-WX(IND-4)
c     ENDIF
c     GOTO 1623
c1622 IF(MPAR.EQ.1) THEN
c     UY(IND)=UY(IND-2)
c     WY(IND)=-WY(IND-2)
c     ELSE
c     UX(IND)=UX(IND-2)
c     WX(IND)=-WX(IND-2)
c     ENDIF
c1623 CONTINUE
c1624 CONTINUE
C
C   Symmetrical left side     ( GOOD FOR INCIDENT S WAVE ONLY !!!!!!!!)
      DO 1620 K=1,2
      DO 1610 L=1,NZ
      IND=L+(K-1)*NZ
      IF(K.EQ.2) GOTO 1605
      IF(MPAR.EQ.1) THEN
      UY(IND)=UY(IND+4*NZ)
      WY(IND)=-WY(IND+4*NZ)
      ELSE
      UX(IND)=UX(IND+4*NZ)
      WX(IND)=-WX(IND+4*NZ)
      ENDIF
      GOTO 1610
 1605 IF(MPAR.EQ.1) THEN
      UY(IND)=UY(IND+2*NZ)
      WY(IND)=-WY(IND+2*NZ)
      ELSE
      UX(IND)=UX(IND+2*NZ)
      WX(IND)=-WX(IND+2*NZ)
      ENDIF
 1610 CONTINUE
 1620 CONTINUE
c improved left side symmetry    GOOD FOR INCIDENT S WAVE ONLY !
c                          For S incidence change u=0 for w=0 !!!!!!!!!!
      DO 1627 L=NZ-1,NZ
      K=3
      IND=L+(K-1)*NZ
      IF(MPAR.EQ.1) THEN
      WY(IND)=0.
      ELSE
      WX(IND)=0.
      ENDIF
 1627 CONTINUE
C
C  Symmetrical right side      (GOOD FOR INCIDENT S WAVE ONLY !!!!!!!!)
c     DO 1640 K=NX-1,NX
c     DO 1630 L=1,NZ
c     IND=L+(K-1)*NZ
c     IF(K.EQ.NX-1) GOTO 1625
c     IF(MPAR.EQ.1) THEN
c     UY(IND)=UY(IND-4*NZ)
c     WY(IND)=-WY(IND-4*NZ)
c     ELSE
c     UX(IND)=UX(IND-4*NZ)
c     WX(IND)=-WX(IND-4*NZ)
c     ENDIF
c     GOTO 1630
c1625 IF(MPAR.EQ.1) THEN
c     UY(IND)=UY(IND-2*NZ)
c     WY(IND)=-WY(IND-2*NZ)
c     ELSE
c     UX(IND)=UX(IND-2*NZ)
c     WX(IND)=-WX(IND-2*NZ)
c     ENDIF
c1630 CONTINUE
c1640 CONTINUE
c improved right side symmetry    GOOD FOR INCIDENT S WAVE ONLY !
c                             For S incid change u=0 for w=0 !!!!!!!!!!
c     DO 1628 L=NZ-1,NZ
c     K=NX-2
c     IND=L+(K-1)*NZ
c     IF(MPAR.EQ.1) THEN
c     WY(IND)=0.
c     ELSE
c     WX(IND)=0.
c     ENDIF
c1628 CONTINUE
C
C
 4300 CONTINUE
C
C      Line source (a body force equivalent)
C
      IF(ISOURC.NE.2) GOTO 1700
      IF(T.GT.(2.*TP)) GOTO 1700
      TIMF=TIMFUN(T,0.,TP)
      IF(MPAR.EQ.1) THEN
C
C     Unit seismic force
*      CA=(DT/DX)**2/DEN(ICSR)
C     Single horizontal force
*      UY(ICSR)=UY(ICSR)-CA*TIMF
C     Single vertical force
*      WY(ICSR)=WY(ICSR)-CA*TIMF
C
      DO 1650 I=1,6
CC     Centre of compression (explosive line source)
       UY(INU(I))=UY(INU(I))-UP(I)*TIMF
       WY(INW(I))=WY(INW(I))-WP(I)*TIMF
CC     Centre of rotation
CC      UY(INW(I))=UY(INW(I))-WP(I)*TIMF
CC      WY(INU(I))=WY(INU(I))+UP(I)*TIMF
CC     Double couple
CC      UY(INW(I))=UY(INW(I))-WP(I)*TIMF
CC      WY(INU(I))=WY(INU(I))-UP(I)*TIMF
CC     Single horizontal force (extended)
CC      UY(INU(I))=UY(INU(I))-UP(I)*TIMF
 1650 CONTINUE
      ELSE
C
C      Single horizontal force
*       UX(ICSR)=UX(ICSR)-CA*TIMF
C      Single vertical force
*       WX(ICSR)=WX(ICSR)-CA*TIMF
C
       DO 1660 I=1,6
CC     Centre of compression
       UX(INU(I))=UX(INU(I))-UP(I)*TIMF
       WX(INW(I))=WX(INW(I))-WP(I)*TIMF
CC     Centre of rotation
CC      UX(INW(I))=UX(INW(I))-WP(I)*TIMF
CC      WX(INU(I))=WX(INU(I))+UP(I)*TIMF
CC     Double couple
CC      UX(INW(I))=UX(INW(I))-WP(I)*TIMF
CC      WX(INU(I))=WX(INU(I))-UP(I)*TIMF
CC     Single horizontal force (extended)
CC      UX(INU(I))=UX(INU(I))-UP(I)*TIMF
 1660 CONTINUE
      ENDIF
 1700 CONTINUE
C
C
C     Saving results for receivers:
C
      IF(NREC.GT.0) THEN
        DO 1800 I=1,NREC
          J=(I-1)*NT+M
          IF(MPAR.EQ.1) THEN
            USS(J)=UY(INREC(I))
            WSS(J)=WY(INREC(I))
          ELSE
            USS(J)=UX(INREC(I))
            WSS(J)=WX(INREC(I))
          END IF
          IF(D1.LT.0.) THEN
            USS(J)=-USS(J)
          END IF
          IF(D3.LT.0.) THEN
            WSS(J)=-WSS(J)
          END IF
 1800   CONTINUE
        if(mod(m,4).eq.0) then
C         Saving results for receivers at each 4th time level
          DO 1801 I=1,NREC
            IF(MPAR.EQ.1) THEN
              ULEV(I)=UY(INREC(I))
              WLEV(I)=WY(INREC(I))
            ELSE
              ULEV(I)=UX(INREC(I))
              WLEV(I)=WX(INREC(I))
            END IF
 1801     CONTINUE
          IF(FOUTU.NE.' ') THEN
            WRITE(LU3,'(8E10.4)') (ULEV(I),I=1,NREC)
          ENDIF
          IF(FOUTW.NE.' ') THEN
            WRITE(LU4,'(8E10.4)') (WLEV(I),I=1,NREC)
          ENDIF
        endif
      ENDIF
C
C     Writing a snapshot:
      I4=NINT((FLOAT(M-1)*DT-O4)/D4)+1
      IT=NINT((FLOAT(I4-1)*D4+O4)/DT)+1
C     I4 is the index of the snapshot nearest to the M-th time level,
C     IT is the time level nearest to the I4-th snapshot.
      IF(M.EQ.IT.AND.I4.LE.N4) THEN
        IF(FILNA1.NE.' ') THEN
          IF(MPAR.EQ.1) THEN
            CALL FDOUT1(LU1,NX,NZ,' ',UY,WORK)
          ELSE
            CALL FDOUT1(LU1,NX,NZ,' ',UX,WORK)
          END IF
        END IF
        IF(FILNA2.NE.' ') THEN
          IF(MPAR.EQ.1) THEN
            CALL FDOUT1(LU2,NX,NZ,' ',WY,WORK)
          ELSE
            CALL FDOUT1(LU2,NX,NZ,' ',WX,WORK)
          END IF
        END IF
      END IF
C
C      Artificial dumpers at right and bottom
C
c     IF(T.LT.TSWI2) GOTO 2000
      DO 1950 IA=1,NSIZE
      UY(IA)=UY(IA)*DUMP(IA)
      UX(IA)=UX(IA)*DUMP(IA)
      WY(IA)=WY(IA)*DUMP(IA)
      WX(IA)=WX(IA)*DUMP(IA)
 1950 CONTINUE
C
 2000 CONTINUE
c     call gettim(ihod2,imin2,isec2,idum)
c     write(*,*) ihod1,imin1,isec1
c     write(*,*) ihod2,imin2,isec2
C
      IF(FILNA1.NE.' ') THEN
        CLOSE(LU1)
      END IF
      IF(FILNA2.NE.' ') THEN
        CLOSE(LU2)
      END IF
C
C     Writing synthetic seismograms:
      IF(NREC.GT.0) THEN
        CALL RSEP3T('SS',OUTSS,'fd2d.gse')
        IF(OUTSS.NE.' ') THEN
          OPEN(LU1,FILE=OUTSS)
          CALL WGSE1(LU1,' ')
          DO 9000 I=1,NREC
            X=O1+D1*FLOAT(KRE(I)-1)
            Y=0.
            Z=O3+D3*FLOAT(LRE(I)-1)
            TL=0.
            CALL WGSE2(LU1,REC(I),' ',1,X,Y,Z,TL,DT,NT,USS(NT*(I-1)+1))
            CALL WGSE2(LU1,REC(I),' ',3,X,Y,Z,TL,DT,NT,WSS(NT*(I-1)+1))
 9000     CONTINUE
          CALL WGSE3(LU1)
          CLOSE(LU1)
        END IF
      END IF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      FUNCTION CDUMP(I)
C
C     Purpose:
C     Spatial variation of artificial dumpers
C
      COMMON /DMP/ NZONE
      PI=3.141593
      A=PI*FLOAT(I)/FLOAT(NZONE)
      CDUMP=0.98+0.01*(1.-COS(A))
C     CDUMP=0.94+0.03*(1.-COS(A))
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE AUXL1(KLMAX,AX,AY,BX,BY)
C
C      Purpose:
C      Transfer of the displacements from arrays UX,UY,etc.
C      into auxiliary variables UOLD,U1,U2,etc.
C      (Selecting values from vicinity of a grid point al line L=1.)
C
      DIMENSION II(25),IJ(25)
      DIMENSION AX(KLMAX),AY(KLMAX),BX(KLMAX),BY(KLMAX)
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      COMMON /DSPLC/ UOLD,WOLD,
     1           U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14,
     2           U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25,
     3           W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14,
     4           W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25
C
      UOLD=AY(ICEN)
      U3=AX(II(3))
      U5=AX(II(5))
      U8=AX(II(8))
      U9=AX(II(9))
      U13=AX(II(13))
      U14=AX(II(14))
      U15=AX(II(15))
      U18=AX(II(18))
      U19=AX(II(19))
      U23=AX(II(23))
      U25=AX(II(25))
      WOLD=BY(ICEN)
      W3=BX(II(3))
      W5=BX(II(5))
      W8=BX(II(8))
      W9=BX(II(9))
      W13=BX(II(13))
      W14=BX(II(14))
      W15=BX(II(15))
      W18=BX(II(18))
      W19=BX(II(19))
      W23=BX(II(23))
      W25=BX(II(25))
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE AUXL2(KLMAX,AX,AY,BX,BY)
C
C      Purpose:
C      Transfer of the displacements from arrays UX,UY,etc.
C      into auxiliary variables UOLD,U1,U2,etc.
C      (Selecting values from vicinity of a grid point at line L=2.)
C
      DIMENSION II(25),IJ(25)
      DIMENSION AX(KLMAX),AY(KLMAX),BX(KLMAX),BY(KLMAX)
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      COMMON /DSPLC/ UOLD,WOLD,
     1           U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14,
     2           U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25,
     3           W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14,
     4           W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25
C
      UOLD=AY(ICEN)
      U3=AX(II(3))
      U5=AX(II(5))
      U7=AX(II(7))
      U8=AX(II(8))
      U9=AX(II(9))
      U12=AX(II(12))
      U13=AX(II(13))
      U14=AX(II(14))
      U15=AX(II(15))
      U17=AX(II(17))
      U18=AX(II(18))
      U19=AX(II(19))
      U23=AX(II(23))
      U25=AX(II(25))
      WOLD=BY(ICEN)
      W3=BX(II(3))
      W5=BX(II(5))
      W7=BX(II(7))
      W8=BX(II(8))
      W9=BX(II(9))
      W12=BX(II(12))
      W13=BX(II(13))
      W14=BX(II(14))
      W15=BX(II(15))
      W17=BX(II(17))
      W18=BX(II(18))
      W19=BX(II(19))
      W23=BX(II(23))
      W25=BX(II(25))
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,AX,AY,BX,BY,
     1           T,TIMEXC,UCOM,WCOM)
C
C      Purpose:
C      Transfer of the displacements from arrays UX,UY,etc.
C      into auxiliary variables UOLD,U1,U2,etc.
C      (Selecting values from vicinity of a grid point at line L>2.)
C
      DIMENSION II(25),IJ(25)
      DIMENSION AX(KLMAX),AY(KLMAX),BX(KLMAX),BY(KLMAX)
      DIMENSION EX(4,KMAX)
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      COMMON /DSPLC/ UOLD,WOLD,
     1           U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14,
     2           U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25,
     3           W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14,
     4           W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25
C
      UOLD=AY(ICEN)
      U1=AX(II(1))
      U3=AX(II(3))
      U5=AX(II(5))
      U7=AX(II(7))
      U8=AX(II(8))
      U9=AX(II(9))
      U11=AX(II(11))
      U12=AX(II(12))
      U13=AX(II(13))
      U14=AX(II(14))
      U15=AX(II(15))
      U17=AX(II(17))
      U18=AX(II(18))
      U19=AX(II(19))
      U21=AX(II(21))
      U23=AX(II(23))
      U25=AX(II(25))
      WOLD=BY(ICEN)
      W1=BX(II(1))
      W3=BX(II(3))
      W5=BX(II(5))
      W7=BX(II(7))
      W8=BX(II(8))
      W9=BX(II(9))
      W11=BX(II(11))
      W12=BX(II(12))
      W13=BX(II(13))
      W14=BX(II(14))
      W15=BX(II(15))
      W17=BX(II(17))
      W18=BX(II(18))
      W19=BX(II(19))
      W21=BX(II(21))
      W23=BX(II(23))
      W25=BX(II(25))
      IF(L.LT.LA1.OR.L.GT.LA4) GOTO 500
      IF(T.GT.TIMEXC) GOTO 500
      IF(L.EQ.LA1) GOTO 100
      IF(L.EQ.LA2) GOTO 200
      IF(L.EQ.LA3) GOTO 300
      IF(L.EQ.LA4) GOTO 400
  100 DIF=EX(3,K-2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U5=U5+DIFU
      W5=W5+DIFW
      DIF=EX(3,K)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U15=U15+DIFU
      W15=W15+DIFW
      DIF=EX(3,K+2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U25=U25+DIFU
      W25=W25+DIFW
      GOTO 500
  200 DIF=EX(3,K-1)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U9=U9+DIFU
      W9=W9+DIFW
      DIF=EX(3,K)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U14=U14+DIFU
      W14=W14+DIFW
      DIF=EX(3,K+1)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U19=U19+DIFU
      W19=W19+DIFW
      DIF=EX(4,K-2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U5=U5+DIFU
      W5=W5+DIFW
      DIF=EX(4,K)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U15=U15+DIFU
      W15=W15+DIFW
      DIF=EX(4,K+2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U25=U25+DIFU
      W25=W25+DIFW
      GOTO 500
  300 DIF=EX(1,K-2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U1=U1-DIFU
      W1=W1-DIFW
      DIF=EX(1,K)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U11=U11-DIFU
      W11=W11-DIFW
      DIF=EX(1,K+2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U21=U21-DIFU
      W21=W21-DIFW
      DIF=EX(2,K-1)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U7=U7-DIFU
      W7=W7-DIFW
      DIF=EX(2,K)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U12=U12-DIFU
      W12=W12-DIFW
      DIF=EX(2,K+1)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U17=U17-DIFU
      W17=W17-DIFW
      GOTO 500
  400 DIF=EX(2,K-2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U1=U1-DIFU
      W1=W1-DIFW
      DIF=EX(2,K)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U11=U11-DIFU
      W11=W11-DIFW
      DIF=EX(2,K+2)
      DIFU=DIF*UCOM
      DIFW=DIF*WCOM
      U21=U21-DIFU
      W21=W21-DIFW
  500 CONTINUE
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE INDEX(K,L)
C
C      Purpose:
C      Computing 1-D index ICEN of a grid point K,L,
C      computing indices II(1)-II(25) of the grid-point neighbours,
C      computing indices isa,isb,...ish of the neighbouring small legs
c      (h/2),
C      computing indices ila,ilb,...ilh of the neighbouring large legs
C      (h).
C
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
      DIMENSION II(25),IJ(25)
c
c     ....f.......g....
c     .   f   .   g   .
c     .   f   .   g   .
c     ........J.I C E N
c     .   c   C   b   .
c     .   c   E   b   .
c     ....c...N...b....
c
c
c     .................
c     .       .       .
c     eeeeeee . hhhhhhh
c     ........J.I C E N
c     .       C       .
c     ddddddd E aaaaaaa
c     ........N........
c
c
C
      ICEN=L+(K-1)*NZ
      ICEN2=ICEN-2
      JCEN=ICEN-IFIX(FLOAT(ICEN)/FLOAT(NZ))
C
      DO 100 IPOM=1,25
      IR=IFIX((FLOAT(IPOM)-0.00001)/5.)
      NZMUL=NZ5*IR-NZ6+IPOM
      II(IPOM)=ICEN2+NZMUL
  100 CONTINUE
C
C     Indices of the grid legs, see Fig. in PARL
c         ila, ilb, etc. = indices for legs at grid lines
C
      ila=icen
      ilb=jcen
      ild=icen-nz
      ilf=jcen-1
c
c     isa, isb, etc. = indices for INTERlegs between grid lines
c
      isa=icen-(k-1)
      isd=isa-(nz-1)
      ise=isd-1
      ish=isa-1
C
      isb=jcen
      isc=jcen-(nz-1)
      isf=isc-1
      isg=isb-1
c
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PARL1(KLMAX,P2,Q2,R2,S2,P4,Q4,R4,S4,D4)
C
C      Purpose:
C      Transfer of the effective parameters from arrays AHL,AHS,etc.
C      into auxiliary variables ........
C      (selecting values from vicinity of a grid point at line L=1)
C
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      DIMENSION P2(KLMAX),Q2(KLMAX),R2(KLMAX),S2(KLMAX)
      DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX),D4(KLMAX)
      DIMENSION II(25),IJ(25)
C
c #####################################################
c Density for the central point
      dnst=d4(icen)
c #####################################################
c
c The parameters for the non-mixed  derivatives Dx(p Dx)
c        horizontal legs de and ah
      pdes=p4(ild)
      pahs=p4(ila)
c
c The parameters for the non-mixed  derivatives Dz(r Dz)
c        vertical   legs fg and bc
      rfgs=0.
      rbcs=r4(ilb)
c
c The parameters for the mixed  derivatives Dx(s Dz)
c       horizontal (inter)legs h,a,e,d
      shs=0.
      sas=s2(isa)
      ses=0.
      sds=s2(isd)
c
c equating legs h-a and d-e
c     shs=s4(ila)
c     sas=shs
c     sds=s4(ild)
c     ses=sds
c
c The parameters for the mixed  derivatives Dz(q Dx)
c       vertical  (inter)legs f,c,g,b
      qfs=0.
      qcs=q2(isc)
      qgs=0.
      qbs=q2(isb)
c
c equating legs f-g and b-c
c     qfs=0.
c     qgs=0.
c     qbs=q4(ilb)
c     qcs=qbs
c
c ######################################################
c The parameters for the outer formula (4th order)
c           not yet OK
c ######################################################
c
c
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PARL2(KLMAX,P2,Q2,R2,S2,P4,Q4,R4,S4,D4)
C
C      Purpose:
C      Transfer of the effective parameters from arrays AHL,AHS,etc.
C      into auxiliary variables ........
C      (selecting values from vicinity of a grid point at line L=2)
C
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      DIMENSION P2(KLMAX),Q2(KLMAX),R2(KLMAX),S2(KLMAX)
      DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX),D4(KLMAX)
      DIMENSION II(25),IJ(25)
C
c #####################################################
c Density for the central point
      dnst=d4(icen)
c ##################################################
c
c The parameters for the non-mixed  derivatives Dx(p Dx)
c        horizontal legs de and ah
      pdes=p4(ild)
      pahs=p4(ila)
c
c The parameters for the non-mixed  derivatives Dz(r Dz)
c        vertical   legs fg and bc
      rfgs=r4(ilf)
      rbcs=r4(ilb)
c
c The parameters for the mixed  derivatives Dx(s Dz)
c       horizontal (inter)legs h,a,e,d
c     shs=s2(ish)
c     sas=s2(isa)
c     ses=s2(ise)
c     sds=s2(isd)
c equating legs h-a and d-e
      shs=s4(ila)
c     sas=shs
      sds=s4(ild)
c     ses=sds
c
c The parameters for the mixed  derivatives Dz(q Dx)
c       vertical  (inter)legs f,c,g,b
c     qfs=q2(isf)
c     qcs=q2(isc)
c     qgs=q2(isg)
c     qbs=q2(isb)
c equating legs f-g and b-c
      qfs=q4(ilf)
c     qgs=qfs
      qbs=q4(ilb)
c     qcs=qbs
c
c ######################################################
c The parameters for the outer (large) square - 4th order
c   not yet OK
c ######################################################
c
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PARL(KLMAX,P2,Q2,R2,S2,P4,Q4,R4,S4,D4)
C
C      Purpose:
C      Transfer of the effective parameters from arrays AH,BH,etc.
C      into auxiliary variables ........
C      (selecting values from vicinity of a grid point at line L>2)
C
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      DIMENSION P2(KLMAX),Q2(KLMAX),R2(KLMAX),S2(KLMAX)
      DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX),D4(KLMAX)
      DIMENSION II(25),IJ(25)
C
c #####################################################
c Density for the central point, formally written
c in the array indexed as 'large' horizontal parameters
c and/or the array of the displacements
      dnst=d4(icen)
c #####################################################
c
c The parameters for the inner formula (2nd order)
c   - denotation of the used arrays: p4,q4,r4,s4
c        - they correspond to the legs
c          coinciding with grid lines
c        - numbering of the legs by indices
c          ila, ilb, ild, ilf
c    - denotation of the used arrays: p2,q2,r2,s2
c        - they correspond to the legs
c          half a way
c          between the grid lines
c        - numbering of the legs by indices
c          isa, isb, ... ish
c ##################################################
c
c The parameters for the non-mixed  derivatives Dx(p Dx)
c        horizontal legs de and ah
      pdes=p4(ild)
      pahs=p4(ila)
c
c The parameters for the non-mixed  derivatives Dz(r Dz)
c        vertical   legs fg and bc
      rfgs=r4(ilf)
      rbcs=r4(ilb)
c
c The parameters for the mixed  derivatives Dx(s Dz)
c       horizontal (inter)legs h,a,e,d
c     shs=s2(ish)
c     sas=s2(isa)
c     ses=s2(ise)
c     sds=s2(isd)
c equating legs h-a and d-e
      shs=s4(ila)
c     sas=shs
      sds=s4(ild)
c     ses=sds
c
c The parameters for the mixed  derivatives Dz(q Dx)
c       vertical  (inter)legs f,c,g,b
c     qfs=q2(isf)
c     qcs=q2(isc)
c     qgs=q2(isg)
c     qbs=q2(isb)
c equating legs f-g and b-c
      qfs=q4(ilf)
c     qgs=qfs
      qbs=q4(ilb)
c     qcs=qbs
c ######################################################
c The parameters for the outer formula (4th order)
c               for EQUATED legs only !!!!
c ######################################################
c
c The parameters for the non-mixed  derivatives
c       double legs de and ah
c     pdel=1./(1./p4(ild)+1./p4(ild-nz))
c     pahl=1./(1./p4(ila)+1./p4(ila+nz))
c
c The parameters for the non-mixed  derivatives
c       double legs fg and bc
c     rfgl=1./(1./r4(ilf)+1./r4(ilf-1))
c     rbcl=1./(1./r4(ilb)+1./r4(ilb+1))
c
c The parameters for the mixed  derivatives
c         legs h,a,e,d
c     shl=1./(1./s4(ila)+1./s4(ila+nz))
cc    sal=shl
c     sdl=1./(1./s4(ild)+1./s4(ild-nz))
cc    sel=sdl
c
c The parameters for the mixed  derivatives
c         legs f,c,g,b
c     qfl=1./(1./q4(ilf)+1./q4(ilf-1))
cc    qgl=qfl
c     qbl=1./(1./q4(ilb)+1./q4(ilb+1))
cc    qcl=qbl
c
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE PARHOM(KLMAX,P4,Q4,R4,S4,D4)
C
C      Purpose:
C      Analogous to PARL but for a grid point with a completely
C      homogeneous vicinity
C
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /INDI/ ICEN,II,IJ,
     1              isa,isb,isc,isd,ise,isf,isg,ish,
     2              ila,ilb,ilc,ild,ile,ilf,ilg,ilh
      DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX)
      DIMENSION D4(KLMAX)
      DIMENSION II(25),IJ(25)
C
      dnst=d4(icen)
c
      phom=p4(ila)
      qhom=q4(ilb)
      rhom=r4(ilb)
      shom=s4(ila)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE DIFL1 (AOLD,A3,A8,A13,A14,A15,
     1                  A18,A23,
     2                  B3,B5,B8,B9,B13,B14,B15,
     3                  B18,B19,B23,B25,
     4                  DISNEW)
C
C      Purpose:
C      The finite-difference scheme for a grid point at the free surface
C      at line L=1
c
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
c inner part (step h in displ.): non mixed + 0.25*(mixed)
c
c   Attention: Zero-valued parameters have been defined in PARL1.
c              They appear here only formally.
c
c
      pintp=pahs*(a18-a13)-pdes*(a13-a8)
      pintr=rbcs*(a14-a13)+rfgs*(0.-a13)
      pints=0.25*( shs*(b13+b18-0. -0. )+sas*(-b13-b18+b19+b14)-
     3             ses*(b8+b13-0.-0. )  -sds*(-b8-b13+b14+b9) )
      pintq=0.25*(-qgs*(0. +b18-0. -b13)-qfs*(-b8-0.+0. +b13)+
     5             qbs*(b18+b19-b13-b14)+qcs*(-b9-b8+b13+b14) )
c
c outer part (step 2h in displ.): non mixed + 0.25*(mixed)
c         (FOR 4th order ONLY)
c
c     poutp=pahl*(a23-a13)-pdel*(a13-a3)
c     poutr=rbcl*(a15-a13)+rfgl*(0.-a13)
c     pouts=0.25*( shl*(b13+b23-0. -0. )+sal*(-b13-b23+b25+b15)-
c    3             sel*(b3+b13-0.-0. )  -sdl*(-b3-b13+b15+b5) )
c     poutq=0.25*(-qgl*(0. +b23-0. -b13)-qfl*(-b3-0.+0. +b13)+
c    5             qbl*(b23+b25-b13-b15)+qcl*(-b5-b3+b13+b15) )
c
c  4th order scheme
c     sleft=16.*pint-pout
c
c  2nd order    (POUT shouldn't be computed in this case)
c     sleft=12.*pint
c
      sleftp=12.*pintp
      sleftr=12.*pintr
      slefts=12.*pints
      sleftq=12.*pintq
c
      sleft=sleftp+sleftr+slefts+sleftq
c
c     DISNEW=2.*CNST*sleft+2.*A13-AOLD
      DISNEW=0.
      IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE DIFL2 (AOLD,A3,A8,A12,A13,A14,A15,
     1                  A18,A23,
     2                  B3,B5,B7,B8,B9,B12,B13,B14,B15,B17,
     3                  B18,B19,B23,B25,
     4                  DISNEW)
C
C      Purpose:
C      The finite-difference scheme for a grid point in a heterogeneous
C      medium at line L=2
C
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
c inner part (step h): non mixed + 0.25*(mixed)
      pintp=pahs*(a18-a13)-pdes*(a13-a8)
      pintr=rbcs*(a14-a13)+rfgs*(a12-a13)
c     pints=0.25*( shs*(b13+b18-b12-b17)+sas*(-b13-b18+b19+b14)-
c    3             ses*(b8+b13-b7-b12)  -sds*(-b8-b13+b14+b9) )
c     pintq=0.25*(-qgs*(b17+b18-b12-b13)-qfs*(-b8-b7+b12+b13)+
c    5             qbs*(b18+b19-b13-b14)+qcs*(-b9-b8+b13+b14) )
c  mixed derivative with equated legs
      pints=0.25*( shs*(-b12-b17+b19+b14)-
     3             sds*(-b7-b12+b14+b9) )
      pintq=0.25*(-qfs*(b17+b18-b8-b7)+
     5             qbs*(b18+b19-b9-b8) )
c outer part (step 2h): non mixed + 0.25*(mixed)
c         (FOR 4th order ONLY)
c
c  !!!!!!!!!! Attention: Here rfgl, qgl, qfl are problematic !!!!!!!!!!!!
c  !!!!!!!!!! because of the coincidence of legs F and G with the surface
c  !!!!!!!!!! Possible actions: redefine them in PARL2, !!!!!!!!!!!!!!!!!
c  !!!!!!!!!! or use only the 2nd order formula here !!!!!!!!!!!!!!!!!!!!
c
c     poutp=pahl*(a23-a13)-pdel*(a13-a3)
c     poutr=rbcl*(a15-a13)+rfgl*(0.-a13)
c     pouts=0.25*( shl*(b13+b23-0. -0. )+sal*(-b13-b23+b25+b15)-
c    3             sel*(b3+b13-0.-0. )  -sdl*(-b3-b13+b15+b5) )
c     poutq=0.25*(-qgl*(0. +b23-0. -b13)-qfl*(-b3-0.+0. +b13)+
c    5             qbl*(b23+b25-b13-b15)+qcl*(-b5-b3+b13+b15) )
c
c  4th order scheme
c     sleft=16.*pint-pout
c
c  2nd order    (POUT shouldn't be computed in this case)
c     sleft=12.*pint
c
c
      sleftp=12.*pintp
      sleftr=12.*pintr
      slefts=12.*pints
      sleftq=12.*pintq
c
      sleft=sleftp+sleftr+slefts+sleftq
c
c     DISNEW=CNST*sleft+2.*A13-AOLD
      DISNEW=0.
      IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE DIFL  (AOLD,A3,A8,A11,A12,A13,A14,A15,
     1                  A18,A23,
     2                  B1,B3,B5,B7,B8,B9,B11,B12,B13,B14,B15,B17,
     3                  B18,B19,B21,B23,B25,
     4                  DISNEW,NOUT)
C
C      Purpose:
C      The finite-difference scheme for a grid point in a heterogeneous
C      medium at line L>2
c
c   Displacement numbering when updating the central point (#13):
c   1  6 11 16 21
c   2  7 12 17 22
c   3  8 13 18 23
c   4  9 14 19 24           A1, A2, etc. for non mixed derivatives,
c   5 10 15 20 25           B1, B2, etc. for mixed derivatives.
c
c
c
c    General denotation    p     q      r     s     A    B
c    means
c    when updating u:     ah    bv     bv    ch     u    w
c    when updating w:     bh    cv     av    bh     w    u
c
c    here 1st char.: a=lambda+2mu, b=mu, c=lambda
c         2nd char.: v=vertical legs, h=horizontal legs
c
c
c
c   Denotation of the effective parameters for non-mixed derivatives:
c
c   |---------|---------|    1st char.: the parameter (p or r)
c   |         r         |    2nd and 3rd char.: double legs (f and g,
c   |         f         |                                     a and h,
c   |         g         |                                     etc.)
c   ---pde----|---pah----
c   |         r         |        (p Ax)x, (r Az)z
c   |         b         |
c   |         c         |
c   |---------|---------|
c
c
c
c
c   Denotation of the effective parameters for mixed derivatives:
c
c   |---------|---------|    1st char.: the parameter (q or s)
c   |    qf   |    qg   |    2nd char.: the leg (a,b,...g,h)
c   |    qf   |    qg   |
c   |    qf   |    qg   |
c   |---------|---------|
c   |    qc   |    qb   |        (q Bx)z
c   |    qc   |    qb   |
c   |    qc   |    qb   |
c   |---------|---------|
c
c   Denotation of the effective parameters for mixed derivatives:
c
c   |---------|---------|    1st char.: the parameter (q or s)
c   |         |         |    2nd char.: the leg (a,b,...g,h)
c   | se se   |  sh sh  |
c   |         |         |
c   |---------|---------|
c   |         |         |        (s Bz)x
c   | sd sd   |  sa sa  |
c   |         |         |
c   |---------|---------|
c
c
c   When using the parameters in the 'inner' formulas (step h)
c   they are denoted shs, sas, etc.
c
c   When using the parameters in the 'outer' formulas (step 2h)
c   they are denoted shl, sal, etc.
c
c
c
c       configuration of the grid points (for the displacements)
c                    and the grid legs   (for the parameters):
c
c
c   1......6......11......16......21
c   .      f       .       g       .
c   .      f       .       g       .
c   .      f       .       g       .
c   2 eeee 7 eeee 12 hhhh 17 hhhh 22       Grid legs a,b,...g,h
c   .      f   f   .   g   g       .       for  the 'small'
c   .      f ee.eee.hh..hh g       .       and 'large' squares,
c   .      f   f   .   g   g       .       used in the inner and
c   3......8......13......18......23       outer formulas, resp.
c   .      d   c   .   b   b       .
c   .      d dd.ddd.aa..aa b      .
c   .      d   c   .   b   b       .
c   4 dddd 9 dddd 14 aaaa 19 aaaa 24
c   .      d       .       b       .
c   .      d       .       b       .
c   .      d       .       b       .
c   5.....10......15......20......25
c
c
c
c    CNST=dt*dt/dx*dx/12
c
c
c
c
C
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
c inner part (step h): non mixed + 0.25*(mixed)
      pintp=pahs*(a18-a13)-pdes*(a13-a8)
      pintr=rbcs*(a14-a13)+rfgs*(a12-a13)
c     pints=0.25*( shs*(b13+b18-b12-b17)+sas*(-b13-b18+b19+b14)-
c    3             ses*(b8+b13-b7-b12)  -sds*(-b8-b13+b14+b9) )
c     pintq=0.25*(-qgs*(b17+b18-b12-b13)-qfs*(-b8-b7+b12+b13)+
c    5             qbs*(b18+b19-b13-b14)+qcs*(-b9-b8+b13+b14) )
c  mixed derivative with equated legs
      pints=0.25*( shs*(-b12-b17+b19+b14)-
     3             sds*(-b7-b12+b14+b9) )
      pintq=0.25*(-qfs*(b17+b18-b8-b7)+
     5             qbs*(b18+b19-b9-b8) )
c outer part (step 2h): non mixed + 0.25*(mixed)
c         (FOR 4th order ONLY)
c
c     poutp=pahl*(a23-a13)-pdel*(a13-a3)
c     poutr=rbcl*(a15-a13)+rfgl*(a11-a13)
cc    pouts=0.25*( shl*(b13+b23-b11-b21)+sal*(-b13-b23+b25+b15)-
cc   3             sel*(b3+b13-b1-b11)  -sdl*(-b3-b13+b15+b5) )
cc    poutq=0.25*(-qgl*(b21+b23-b11-b13)-qfl*(-b3-b1+b11+b13)+
cc   5             qbl*(b23+b25-b13-b15)+qcl*(-b5-b3+b13+b15) )
cc mixed derivative with equated legs
c     pouts=0.25*( shl*(-b12-b17+b19+b14)-
c    3             sdl*(-b7-b12+b14+b9) )
c     poutq=0.25*(-qfl*(b17+b18-b8-b7)+
c    5             qbl*(b18+b19-b9-b8) )
cc
cc 4th order scheme
cc    sleft=16.*pint-pout
cc
cc 2nd order    (POUT shouldn't be computed in this case)
cc    sleft=12.*pint
cc
c
      sleftp=12.*pintp
      sleftr=12.*pintr
      slefts=12.*pints
      sleftq=12.*pintq
c
      sleft=sleftp+sleftr+slefts+sleftq
c
c     DISNEW=CNST*sleft+2.*A13-AOLD
      DISNEW=0.
      IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE DIFHOM(AOLD,A3,A8,A11,A12,A13,A14,A15,A18,A23,
     2                  B1,B5,B7,B9,B17,B19,B21,B25,
     4                  DISNEW)
C
C      Purpose:
C      The finite-difference scheme for a grid point in a homogeneous
C      medium
C
      COMMON /PPAR/ pahs,pdes,rbcs,rfgs,
     1              shs,sas,ses,sds,qgs,qfs,qbs,qcs,
     2              pahl,pdel,rbcl,rfgl,
     3              shl,sal,sel,sdl,qgl,qfl,qbl,qcl,
     4              phom,qhom,rhom,shom,dnst
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
c  inner part
      pint=phom*(A18-2.*A13+A8)+rhom*(A14-2.*A13+A12)
     1    +0.25*(shom+qhom)*(b19-b17+b7-b9)
c  outer part  (for 4th order only)
c     pout=phom*(A23-2.*A13+A3)+rhom*(A15-2.*A13+A11)
c    1    +0.25*(shom+qhom)*(b25-b21+b1-b5)
c
c 4th order
c     sleft=16.*pint-pout
c 2nd order (POUT shouldn't be computed in this case)
      sleft=12.*pint
c
c     DISNEW=CNST*sleft+2.*A13-AOLD
      DISNEW=0.
      IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE FLERIB(KLMAX,LMAX,LIM1,LIM2,IPT1,IPT2,IPT3,
     1                  ULINA,ULINB,WLINA,WLINB,
     2                  UX,UY,WX,WY,AH,BH,DENS)
C
C      Purpose:
C      Left and right nonreflecting boundaries
C
      DIMENSION AH(KLMAX),BH(KLMAX),DENS(KLMAX)
      DIMENSION UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX)
      DIMENSION ULINA(LMAX),ULINB(LMAX),WLINA(LMAX),WLINB(LMAX)
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
      DO 500 L=LIM1,LIM2
C
c     ALFA=SQRT(AH(IPT3+L))
      ALFA=0.
      IF(DENS(IPT3+L).NE.0.) ALFA=SQRT(AH(IPT3+L)/DENS(IPT3+L))
      C2=CN*ALFA
      IF(L.EQ.1) C2=C2*1.414
      C3=1./(1.+C2)
C
c     BETA=SQRT(BH(IPT3+L))
      BETA=0.
      IF(DENS(IPT3+L).NE.0.) BETA=SQRT(BH(IPT3+L)/DENS(IPT3+L))
      C4=CN*BETA
      IF(L.EQ.1) C4=C4*1.414
      C5=1./(1.+C4)
C
cc    C6=(BETA-ALFA)/ALFA
cc    C7=(BETA-ALFA)/BETA
C
C
      IF(MPAR.EQ.1) THEN
       PAR1=UY(IPT1+L)
       PAR2=UX(IPT2+L)
       PAR3=UX(IPT1+L)
       PAR4=UX(IPT2+L+1)
       RAR1=WY(IPT1+L)
       RAR2=WX(IPT2+L)
       RAR3=WX(IPT1+L)
       RAR4=WX(IPT2+L+1)
      ELSE
       PAR1=UX(IPT1+L)
       PAR2=UY(IPT2+L)
       PAR3=UY(IPT1+L)
       PAR4=UY(IPT2+L+1)
       RAR1=WX(IPT1+L)
       RAR2=WY(IPT2+L)
       RAR3=WY(IPT1+L)
       RAR4=WY(IPT2+L+1)
      ENDIF
C
C     Emerman, Stephen
      UAUX=C2*(PAR1-ULINA(L)+ULINB(L))-
     1        (PAR1-2.*(PAR2+PAR3)+ULINB(L)+ULINA(L))
      WAUX=C4*(RAR1-WLINA(L)+WLINB(L))-
     1        (RAR1-2.*(RAR2+RAR3)+WLINB(L)+WLINA(L))
C
C     Stacey P3
C     C3=1.
C     C5=1.
C     UAUX=-1.*C2*(PAR2-PAR3)+C6*C2*(RAR4-RAR2)+PAR2
C     WAUX=-1.*C4*(RAR2-RAR3)+C7*C4*(PAR4-PAR2)+RAR2
C
C     Clayton, Engquist
C     C3=1.
C     C5=1.
C     UAUX=-1.*C2*(PAR2-PAR3)+PAR2
C     WAUX=-1.*C4*(RAR2-RAR3)+RAR2
C
C     Reynolds
C  Attention: CHANGE ALSO CALL FLERIB AND IZ4 IN SAVER !!!!!!!!!!!!
C     C3=1.
C     C5=1.
C     UAUX=-1.*C2*(PAR2-PAR3-ULINB(L)+ULINA(L))+(PAR2+PAR3-ULINB(L))
C     WAUX=-1.*C4*(RAR2-RAR3-WLINB(L)+WLINA(L))+(RAR2+RAR3-WLINB(L))
C
      IF(MPAR.EQ.1) THEN
       UY(IPT2+L)=C3*UAUX
       WY(IPT2+L)=C5*WAUX
      ELSE
       UX(IPT2+L)=C3*UAUX
       WX(IPT2+L)=C5*WAUX
      ENDIF
  500 CONTINUE
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE FBOB(KLMAX,KMAX,KIM1,KIM2,IPT1,IPT2,
     1                ULINA,ULINB,WLINA,WLINB,
     2                UX,UY,WX,WY,AV,BV,DENS)
C
C      Purpose:
C      Bottom nonreflecting boundaries
C
      DIMENSION AV(KLMAX),BV(KLMAX),DENS(KLMAX)
      DIMENSION UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX)
      DIMENSION ULINA(KMAX),ULINB(KMAX),WLINA(KMAX),WLINB(KMAX)
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
      DO 500 K=KIM1,KIM2
      IPT=IPT1+(K-1)*NZ
      IPTV=K*(NZ-1)-IPT2
      NZ99=K*NZ-IPT2
C
c     C2=CN*SQRT(BV(IPTV))
      C2=0.
      IF(DENS(NZ99).NE.0.) C2=CN*SQRT(BV(IPTV)/DENS(NZ99))
      C3=1./(1.+C2)
c     C4=CN*SQRT(AV(IPTV))
      C4=0.
      IF(DENS(NZ99).NE.0.) C4=CN*SQRT(AV(IPTV)/DENS(NZ99))
      C5=1./(1.+C4)
C
      IF(MPAR.EQ.1) THEN
       PAR1=UY(IPT)
       PAR2=UX(IPT+1)
       PAR3=UX(IPT)
       RAR1=WY(IPT)
       RAR2=WX(IPT+1)
       RAR3=WX(IPT)
      ELSE
       PAR1=UX(IPT)
       PAR2=UY(IPT+1)
       PAR3=UY(IPT)
       RAR1=WX(IPT)
       RAR2=WY(IPT+1)
       RAR3=WY(IPT)
      ENDIF
      UAUX=C2*(PAR1-ULINA(K)+ULINB(K))-
     1        (PAR1-2.*(PAR2+PAR3)+ULINB(K)+ULINA(K))
      WAUX=C4*(RAR1-WLINA(K)+WLINB(K))-
     1        (RAR1-2.*(RAR2+RAR3)+WLINB(K)+WLINA(K))
      IF(MPAR.EQ.1) THEN
       UY(IPT+1)=C3*UAUX
       WY(IPT+1)=C5*WAUX
      ELSE
       UX(IPT+1)=C3*UAUX
       WX(IPT+1)=C5*WAUX
      ENDIF
  500 CONTINUE
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE FCORNB(KLMAX,NIND,NDIF1,NDIF2,NLFT,
     1                  UX,UY,WX,WY,AH,BH,DENS)
C
C        Purpose:
C      Left and right bottom corners of the nonreflecting boundaries
C           The 'corner' means 3 points:   x    , or      x
C                                          x x          x x
C
      DIMENSION AH(KLMAX),BH(KLMAX),DENS(KLMAX)
      DIMENSION UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX)
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
c     ALFA=SQRT(AH(NIND-NDIF1))
c     BETA=SQRT(BH(NIND-NDIF1))
      ALFA=0.
      NZ99=NIND-NDIF1
      IF(DENS(NZ99).NE.0.) ALFA=SQRT(AH(NIND-NDIF1)/DENS(NZ99))
      BETA=0.
      IF(DENS(NZ99).NE.0.) BETA=SQRT(BH(NIND-NDIF1)/DENS(NZ99))
      C4=ALFA+BETA
      IF(NLFT.EQ.1) THEN
       C5=BETA-ALFA
      ELSE
       C5=ALFA-BETA
      ENDIF
      C6=CN/1.4142135/2.
C**************************
      DO 500 I=1,3
      IF(I.EQ.1) IND=NIND-1
      IF(I.EQ.2) IND=NIND-NDIF2
      IF(I.EQ.3) IND=NIND
C**************************
C     IND=NIND
C**************************
      IF(MPAR.EQ.1) THEN
       UAUX=C4*(UX(IND-1)-2.*UX(IND)+UX(IND-NDIF2))+
     1      C5*(WX(IND-1)-2.*WX(IND)+WX(IND-NDIF2))
       WAUX=C4*(WX(IND-1)-2.*WX(IND)+WX(IND-NDIF2))+
     1      C5*(UX(IND-1)-2.*UX(IND)+UX(IND-NDIF2))
       UY(IND)=C6*UAUX+UX(IND)
       WY(IND)=C6*WAUX+WX(IND)
      ELSE
       UAUX=C4*(UY(IND-1)-2.*UY(IND)+UY(IND-NDIF2))+
     1      C5*(WY(IND-1)-2.*WY(IND)+WY(IND-NDIF2))
       WAUX=C4*(WY(IND-1)-2.*WY(IND)+WY(IND-NDIF2))+
     1      C5*(UY(IND-1)-2.*UY(IND)+UY(IND-NDIF2))
       UX(IND)=C6*UAUX+UY(IND)
       WX(IND)=C6*WAUX+WY(IND)
      ENDIF
C**********************
  500 CONTINUE
C***********************
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE SAVER(KMAX,LMAX,KLMAX,A,B,
     1                 UK1,UK2,UK3,UNX2,UNX1,UNX,UNZ2,UNZ1,UNZ,
     2                 WK1,WK2,WK3,WNX2,WNX1,WNX,WNZ2,WNZ1,WNZ)
C
C      Purpose:
C      Saving 'very old' values (from time level M-1) at boundary
C      lines for subsequent computations at nonreflecting boundaries
C
      DIMENSION UK1(LMAX),UK2(LMAX),UK3(LMAX)
      DIMENSION UNX2(LMAX),UNX1(LMAX),UNX(LMAX)
      DIMENSION UNZ2(KMAX),UNZ1(KMAX),UNZ(KMAX)
      DIMENSION WK1(LMAX),WK2(LMAX),WK3(LMAX)
      DIMENSION WNX2(LMAX),WNX1(LMAX),WNX(LMAX)
      DIMENSION WNZ2(KMAX),WNZ1(KMAX),WNZ(KMAX)
      DIMENSION A(KLMAX),B(KLMAX)
      COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST
C
      DO 100 I=1,NZ
      UK1(I)=A(I)
      WK1(I)=B(I)
      UK2(I)=A(NZ+I)
      WK2(I)=B(NZ+I)
      IZ1=2*NZ+I
      UK3(I)=A(IZ1)
      WK3(I)=B(IZ1)
      IZ2=NSIZE-3*NZ+I
      UNX2(I)=A(IZ2)
      WNX2(I)=B(IZ2)
      IZ3=IZ2+NZ
      UNX1(I)=A(IZ3)
      WNX1(I)=B(IZ3)
C For Emerman-Stephen:
      IZ4=IZ3+NZ
C For Reynolds:
C     IZ4=NSIZE-4*NZ+I
      UNX(I)=A(IZ4)
      WNX(I)=B(IZ4)
  100 CONTINUE
      DO 200 I=1,NX
      IZ5=NZ-2+(I-1)*NZ
      UNZ2(I)=A(IZ5)
      WNZ2(I)=B(IZ5)
      IZ6=IZ5+1
      UNZ1(I)=A(IZ6)
      WNZ1(I)=B(IZ6)
      IZ7=IZ6+1
      UNZ(I)=A(IZ7)
      WNZ(I)=B(IZ7)
  200 CONTINUE
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      FUNCTION TIMFUN(T,TAU,TP)
C
C      Purpose:
C      The excitation signal, which is the half-th integral of a seismic
C      force or the half-th derivative of a seismic moment
C
      TIMFUN=0.
      IF(T.LT.TAU.OR.T.GT.TAU+2.*TP) RETURN
      PI=3.141593
C
C     Type of the signal
      CALL RSEP3I('KSIG',KSIG,3)
C     Reference frequency SIGF=1/(2*TP)
C     Reference time      SIGT=TP
C     Reference amplitude
      CALL RSEP3R('SIGA',SIGA,1.)
C     Reference width in reference half-periods:
      CALL RSEP3R('SIGW',SIGW,1.)
C     Phase
      CALL RSEP3R('SIGPH',SIGPH,1.)
C
      GO TO (1,1,10,1,30,1,1,1) KSIG+2
    1 CONTINUE
C     FD2D-06
      CALL ERROR('FD2D-06: Only signals KSIG=1,3 allowed')
C
C    Gabor:
C    [Hron: frequency SIGF=1, width SIGW=4, phase SIGPH=-pi/2]
   10 CONTINUE
      FNUL=1.
      ARA=2.*PI*FNUL*(T-TAU-TP)
      ARAG=ARA/SIGW
      ARAG=-1.*ARAG*ARAG
      TIMFUN=SIGA*SIN(ARA)*EXP(ARAG)
      RETURN
C
C    [Jih: frequency SIGF=0.45, width SIGW=0.05, phase SIGPH=0]
   11 CONTINUE
      FNUL=0.45
      GAMA=0.05
      ARA=2.*PI*FNUL*(T-TAU-TP)
      ARAG=ARA/GAMA
      ARAG=-1.*ARAG*ARAG
      TIMFUN=COS(ARA)*EXP(ARAG)
      RETURN
C
C    [Bard, model B:
C     frequency SIGF=2.416666, width SIGW=7.5, phase SIGPH=-pi]
   12 CONTINUE
      FNUL=2.416666
      GAMA=7.5
      ARA=2.*PI*FNUL*(T-TAU-TP)
      ARAG=ARA/GAMA
      ARAG=-1.*ARAG*ARAG
      TIMFUN=-1.*COS(ARA)*EXP(ARAG)
      RETURN
C
C    Hermite-Gaussian (Ricker) [Gauss: order SIGW=1]:
   21 CONTINUE
      GAMA=40.
      ARB=T-TAU-TP
      ARD=-2.*GAMA*ARB
      ARE=-1.*GAMA*ARB*ARB
      TIMFUN=ARD*EXP(ARE)
      RETURN
C
C    Hermite-Gaussian (Ricker) [Ricker: order SIGW=2]:
   22 CONTINUE
      PI2=PI*PI
      A=PI2*(T-TAU-TP)*(T-TAU-TP)/TP/TP
      XMAX=TP*(PI2-0.5)/SQRT(PI)
      TIMFUN=(PI2-A)*EXP(-1.*A)/XMAX
      RETURN
C
C    Hermite-Gaussian (Ricker) [Bard, models E,F: SIGW=2]
   23 CONTINUE
c     FNUL=0.862070  for test E
c     FNUL=0.714286  for test F
      ARG=PI*FNUL*(T-TAU-TP)
      ARG=ARG*ARG
      TIMFUN=-0.707107*(ARG-0.5)*EXP(-1.*ARG)
      RETURN
C
C    Hermite-Gaussian (Ricker) [Kawase: SIGW=2]
   24 CONTINUE
      FC=0.25
      ARG=PI*PI*FC*FC*(T-TAU-TP)*(T-TAU-TP)
      TIMFUN=(2.*ARG-1.)*EXP(-1.*ARG)
      RETURN
C
C    Kuepper (Mueller) (duration=2TP):
   30 CONTINUE
C     Reference amplitude is the maximum amplitude
      AMPL=SIGA*4./SQRT(27.)
      DPI=PI+PI
      ARGMN1=DPI*(T-TAU)/(2.*TP)
      ARGMN2=2.*ARGMN1
      TIMFUN=AMPL*(SIN(ARGMN1)-0.5*SIN(ARGMN2))
      RETURN
C
c    Priolo  (max. frequency 22 Hz)
      fmax=22.
      eta=0.5
      eps=1.
      s1=-fmax*fmax*(t-tp)*(t-tp)*eta
      s2=eps*pi*fmax*(t-tp)
cc  sign. changed in contrast to Priolo definition due to
cc  opposite sign. in my source (-exp  changed to +exp)
      TIMFUN=+exp(s1)*fmax*
     1       (2.*fmax*eta*cos(s2)*(t-tp)+eps*pi*sin(s2))
      RETURN
c
C    Aboudi
      TIMFUN=0.
      X=T-TAU
      IF(X.LT.0..OR.X.GT.2.*TP) RETURN
      A=ABS(X-TP)
      B=2./TP/TP
      IF(A.GE.0..AND.A.LE.TP/2.) GOTO 100
      IF(A.GE.TP/2..AND.A.LE.TP) GOTO 200
  100 TIMFUN=1.-B*A*A
      RETURN
  200 C=A-TP/2.
      TIMFUN=1.-B*(A*A-2.*C*C)
      RETURN
C
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'fdaux.for'
C     fdaux.for
      INCLUDE 'gmtra.for'
C     gmtra.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'gse.for'
C     gse.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C