C 
C Program FDMOD to generate finite-difference effective parameters
C (harmonic averages of elastic parameters) in 2-D models specified in
C terms of homogeneous convex polygons
C
C Version: 5.30
C Date: 1999, April 2
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 Data input and output modified by Ludek Klimes
C
C.......................................................................
C
C Contents:
C     Description of the polygonal model
C     Description of data files
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 the polygonal model:
C
C The medium is composed of homogeneous polygonal blocks. The
C blocks do not overlap each other and cover the computational region.
C Convex blocks only are allowed (internal angles.le.180 degrees).
C For this  reason even some physically homogeneous parts
C of the model must be subdivided in several (formal) blocks.
C Numbering of the blocks is arbitrary.
C Each block is characterized by the physical coordinates of its
C nodes (in meters), S-wave velocity (m/s), P/S velocity ratio,
C and density (kg/m**3). The nodes are numbered counterclockwise.
C Also specified for each block segment is the serial number of the
C neighbouring block.
C See the accompanying example, or ask the author of the program
C for more details.
C
C Description and remarks on FD method
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 the filename of the polygonal model,
C             dimensions of the grid of N1*N2*N3 gridpoints, the form
C             and names of the files with effective elastic parameters,
C             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 Data specifying the model by means of the MODEL package:
C     FDMOD='string'... Name of the input data file describing a 2-D
C             model composed of homogeneous convex polygons.
C             The data file contains the name of the model, the number
C             of blocks, and for each block:
C             number of edges of the block, shear-wave velocity
C             (e.g., in m/s) inside the block, P/S velocity ratio inside
C             the block, density (e.g., in kg/m**3) inside the block,
C             x and z coordinates (e.g., in m) of the initial point of
C             segment, sequential number of the neighbouring block
C             (outside put 0).
C             Parameter FDMOD is required by this program (but no other
C             one).
C             Example of data set FDMOD
C             Default: FDMOD=' '
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.  Default: N1=1
C     N3=positive integer... Number of gridpoints of the basic grid
C             along the X3 axis.  Default: N3=1
C     O1=real... X1 coordinate of the origin of the grid. 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.  Default: D1=1.
C     D3=real... Grid spacing along the X3 axis.  Default: D3=1.
C             The grid intervals may also be negative.
C And the following parameters specific to effective-parameter FD codes:
C     FORM='string'... Form of the output files:
C             FORM='formatted': Formatted ASCII output files,
C             FORM='unformatted': Unformatted output files.
C     A11='string', B11='string', C11='string'... Strings in apostrophes
C             containing the names of the output files with (N1-1)*N2*N3
C             values of the elastic parameters averaged in the direction
C             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 output files with N1*N2*(N3-1)
C             values of the elastic parameters averaged in the direction
C             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 output files with
C             (N1-1)*N2*(N3-1) values of the elastic parameters averaged
C             in the direction the X1 axis, along the gridlines shifted
C             half a grid interval in the direction of the X3 axis.
C             Defaults: A13=' ', B13=' ', C13=' '.
C     A31='string', B31='string', C31='string'... Strings in apostrophes
C             containing the names of the output files with
C             (N1-1)*N2*(N3-1) values of the elastic parameters averaged
C             in the direction the X3 axis, along the gridlines shifted
C             half a grid interval in the direction of the X1 axis.
C             Defaults: A31=' ', B31=' ', C31=' '.
C     DEN='string'... String in apostrophes containing the name of the
C             output ASCII file with N1*N2*N3 grid values of the
C             density.
C             Default: DEN=' '.
C     QFC='string'... String in apostrophes containing the name of the
C             output ASCII file with N1*N2*N3 grid values of the P-wave
C             quality factor.
C             Default: QFC=' '
C Example of data set FDSEP
C
C.......................................................................
C                                                
C This file consists of the following external procedures:
C     FDMOD...Main program allocating the memory for FDPAR.
C             FDMOD
C     FDPAR...Subroutine to generate the parameters (originally the main
C             program).
C             FDPAR
C     XLEG... Subroutine for computation of effective (mean) squared
C             velocities for a grid leg.
C             XLEG
C     GRPT... Subroutine for location of grid point (XPT,ZPT) with
C             respect to block IBL.
C             GRPT
C     TWOPT...Subroutine for location of a grid leg with respect to
C             boundary of block IB.
C             TWOPT
C
C=======================================================================
C
C     
C
      PROGRAM FDMOD
C
C Main program allocating the memory.
C
C Subroutines and external functions required:
      EXTERNAL FDPAR,RSEP1,RSEP3I
C     FDPAR...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
C
C-----------------------------------------------------------------------
C
C     Input data filenames and logical unit numbers:
      INTEGER LU
      PARAMETER (LU=1)
      CHARACTER*80 FSEP
C     LU...   Logical unit number to be used to read and write all data
C             files.
C     FSEP... The name of the input data file with FD parameters in SEP
C             format.
C
      INTEGER N1,N3,N,NR,NS,NGRID,I1,I2,I3,I4
C
C.......................................................................
C
C     Main input data:
C     Default:
      FSEP='fd.h'
C     Reading main input data:
      WRITE(*,'(A)') ' FDMOD: Enter name of input SEP file [''fd.h'']: '
      READ (*,*) FSEP
      WRITE(*,'(A)') '+FDMOD: Working...                               '
C
C     Reading all the data from input SEP parameter file:
      CALL RSEP1(LU1,FSEP)
C
C     Data specifying numbers of gridpoints:
      CALL RSEP3I('N1'  ,N1  ,1)
      CALL RSEP3I('N2'  ,N3  ,1)
C
C     Allocating the memory in array RAM:
      NGRID=15*N1*N3
      IF(NGRID.GT.MRAM) THEN
C       FDMOD-01
        CALL ERROR('FDMOD-01: Small dimension MRAM in ram.inc')
C       Dimension MRAM must be at least
C         15*N1*N3 ,
C       where N1 and N3 are parameters in the main input file.
      END IF
      N =N1*N3
C
      CALL FDPAR(LU,N,
     *        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))
C
      WRITE(*,'(A)') '+FDMOD: Done.                                    '
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE FDPAR(LU,KLMAX,
     *                 AHS,BHS,CHS,AVS,BVS,CVS,
     *                 AHL,BHL,CHL,AVL,BVL,CVL,DEN,QFC,WORK)
      INTEGER LU,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),QFC(KLMAX),WORK(KLMAX)
C
C Subroutine (originally main program) to read the input data and
C generate finite-difference effective parameters.
C
C Input:
C     LU...   Logical unit number to be used to read and write data
C             files.
C     KLMAX...Array dimension.
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
C No output.
C
C Subroutines and external functions required:
      EXTERNAL XLEG,GRPT,TWOPT
      EXTERNAL RSEP3T,FDGRID,FDOUT1
C     XLEG... This file.
C     GRPT... This file.
C     TWOPT.. This file.
C     RSEP3T..File 'sep.for'.
C     FDGRID..File 'fdaux.for'.
C     FDOUT1..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
C     Storage locations:
C
C     Maximum array dimensions:
      PARAMETER (MAXBL=50,MAXEDG=10)
C     MAXBL.. Max number of polygonal blocks.
C     MAXEDG..Max number of edges for any block.
C
C     Arrays:
      DIMENSION NOED(0:MAXBL),VS(0:MAXBL),PSRAT(0:MAXBL),HU(0:MAXBL)
      DIMENSION XED(MAXEDG,MAXBL),ZED(MAXEDG,MAXBL),NOVIC(MAXEDG,MAXBL)
C
      COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU
      COMMON /BLO2/ XED,ZED,NOVIC
C     NOBL...  Number of blocks.
C     NOED...  Number of edges of the block.
C     VS...    Shear-wave velocity (m/s) inside the block.
C     PSRAT..  P/S velocity ratio inside the block.
C     HU...    Density (kg/m**3) inside the block.
C     XED,ZED. X and Z coordinates (m) of the initial point of segment.
C     NOVIC... Serial number of the neigh. block (outside put 0).
C
      CHARACTER*80 FMOD
      CHARACTER *5 NAMDUM
C     FMOD... The name of the input data file with model parameters in
C             list directed (free) format.
C     NAMDUM... Auxiliary storage location.
C
C.......................................................................
C
C     Data specifying grid dimensions
      CALL FDGRID(KLMAX,KLMAX,KLMAX,NX,NZ,DX)
C
C     Reading data describing a model composed of homogeneous polygons:
      CALL RSEP3T('FDMOD',FMOD,' ')
      OPEN(LU,FILE=FMOD,STATUS='OLD')
C
      READ(1,'(1X,A5)') NAMDUM
C     READ(1,'(25X,I5)') NSYM
      NSYM=0
C
      READ(1,'(25X,I5)') NOBL
      DO 150 I=1,NOBL
        READ(1,'(25X,I5,3F10.0)') NOED(I),VS(I),PSRAT(I),HU(I)
        DO 100 J=1,NOED(I)
          READ(1,'(25X,2F10.0,I5)') XED(J,I),ZED(J,I),NOVIC(J,I)
  100   CONTINUE
  150 CONTINUE
C
      WRITE(*,'(A)') '+FDMOD: Computing effective parameters           '
C
C     Auxiliary parameters.
      LACC=0
C
      DO 310 I=1,NOBL
C       Attention: Now VS no more contains shear velocity but rigidity
C                  and PSRAT is squared Vp/Vs ratio.  The corresponding
C                  change was made in FLERIB,FBOB,FCORNB.
C       VS(I)=VS(I)*VS(I)
        VS(I)=VS(I)*VS(I)*HU(I)
        PSRAT(I)=PSRAT(I)*PSRAT(I)
  310 CONTINUE
C
      NOED(0)=0
      VS(0)=0.
      PSRAT(0)=1.
C
      DXA=DX
C
C     Vertical and horizontal parameters of medium for INTERlegs
C     (i.e. for legs between grid lines)
      NPUP=(NX-1)*(NZ-1)
      DO 430 IPU=1,NPUP
        NVER=NZ-1
        NP=IFIX((IPU-0.01)/NVER)
        XSOUR=DX/2.+FLOAT(NP)*DX
        NPOR=IPU-NP*NVER
        ZSOUR=FLOAT(NPOR-1)*DX
        CALL XLEG(0,XSOUR,ZSOUR,DXA,AL,BE,RH)
        AVS(IPU)=AL
        BVS(IPU)=BE
        CVS(IPU)=AL-2.*BE
        NHOR=NZ-1
        NP=IFIX((IPU-0.01)/NHOR)
        XSOUR=FLOAT(NP)*DX
        NPOR=IPU-NP*NHOR
        ZSOUR=DX/2.+FLOAT(NPOR-1)*DX
        CALL XLEG(1,XSOUR,ZSOUR,DXA,AL,BE,RH)
        AHS(IPU)=AL
        BHS(IPU)=BE
        CHS(IPU)=AL-2.*BE
  430 CONTINUE
C
C     Vertical parameters of medium for LARGE squares (= grid legs)
C     (for ignoring detailed computations, choose LACC=0)
      DO 450 K=2,NX-1
        XOLD=FLOAT(K-1)*DX
        LAU=(K-1)*NZ
        XNZ=FLOAT(NZ)
        DO 400 L=1,NZ-1
          ICEN=L+LAU
          JCEN=ICEN-IFIX(FLOAT(ICEN)/XNZ)
          ZOLD=FLOAT(L-1)*DX
          IF(L.LT.LACC) THEN
            CALL XLEG(0,XOLD,ZOLD,DXA/2.,AL1,BE1,RH1)
            ZNEW=ZOLD+DX/2.
            CALL XLEG(0,XOLD,ZNEW,DXA/2.,AL2,BE2,RH2)
            AL=(AL1+AL2)/2.
            BE=(BE1+BE2)/2.
          ELSE
            CALL XLEG(0,XOLD,ZOLD,DXA,AL,BE,RH)
          ENDIF
          AVL(JCEN)=AL
          BVL(JCEN)=BE
          CVL(JCEN)=AL-2.*BE
  400   CONTINUE
  450 CONTINUE
C
C     Horizontal parameters of medium  for LARGE squares (=grid legs)
C     (incl. density)
      DO 550 K=1,NX-1
        XOLD=FLOAT(K-1)*DX
        LAU=(K-1)*NZ
        DO 500 L=1,NZ-1
          ICEN=L+LAU
          ZOLD=FLOAT(L-1)*DX
          IF(L.LE.LACC) THEN
            CALL XLEG(1,XOLD,ZOLD,DXA/2.,AL1,BE1,RH1)
            XNEW=XOLD+DX/2.
            CALL XLEG(1,XNEW,ZOLD,DXA/2.,AL2,BE2,RH2)
            AL=(AL1+AL2)/2.
            BE=(BE1+BE2)/2.
          ELSE
            CALL XLEG(1,XOLD,ZOLD,DXA,AL,BE,RH)
          ENDIF
          AHL(ICEN)=AL
          BHL(ICEN)=BE
          CHL(ICEN)=AL-2.*BE
          DEN(ICEN)=RH
  500   CONTINUE
  550 CONTINUE
C
C     Vertical and horizontal parameters for legs in the strip between
C     the inner and outer nonreflecting boundaries (a continuation)
C     (incl. density for horizontal legs)
      DO 600 L=1,NZ-1
        ICEN=L
        ICEN2=ICEN+NZ
        JCEN=L
        JCEN2=JCEN+NZ-1
        AVL(JCEN)=AVL(JCEN2)
        BVL(JCEN)=BVL(JCEN2)
        CVL(JCEN)=CVL(JCEN2)
  600 CONTINUE
      LAU=(NX-1)*NZ
      DO 610 L=1,NZ-1
        ICEN=L+LAU
        JCEN=ICEN-IFIX(FLOAT(ICEN)/FLOAT(NZ))
        JCEN2=JCEN-NZ+1
        AVL(JCEN)=AVL(JCEN2)
        BVL(JCEN)=BVL(JCEN2)
        CVL(JCEN)=CVL(JCEN2)
  610 CONTINUE
      DO 620 K=1,NX-1
        ICEN=K*NZ
        ICEN2=ICEN-1
        AHL(ICEN)=AHL(ICEN2)
        BHL(ICEN)=BHL(ICEN2)
        CHL(ICEN)=CHL(ICEN2)
        DEN(ICEN)=DEN(ICEN2)
  620 CONTINUE
C
C     Halving the density at the free surface IZ=1:
      DO 70 I=1,NX*NZ,NZ
        DEN(I)=0.5*DEN(I)
   70 CONTINUE
C
C     No attenuation is considered in this version
      DO 80 I=1,NX*NZ
        QFC(I)=0.
   80 CONTINUE
C
C     Writing the gridded effective elastic parameters:
      CALL FDOUT1(LU,NX-1,NZ  ,1.,'A11',AHL,WORK)
      CALL FDOUT1(LU,NX-1,NZ  ,1.,'B11',BHL,WORK)
      CALL FDOUT1(LU,NX-1,NZ  ,1.,'C11',CHL,WORK)
      CALL FDOUT1(LU,NX  ,NZ-1,1.,'A33',AVL,WORK)
      CALL FDOUT1(LU,NX  ,NZ-1,1.,'B33',BVL,WORK)
      CALL FDOUT1(LU,NX  ,NZ-1,1.,'C33',CVL,WORK)
      CALL FDOUT1(LU,NX-1,NZ-1,1.,'A13',AHS,WORK)
      CALL FDOUT1(LU,NX-1,NZ-1,1.,'B13',BHS,WORK)
      CALL FDOUT1(LU,NX-1,NZ-1,1.,'C13',CHS,WORK)
      CALL FDOUT1(LU,NX-1,NZ-1,1.,'A31',AVS,WORK)
      CALL FDOUT1(LU,NX-1,NZ-1,1.,'B31',BVS,WORK)
      CALL FDOUT1(LU,NX-1,NZ-1,1.,'C31',CVS,WORK)
      CALL FDOUT1(LU,NX  ,NZ  ,1.,'DEN',DEN,WORK)
      CALL FDOUT1(LU,NX  ,NZ  ,1.,'QFC',QFC,WORK)
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE XLEG(KEY,XPT,ZPT,DEL,ALF,BET,RHO)
C
C      Purpose:
C      Effective (mean) squared velocities for a grid leg
C
C      Input:
C      KEY=0.....vertical leg
C        =1.....horizontal leg
C     XPT,ZPT...x,z coordinates of the grid point
C     DEL.......length of the leg
C
C      Output:
C     ALF,BET...P and S squared velocities for the leg
C
C
      DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50)
      DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50)
      COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU
      COMMON /BLO2/ XED,ZED,NOVIC
C
      RHO=0.
C
C      An artificial shortening of the leg to avoid numerical problems
C      !!!!!!Some models possible with 0.01 instead of this 0.03!!!!!!
C
      DELSA=DEL-0.03*DEL
cc    DELSA=DEL
cc    DELSA=DEL-0.0001*DEL
C
C      Location of the grid point with respect to a block
C
      DO 100 J=1,NOBL
      DIST=0.
      CALL GRPT(XPT,ZPT,J,IPRES,IVIC,DIST,ISEG,IBLVC,IORT)
      DIST=DIST/DELSA
      IF(IPRES.EQ.0) GOTO 100
      IBL=J
      GOTO 200
  100 CONTINUE
      WRITE(*,'(A,2F10.3/A)')
     *  '+FDMOD: Gridpoint not located successfully at x,z=',XPT,ZPT,' '
  200 CONTINUE
      IF(IPRES.EQ.2) GOTO 1000
C
C       The grid point inside block IBL
C
      RHO=HU(IBL)
      CALL TWOPT(XPT,ZPT,IBL,DELSA,KEY,0,ISG,IVC,PART)
      IF(ISG.NE.0) GOTO 400
C
C      No intersection of the leg with any boundary segment
C
  300 BET=VS(IBL)
      ALF=BET*PSRAT(IBL)
      RETURN
C
C      The leg is intersecting segment ISG
C      Locating the end point of the leg with respect to block IVC
C
  400 IF(KEY.EQ.0) THEN
       XEND=XPT
       ZEND=ZPT+DELSA
      ELSE
       XEND=XPT+DELSA
       ZEND=ZPT
      ENDIF
      DIS2=0.
      CALL GRPT(XEND,ZEND,IVC,IPR2,IVI2,DIS2,ISE2,IBL2,IOR2)
      DIS2=DIS2/DELSA
      IF(IPR2.EQ.0) GOTO 600
      IF(IPR2.EQ.1) GOTO 450
      IF(IPR2.EQ.2.AND.KEY.EQ.0.AND.IOR2.EQ.0) GOTO 500
      IF(IPR2.EQ.2.AND.KEY.EQ.1.AND.IOR2.EQ.1) GOTO 500
C
C      The end point located inside block IVC
C      (or on its boundary segment not coinciding with the leg)
C
  450 VSIN=VS(IBL)
      VSOU=VS(IVC)
      VPIN=VSIN*PSRAT(IBL)
      VPOU=VSOU*PSRAT(IVC)
c     BET=PART*VSIN+(1.-PART)*VSOU
c     ALF=PART*VPIN+(1.-PART)*VPOU
      BET=1./(PART/VSIN+(1.-PART)/VSOU)
      ALF=1./(PART/VPIN+(1.-PART)/VPOU)
      RETURN
C
C      The end point located on a vertical/horizontal boundary
C      between blocks IVC and IBL2
C
  500 VSIN=VS(IBL)
      VSOU1=VS(IVC)
      VSOU2=VS(IBL2)
      VSOU=(VSOU1+VSOU2)/2.
      VPIN=VSIN*PSRAT(IBL)
      VPOU=(VSOU1*PSRAT(IVC)+VSOU2*PSRAT(IBL2))/2.
c     BET=PART*VSIN+(1.-PART)*VSOU
c     ALF=PART*VPIN+(1.-PART)*VPOU
      BET=1./(PART/VSIN+(1.-PART)/VSOU)
      ALF=1./(PART/VPIN+(1.-PART)/VPOU)
      RETURN
C
C      The end point located outside block IVC,inside block IBLNEW
C
  600 DO 700 J=1,NOBL
      CALL GRPT(XEND,ZEND,J,IPR3,IVI3,DIS3,ISE3,IBL3,IOR3)
      IF(IPR3.EQ.0) GOTO 700
      IBLNEW=J
      GOTO 800
  700 CONTINUE
      WRITE(*,'(A,2F10.3/A)')
     *'+FDMOD: Gridpoint not located successfully at x,z=',XEND,ZEND,' '
C
C      The end point in block IBLNEW, the leg in blocks IBL, IVC, IBLNEW
C
  800 CALL TWOPT(XEND,ZEND,IBLNEW,DELSA,KEY,1,ISG2,IVC2, PART2)
      VSIN=VS(IBL)
      VSOU1=VS(IVC)
      VSOU2=VS(IBLNEW)
      VPIN=VSIN*PSRAT(IBL)
      VPOU1=VSOU1*PSRAT(IVC)
      VPOU2=VSOU2*PSRAT(IBLNEW)
c     BET=PART*VSIN+(1.-PART-PART2)*VSOU1+PART2*VSOU2
c     ALF=PART*VPIN+(1.-PART-PART2)*VPOU1+PART2*VPOU2
      BET=1./(PART/VSIN+(1.-PART-PART2)/VSOU1+PART2/VSOU2)
      ALF=1./(PART/VPIN+(1.-PART-PART2)/VPOU1+PART2/VPOU2)
      RETURN
C
C      The grid point located on the boundary of block IBL
C      (on segment ISEG, neighbouring to block IBLVC)
C       Locating the end point with respect to block IBL
C
 1000 CONTINUE
      if(ibl.eq.0) rho=hu(iblvc)
      if(iblvc.eq.0) rho=hu(ibl)
      if(ibl.ne.0.and.iblvc.ne.0) RHO=(HU(IBL)+HU(IBLVC))/2.
      IF(KEY.EQ.0) THEN
       XEND=XPT
       ZEND=ZPT+DELSA
      ELSE
       XEND=XPT+DELSA
       ZEND=ZPT
      ENDIF
      DIS2=0.
      CALL GRPT(XEND,ZEND,IBL,IPR2,IVI2,DIS2,ISE2,IBL2,IOR2)
      DIS2=DIS2/DELSA
      IF(IPR2.EQ.0) GOTO 1800
      IF(IPR2.EQ.1) GOTO 1600
      IF(IPR2.EQ.2.AND.KEY.EQ.0.AND.IOR2.EQ.0) GOTO 1700
      IF(IPR2.EQ.2.AND.KEY.EQ.1.AND.IOR2.EQ.1) GOTO 1700
C
C      The end point located inside block IBL
C      (or on its boundary segment not coinciding with the leg)
C
 1600 BET=VS(IBL)
      ALF=BET*PSRAT(IBL)
      RETURN
C
C      The leg is coinciding with a segment
C
 1700 VSIN=VS(IBL)
C     VSOU=VS(IBLVC)
      VSOU=VS(IBL2)
      VPIN=VSIN*PSRAT(IBL)
C     VPOU=VSOU*PSRAT(IBLVC)
      VPOU=VSOU*PSRAT(IBL2)
      BET=(VSIN+VSOU)/2.
      ALF=(VPIN+VPOU)/2.
      RETURN
C
C      The end point located outside block IBL
C     Locating the end point with respect to block IBLVC
C
 1800 DIS3=0.
      CALL GRPT(XEND,ZEND,IBLVC,IPR3,IVI3,DIS3,ISE3,IBL3,IOR3)
      DIS3=DIS3/DELSA
      IF(IPR3.EQ.0) GOTO 2100
      IF(IPR3.EQ.1) GOTO 1850
      IF(IPR3.EQ.2) GOTO 1900
C
C      The end point located inside block IBLVC
C
 1850 BET=VS(IBLVC)
      ALF=BET*PSRAT(IBLVC)
      RETURN
C
C      The end point located on segment ISE3 of block IBLVC
C
 1900 IF(KEY.EQ.0.AND.IOR3.EQ.0) GOTO 2000
      IF(KEY.EQ.1.AND.IOR3.EQ.1) GOTO 2000
C
 1950 BET=VS(IBLVC)
      ALF=BET*PSRAT(IBLVC)
      RETURN
C
C      The end point located on a vertical/horizontal boundary
C      between blocks IBLVC and IBL3
C
 2000 VSIN=VS(IBL)
      VSOU=(VS(IBLVC)+VS(IBL3))/2.
      VPIN=VS(IBL)*PSRAT(IBL)
      VPOU=(VS(IBLVC)*PSRAT(IBLVC)+VS(IBL3)*PSRAT(IBL3))/2.
c     BET=DIST*VSIN+(1.-DIST)*VSOU
c     ALF=DIST*VPIN+(1.-DIST)*VPOU
      BET=1./(DIST/VSIN+(1.-DIST)/VSOU)
      ALF=1./(DIST/VPIN+(1.-DIST)/VPOU)
      RETURN
C
C      The end point located outside block IBLVC, inside block IBLNEW
C
 2100 DO 2150 J=1,NOBL
      DIS5=0.
      CALL GRPT(XEND,ZEND,J,IPR5,IVI5,DIS5,ISE5,IBL5,IOR5)
      DIS5=DIS5/DELSA
      IF(IPR5.EQ.0) GOTO 2150
      IBLNEW=J
      GOTO 2170
 2150 CONTINUE
 2170 CONTINUE
      IF(KEY.EQ.0.AND.IORT.EQ.0) GOTO 2200
      IF(KEY.EQ.1.AND.IORT.EQ.1) GOTO 2200
      GOTO 2400
 2200 IF(IPR5.EQ.1) GOTO 2300
      IF(IPR5.EQ.2) GOTO 2350
C
C      The grid point located on a vertical/horizontal boundary
C      between blocks IBL and IBLVC
C      The end point inside block IBLNEW
C
 2300 VSIN=(VS(IBL)+VS(IBLVC))/2.
      VSOU=VS(IBLNEW)
      VPIN=(VS(IBL)*PSRAT(IBL)+VS(IBLVC)*PSRAT(IBLVC))/2.
      VPOU=VSOU*PSRAT(IBLNEW)
c     BET=DIST*VSIN+(1.-DIST)*VSOU
c     ALF=DIST*VPIN+(1.-DIST)*VPOU
      BET=1./(DIST/VSIN+(1.-DIST)/VSOU)
      ALF=1./(DIST/VPIN+(1.-DIST)/VPOU)
      RETURN
C
C      The grid point located on a vertical/horizontal boundary
C      between blocks IBL and IBLVC
C      The end point located on a vertical/horizontal boundary
C      between blocks IBLNEW and IBL5
C
 2350 VSIN=(VS(IBL)+VS(IBLVC))/2.
      VSOU=(VS(IBLNEW)+VS(IBL5))/2.
      VPIN=(VS(IBL)*PSRAT(IBL)+VS(IBLVC)*PSRAT(IBLVC))/2.
      VPOU=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2.
c     BET=DIST*VSIN+(1.-DIST)*VSOU
c     ALF=DIST*VPIN+(1.-DIST)*VPOU
      BET=1./(DIST/VSIN+(1.-DIST)/VSOU)
      ALF=1./(DIST/VPIN+(1.-DIST)/VPOU)
      RETURN
C
C      The end point located inside block IBLNEW,
C      the leg in blocks IBLNEW and IVC4 (where IVC4=IBL or IBLVC)
C
 2400 IF(IPR5.EQ.1) GOTO 2600
      IF(IPR5.EQ.2.AND.KEY.EQ.0.AND.IOR5.EQ.0) GOTO 2800
      IF(IPR5.EQ.2.AND.KEY.EQ.1.AND.IOR5.EQ.1) GOTO 2800
 2600 CALL TWOPT(XEND,ZEND,IBLNEW,DELSA,KEY,1,ISG4,IVC4,PART4)
      VSIN=VS(IVC4)
      VSOU=VS(IBLNEW)
      VPIN=VSIN*PSRAT(IVC4)
      VPOU=VSOU*PSRAT(IBLNEW)
c     BET=(1.-PART4)*VSIN+PART4*VSOU
c     ALF=(1.-PART4)*VPIN+PART4*VPOU
c  some models with inclined boundaries ending at left or right sides
c  require the following modification:
c     if(vsin.eq.0.) then
c      BET=1./(PART4/VSOU)
c      ALF=1./(PART4/VPOU)
c      RETURN
c     endif
      BET=1./((1.-PART4)/VSIN+PART4/VSOU)
      ALF=1./((1.-PART4)/VPIN+PART4/VPOU)
      RETURN
 2800 CALL TWOPT(XEND,ZEND,IBL,DELSA,KEY,1,ISG5,IVC5,PART5)
      IF(ISG5.NE.0) THEN
       IVXX=IBL
       PARTX=PART5
       GOTO 2900
      ELSE
      CALL TWOPT(XEND,ZEND,IBLVC,DELSA,KEY,1,ISG6,IVC6,PART6)
      IF(ISG6.NE.0) THEN
       IVXX=IBLVC
       PARTX=PART6
       GOTO 2900
      ELSE
c
 2850  BET=(VS(IBLNEW)+VS(IBL5))/2.
       ALF=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2.
       RETURN
      ENDIF
      ENDIF
 2900 VSIN=VS(IVXX)
      VSOU=(VS(IBLNEW)+VS(IBL5))/2.
      VPIN=VSIN*PSRAT(IVXX)
      VPOU=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2.
c     BET=(1.-PARTX)*VSIN+PARTX*VSOU
c     ALF=(1.-PARTX)*VPIN+PARTX*VPOU
      BET=1./((1.-PARTX)/VSIN+PARTX/VSOU)
      ALF=1./((1.-PARTX)/VPIN+PARTX/VPOU)
      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE GRPT(XPT,ZPT,IBL,IPRES,IVIC,DIST,ISEG,IBLVC,IORT)
C
C      Purpose:
C      Location of grid point (XPT,ZPT) with respect to block IBL
C
C      Input:
C      XPT,ZPT......x,z coordinates of the grid point
C      IBL..........the block to be tested
C
C      Output:
C      IPRES=0......the grid point located outside the block
C          =1......the grid point located inside the block
C          =2......the grid point located on the block boundary
C
C      IVIC=0.......the grid point not in vicinity (.LT.DX) of any edge
C      .ne.0.......the grid point in vicinity of edge IVIC
C
C      DIST.........distance from the edge (for IVIC.ne.0)
C
C      ISEG.........the grid point on segment ISEG (for IPRES=2)
C
C      IBLVC........the block neighbouring to ISEG (for IPRES=2)
C
C      IORT=0.......a vertical segment (for IPRES=2)
C         =1.......a horizontal segment (for IPRES=2)
C         =2.......an oblique segment (for IPRES=2)
C
C
      DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50)
      DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50)
      COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU
      COMMON /BLO2/ XED,ZED,NOVIC
C
      I=IBL
      IF(I.EQ.0) THEN
       IPRES=0
       RETURN
      ENDIF
      NEDG=NOED(I)
      DO 200 J=1,NEDG
C
C      Unit position vector (PX,PZ) of grid point (XPT,ZPT)
C      with respect to the origin of J-th segment of I-th block
C
      PX=XPT-XED(J,I)
      PZ=ZPT-ZED(J,I)
      SIZ=SQRT(PX*PX+PZ*PZ)
      IF(ABS(SIZ).LT.0.00001) GOTO 300
      PX=PX/SIZ
      PZ=PZ/SIZ
C
C      Unit vector (BX,BZ) of J-th segment of I-th block
C
      IF(J.LT.NEDG) THEN
       BX=XED(J+1,I)-XED(J,I)
       BZ=ZED(J+1,I)-ZED(J,I)
      ELSE
       BX=XED(1,I)-XED(J,I)
       BZ=ZED(1,I)-ZED(J,I)
      ENDIF
      SIZ=SQRT(BX*BX+BZ*BZ)
      BX=BX/SIZ
      BZ=BZ/SIZ
C
C      Y-component  (=PROD)  of the vector product (PX,PZ)x(BX,BZ)
C
      PROD=PZ*BX-PX*BZ
      IF(PROD.LT.-0.00001) GOTO 200
      IF(ABS(PROD).LT.0.00001) GOTO 300
C
C       The grid point located outside the block
C
  250 IPRES=0
      GOTO 600
C
C      The grid point possibly located on J-th segment
C
  300 IPRES=2
      ISEG=J
      IBLVC=NOVIC(J,I)
C
C      End points of J-th segment, orientation of the segment
C
      X1=XED(J,I)
      Z1=ZED(J,I)
      IF(J.LT.NEDG) THEN
       X2=XED(J+1,I)
       Z2=ZED(J+1,I)
      ELSE
       X2=XED(1,I)
       Z2=ZED(1,I)
      ENDIF
      SEGLEN=SQRT((X2-X1)**2+(Z2-Z1)**2)
      IORT=2
      IF(ABS(X1-X2).LT.0.00001) IORT=0
      IF(ABS(Z1-Z2).LT.0.00001) IORT=1
C
C      Distances of the grid point from the end points
C
      DIST1=SQRT((XPT-X1)**2+(ZPT-Z1)**2)
      DIST2=SQRT((XPT-X2)**2+(ZPT-Z2)**2)
      ROZ1=DIST1-SEGLEN
      ROZ2=DIST2-SEGLEN
C     IF(DIST1.GT.SEGLEN.OR.DIST2.GT.SEGLEN) GOTO 250
      IF(ROZ1.GT.0.0001.OR.ROZ2.GT.0.0001) GOTO 250
C
C      The grid point certainly located on J-th segment
C
      IVIC=0
      IF(DIST1.LT.DX) GOTO 400
      IF(DIST2.LT.DX) GOTO 500
      GOTO 600
  400   IVIC=J
      DIST=DIST1
      GOTO 600
  500      IF(J.LT.NEDG) THEN
       IVIC=J+1
      ELSE
       IVIC=1
      ENDIF
      DIST=DIST2
      GOTO 600
  200   CONTINUE
C
C      The grid point 'on the left' with respect to all segments,
C      i.e., the grid point located inside I-th block
C
      IPRES=1
C
C      Distances from all edges of the block
C      (The program cannot detect proximity of two edges simultan.)
C
      IVIC=0
      DO 550 J=1,NEDG
      DISTX=SQRT((XPT-XED(J,I))**2+(ZPT-ZED(J,I))**2)
      IF(DISTX.GT.DX) GOTO 550
      IVIC=J
      DIST=DISTX
      GOTO 600
  550      CONTINUE
  600      RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
C
C=======================================================================
C
C     
C
      SUBROUTINE TWOPT(X,Z,IB,DELSA,KY,KYOR,ISG,IVC,PART)
C
C      Purpose:
C      Location of a grid leg with respect to boundary of block IB
C
C      Input:
C      X,Z........x,z coordinates of the grid point
C      IB.........the block to which the grid point is belonging
C      DELSA........length of the leg
C      KY=0.......vertical leg
C       =1.......horizontal leg
C      KYOR=0.....leg oriented down, or right, of the grid point
C          =1.....up or left
C
C      Output:
C      ISG=0......no intersection of the leg with the block boundary
C       .ne.0....intersection with ISG-th segment
C      IVC........the neighbouring block (for ISG.ne.0)
C      PART.......a relative distance of the grid point from
C                the intersection (for ISG.ne.0)
C
C
      DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50)
      DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50)
      COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU
      COMMON /BLO2/ XED,ZED,NOVIC
C
      IF(IB.EQ.0) THEN
       ISG=0
       RETURN
      ENDIF
C
C      End point (XX,ZZ) of the leg (the leg from (X,Z) to (XX,ZZ))
C
      IF(KY.EQ.0) GOTO 100
      IF(KY.EQ.1) GOTO 200
  100 XX=X
      IF(KYOR.EQ.0) THEN
       ZZ=Z+DELSA
      ELSE
       ZZ=Z-DELSA
      ENDIF
      GOTO 300
  200 ZZ=Z
      IF(KYOR.EQ.0) THEN
       XX=X+DELSA
      ELSE
       XX=X-DELSA
      ENDIF
C      Intersection of the leg with the block boundary
C      (legs coinciding with boundaries are not allowed)
C
  300 NEDG=NOED(IB)
      DO 500 J=1,NEDG
C
C      The end points (X1,Z1) and (X2,Z2) of J-th segment
C
      X1=XED(J,IB)
      Z1=ZED(J,IB)
      IF(J.LT.NEDG) THEN
       X2=XED(J+1,IB)
       Z2=ZED(J+1,IB)
      ELSE
       X2=XED(1,IB)
       Z2=ZED(1,IB)
      ENDIF
      IF(KY.EQ.0) GOTO 600
      IF(KY.EQ.1) GOTO 700
  600 IF(X1.LT.X2) THEN
       XMIN=X1
       XMAX=X2
      ELSE
       XMIN=X2
       XMAX=X1
      ENDIF
      IF(X.LT.XMIN.OR.X.GT.XMAX) GOTO 500
      ZSEC=Z1+(X-X1)*(Z2-Z1)/(X2-X1)
      IF(KYOR.EQ.0) THEN
       ZBOT=Z
       ZTOP=ZZ
      ELSE
       ZBOT=ZZ
       ZTOP=Z
      ENDIF
      IF(ZSEC.LT.ZBOT.OR.ZSEC.GT.ZTOP) GOTO 500
      ISG=J
      IVC=NOVIC(J,IB)
      PART=ABS(ZSEC-Z)/DELSA
      GOTO 1000
  700 IF(Z1.LT.Z2) THEN
       ZMIN=Z1
       ZMAX=Z2
      ELSE
       ZMIN=Z2
       ZMAX=Z1
      ENDIF
      IF(Z.LT.ZMIN.OR.Z.GT.ZMAX) GOTO 500
      XSEC=X1+(Z-Z1)*(X2-X1)/(Z2-Z1)
      IF(KYOR.EQ.0) THEN
       XBOT=X
       XTOP=XX
      ELSE
       XBOT=XX
       XTOP=X
      ENDIF
      IF(XSEC.LT.XBOT.OR.XSEC.GT.XTOP) GOTO 500
      ISG=J
      IVC=NOVIC(J,IB)
      PART=ABS(XSEC-X)/DELSA
      GOTO 1000
  500 CONTINUE
      ISG=0
 1000 RETURN
      END
C     DEBUG SUBCHK
C     END DEBUG
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 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C