C
C Program MODFD to calculate effective elastic parameters (harmonic
C averages of elastic parameters) for 2-D elastic finite differences
C in models specified by means of the MODEL package
C
C Version: 5.40
C Date: 2000, April 9
C
C Coded by Ludek Klimes
C             Department of Geophysics, Charles University Prague
C             Ke Karlovu 3,  121 16  Praha 2,  Czech Republic
C             E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C Contents:
C     Description of data files
C     Description of effective elast. parameters
C     List of external procedures
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 data for the 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     MODEL='string'... Name of the input file with the data specifying
C             the model.  Must be specified.
C             Description of MODEL
C             Example of data set MODEL
C             Default: MODEL=' '
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     N2=positive integer... Number of gridpoints of the basic grid
C             along the X2 axis.  Default: N2=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     O2=real... X2 coordinate of the origin of the grid. Default: O2=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     D2=real... Grid spacing along the X2 axis.  Default: D2=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     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             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     VP='string'... String in apostrophes containing the name of the
C             output ASCII file with N1*N2*N3 grid values of the P wave
C             velocity.
C             Default: VP=' '.
C     VS='string'... String in apostrophes containing the name of the
C             output ASCII file with N1*N2*N3 grid values of the S wave
C             velocity.
C             Default: VS=' '.
C Example of data set FDSEP
C
C.......................................................................
C                                                 
C Description of effective elastic parameters:
C
C Harmonic averages of the following elastic parameters are evaluated:
C     A = 1 / average( 1/(lambda+2*mu) )
C     B = 1 / average( 1/mu     )
C     C = 1 / average( 1/lambda )
C Elastic parameters A,B,C are averaged between gridpoints along legs
C in direction X1 (yielding parameters A1j,B1j,C1j), direction X2
C (yielding parameters A2j,B2j,C2j), and direction X3 (yielding
C parameters A3j,B3j,C3j).  Parameters Aij,Bij,Cij correspond to the
C grid legs between gridpoints if i=j.  For j different from i, the
C interlegs are considered.  The interlegs Aij,Bij,Cij are centred
C between the grid legs Aii,Bii,Cii with respect to the Xj axis.
C For example, A11,B11,C11 correspond to the grid legs in the direction
C of the X1 axis, whereas A13,B13,C13 correspond to the interlegs in the
C direction of the X1 axis, shifted vertically with respect to the
C former grid legs if the X3 axis is vertical.
C
C The density is discretized at gridpoints of the basic grid N1*N2*N3:
C     DEN = density
C
C The P-wave quality factor is discretized at gridpoints of the basic
C grid N1*N2*N3 to be substituted for the finite-difference quantity
C     QFC = Q/f,
C where Q is the quality factor and f is frequency.
C
C For the legs or gridpoints situated at the structural interfaces, the
C parameters are calculated independently at both the sides of the
C interface.  The results are then averaged arithmetically.
C Please, note that this is the difference with respect to program FDMOD
C by Jiri Zahradnik.  The arithmetic averages should be calculated first
C and harmonically averaged along grid legs.  Moreover, when calculating
C the density at the edges, the slopes of interfaces should be taken
C into account instead of averaging the values from different blocks
C with the same weight.
C
C The number of output values of each average elastic parameter is thus:
C     for A11,B11,C11:              (N1-1)* N2   * N3    values,
C     for A12,B12,C12,A21,B21,C21:  (N1-1)*(N2-1)* N3    values,
C     for A22,B22,C22:               N1   *(N2-1)* N3    values,
C     for A13,B13,C13,A31,B31,C31:  (N1-1)* N2   *(N3-1) values,
C     for A23,B23,C23,A32,B32,C32:   N1   *(N2-1)*(N3-1) values,
C     for A33,B33,C33:               N1   * N2   *(N3-1) values,
C     for DEN,QFC:                   N1   * N2   * N3    values,
C where the dimensions of the basic grid are N1*N2*N3.
C
C.......................................................................
C                                                
C This file consists of the following external procedures:
C     MODFD...Main program.
C             MODFD
C     LEGS... Subroutine to evaluate reciprocal values of elastic
C             parameters A=lambda+2*mu, B=mu and C=lambda along
C             the legs of a given grid.
C             LEGS
C     ELPAR...Subroutine to evaluate reciprocal values of elastic
C             parameters A=lambda+2*mu, B=mu and C=lambda at a
C             given point.
C             ELPAR
C
C=======================================================================
C
C     
C
      PROGRAM MODFD
C
C.......................................................................
C
C     Input and output data files:
      CHARACTER*80 FSEP,FMOD
      CHARACTER*80 FA11,FB11,FC11,FA12,FB12,FC12,FA13,FB13,FC13
      CHARACTER*80 FA21,FB21,FC21,FA22,FB22,FC22,FA23,FB23,FC23
      CHARACTER*80 FA31,FB31,FC31,FA32,FB32,FC32,FA33,FB33,FC33
      CHARACTER*80 FDEN,FQFC,FVP,FVS
      INTEGER LU
      PARAMETER (LU=1)
C
C     Input data and parameters:
      INTEGER N1,N2,N3
      REAL O1,O2,O3,D1,D2,D3,S1,S2,S3
C
C     Auxiliary storage locations:
      INTEGER IOUT,NOUT,NALL
C
C.......................................................................
C
C     Reading main input data:
      FSEP='fd.h'
      WRITE(*,'(A)') ' MODFD: Enter names of SEP file [''fd.h'']: '
      READ (*,*) FSEP
      WRITE(*,'(A)') '+MODFD: Working...                          '

C     Reading all the data from file FSEP to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FSEP)
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
C     Shifting grid origin by half the grid intervals for interlegs
      S1=O1+0.5*D1
      S2=O2+0.5*D2
      S3=O3+0.5*D3
C
C     Reading the data for the model:
      CALL RSEP3T('MODEL',FMOD,' ')
      OPEN(LU,FILE=FMOD,STATUS='OLD')
      CALL MODEL1(LU)
      CLOSE(LU)
C
C     Recalling the output filenames:
      CALL RSEP3T('A11',FA11,' ')
      CALL RSEP3T('B11',FB11,' ')
      CALL RSEP3T('C11',FC11,' ')
      CALL RSEP3T('A12',FA12,' ')
      CALL RSEP3T('B12',FB12,' ')
      CALL RSEP3T('C12',FC12,' ')
      CALL RSEP3T('A13',FA13,' ')
      CALL RSEP3T('B13',FB13,' ')
      CALL RSEP3T('C13',FC13,' ')
      CALL RSEP3T('A21',FA21,' ')
      CALL RSEP3T('B21',FB21,' ')
      CALL RSEP3T('C21',FC21,' ')
      CALL RSEP3T('A22',FA22,' ')
      CALL RSEP3T('B22',FB22,' ')
      CALL RSEP3T('C22',FC22,' ')
      CALL RSEP3T('A23',FA23,' ')
      CALL RSEP3T('B23',FB23,' ')
      CALL RSEP3T('C23',FC23,' ')
      CALL RSEP3T('A31',FA31,' ')
      CALL RSEP3T('B31',FB31,' ')
      CALL RSEP3T('C31',FC31,' ')
      CALL RSEP3T('A32',FA32,' ')
      CALL RSEP3T('B32',FB32,' ')
      CALL RSEP3T('C32',FC32,' ')
      CALL RSEP3T('A33',FA33,' ')
      CALL RSEP3T('B33',FB33,' ')
      CALL RSEP3T('C33',FC33,' ')
      CALL RSEP3T('DEN',FDEN,' ')
      CALL RSEP3T('QFC',FQFC,' ')
      CALL RSEP3T('VP' ,FVP ,' ')
      CALL RSEP3T('VS' ,FVS ,' ')
C
C.......................................................................
C
C     Calculating elastic parameters:
C
C     Preparing screen output with the progress report:
C     Percents of calculated values
      IOUT=0
C     Number of calculated values
      NOUT=0
C     Number of values to be calculated
      NALL=0
      IF(FA11.NE.' '.OR.FB11.NE.' '.OR.FC11.NE.' '
     *              .OR.FDEN.NE.' '.OR.FQFC.NE.' '
     *              .OR.FVP .NE.' '.OR.FVS .NE.' ') THEN
        NALL=NALL+(N1-1)* N2   * N3
      END IF
      IF(FA12.NE.' '.OR.FB12.NE.' '.OR.FC12.NE.' ') THEN
        NALL=NALL+(N1-1)*(N2-1)* N3
      END IF
      IF(FA13.NE.' '.OR.FB13.NE.' '.OR.FC13.NE.' ') THEN
        NALL=NALL+(N1-1)* N2   *(N3-1)
      END IF
      IF(FA21.NE.' '.OR.FB21.NE.' '.OR.FC21.NE.' ') THEN
        NALL=NALL+(N1-1)*(N2-1)* N3
      END IF
      IF(FA22.NE.' '.OR.FB22.NE.' '.OR.FC22.NE.' ') THEN
        NALL=NALL+ N1   *(N2-1)* N3
      END IF
      IF(FA23.NE.' '.OR.FB23.NE.' '.OR.FC23.NE.' ') THEN
        NALL=NALL+ N1   *(N2-1)*(N3-1)
      END IF
      IF(FA31.NE.' '.OR.FB31.NE.' '.OR.FC31.NE.' ') THEN
        NALL=NALL+(N1-1)* N2   *(N3-1)
      END IF
      IF(FA32.NE.' '.OR.FB32.NE.' '.OR.FC32.NE.' ') THEN
        NALL=NALL+ N1   *(N2-1)*(N3-1)
      END IF
      IF(FA33.NE.' '.OR.FB33.NE.' '.OR.FC33.NE.' ') THEN
        NALL=NALL+ N1   * N2   *(N3-1)
      END IF
      NALL=NALL*4
C
C     Initial screen output:
      WRITE(*,'(A,I3,A)') '+',IOUT,
     * ' per cent evaluated                                            '
C
C     Calculating elastic parameters along grid legs in direction X1:
      CALL LEGS(LU,FA11,FB11,FC11,FDEN,FQFC,FVP,FVS,
     *              1,N1,O1,D1,2,N2  ,O2,D2,3,N3  ,O3,D3,IOUT,NOUT,NALL)
      CALL LEGS(LU,FA12,FB12,FC12,' ',' ',' ',' ',
     *              1,N1,O1,D1,2,N2-1,S2,D2,3,N3  ,O3,D3,IOUT,NOUT,NALL)
      CALL LEGS(LU,FA13,FB13,FC13,' ',' ',' ',' ',
     *              1,N1,O1,D1,2,N2  ,O2,D2,3,N3-1,S3,D3,IOUT,NOUT,NALL)
C
C     Calculating elastic parameters along grid legs in direction X2:
      CALL LEGS(LU,FA21,FB21,FC21,' ',' ',' ',' ',
     *              2,N2,O2,D2,1,N1-1,S1,D1,3,N3  ,O3,D3,IOUT,NOUT,NALL)
      CALL LEGS(LU,FA22,FB22,FC22,' ',' ',' ',' ',
     *              2,N2,O2,D2,1,N1  ,O1,D1,3,N3  ,O3,D3,IOUT,NOUT,NALL)
      CALL LEGS(LU,FA23,FB23,FC23,' ',' ',' ',' ',
     *              2,N2,O2,D2,1,N1  ,O1,D1,3,N3-1,S3,D3,IOUT,NOUT,NALL)
C
C     Calculating elastic parameters along grid legs in direction X3:
      CALL LEGS(LU,FA31,FB31,FC31,' ',' ',' ',' ',
     *              3,N3,O3,D3,1,N1-1,S1,D1,2,N2  ,O2,D2,IOUT,NOUT,NALL)
      CALL LEGS(LU,FA32,FB12,FC32,' ',' ',' ',' ',
     *              3,N3,O3,D3,1,N1  ,O1,D1,2,N2-1,S2,D2,IOUT,NOUT,NALL)
      CALL LEGS(LU,FA33,FB33,FC33,' ',' ',' ',' ',
     *              3,N3,O3,D3,1,N1  ,O1,D1,2,N2  ,O2,D2,IOUT,NOUT,NALL)
C
      WRITE(*,'(A)') '+MODFD: Done.                                    '
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE LEGS(LU,FA,FB,FC,FDEN,FQFC,FVP,FVS,
     *               K1,N1,O1,D1,K2,N2,O2,D2,K3,N3,O3,D3,IOUT,NOUT,NALL)
C
      CHARACTER*(*) FA,FB,FC,FDEN,FQFC,FVP,FVS
      INTEGER LU,K1,N1,K2,N2,K3,N3,IOUT,NOUT,NALL
      REAL O1,D1,O2,D2,O3,D3
C
C Subroutine to evaluate reciprocal values of elastic parameters
C A=lambda+2*mu, B=mu and C=lambda along the legs of a given grid.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C     FA,FB,FC,FDEN,FQFC... Output file names.  File of blank names will
C             not be generated.  FDEN and FQFC may be non-blank only if
C             K1=1, K2=2, K3=3.
C     K1...   Index of the coordinate axis parallel with the grid legs.
C     K2,K3...Indices of the other two axes in an ascending order.
C     N1...   Number of grid legs along a gridline increased by 1 (i.e.,
C             the number of grid points in the corresponding direction).
C     N2,N3...Numbers of gridlines.
C     O1,O2,O3... K1,K2,K3-th coordinates of the first gridpoint of the
C             first gridline.
C     D1,D2,D3... Grid intervals along the K1,K2,K3-th coordinate axes.
C     IOUT,NOUT,NALL... Integers to control the screen output with the
C             progress report.
C             IOUT: Percents of calculated values.
C             NOUT: Number of calculated values.
C             NALL: Number of values to be calculated.
C
C Output:
C     IOUT,NOUT,NALL... Updated input values.
C
C Date: 2000, April 9
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/ to allocate large array RAM:
      INCLUDE 'ram.inc'
C     ram.inc
C
C-----------------------------------------------------------------------
C
C     Input data and parameters:
      CHARACTER*80 FORM
      REAL ERR,W1,W2,W3
C
C     Averaged parameters:
      LOGICAL LA,LB,LC
      REAL AOLD,BOLD,COLD,ANEW,BNEW,CNEW,A,B,C,DEN,QFC,VP,VS
C
C     Auxiliary storage locations:
      INTEGER NN1,N1N2N3,J1,J2,J3,I1,I2,I3,IW2,IW3,IA,IB,IC,ID,IE,IF,IG
      INTEGER IH,ITER
C     Indices of blocks, surfaces, etc.:
      INTEGER IY(8),ISBOLD,ICBOLD,ISBNEW,ICBNEW,ISRF,IDUMMY(1)
C     Coordinates:
      REAL XOLD(3),XNEW(3),XTMP(3),X(3),XOLD1,XNEW1,XTMP1,X1,X2,X3
C     Directional vectors etc.:
      REAL DOLD(3),DNEW(3),DTMP(3),D(3),DIST,H,DUMMY(1)
C
C.......................................................................
C
C     Upper error bound for determination of the interface position:
      ERR=0.001000*ABS(D1)
      W1=AMAX1(ABS(O1),ABS(O1+FLOAT(N1-1)*D1))
      W2=AMAX1(ABS(O2),ABS(O2+FLOAT(N2-1)*D2))
      W3=AMAX1(ABS(O3),ABS(O3+FLOAT(N3-1)*D3))
      W1=AMAX1(0.0000001*W1,0.000001*ABS(D1))
      W2=AMAX1(0.0000001*W2,0.000001*ABS(D2))
      W3=AMAX1(0.0000001*W3,0.000001*ABS(D3))
      W1=SIGN(W1,D1)
C
C     Check whether the output files should be generated:
      IF(FA.EQ.' '.AND.FB  .EQ.' '.AND.FC  .EQ.' '
     *            .AND.FDEN.EQ.' '.AND.FQFC.EQ.' '
     *            .AND.FVP .EQ.' '.AND.FVS .EQ.' ') THEN
        RETURN
      END IF
      NN1=N1-1
      N1N2N3=NN1*N2*N3
C
C     Increments of indices in array RAM:
      IF(K1.EQ.1.AND.K2.EQ.2.AND.K3.EQ.3) THEN
        J1=1
        J2=NN1
        J3=0
      ELSE IF(K1.EQ.2.AND.K2.EQ.1.AND.K3.EQ.3) THEN
        J1=N2
        J2=1
        J3=NN1*N2-N2
      ELSE IF(K1.EQ.3.AND.K2.EQ.1.AND.K3.EQ.2) THEN
        J1=N2*N3
        J2=1
        J3=0
      ELSE
        RETURN
      END IF
C
C.......................................................................
C
C     Calculating elastic parameters along grid legs in direction K1:
      IA=1
      IB=IA+N1N2N3
      IC=IB+N1N2N3
      ID=IC+N1N2N3
      IE=ID+N1*N2*N3
      IF=IE+N1*N2*N3
      IF(FVP.EQ.' ') THEN
        IG=IF+1
      ELSE
        IG=IF+N1*N2*N3
      END IF
      IF(FVS.EQ.' ') THEN
        IH=IG+1
      ELSE
        IH=IG+N1*N2*N3
      END IF
      IF(IH.GT.MRAM) THEN
C       MODFD-01
        CALL ERROR('MODFD-01: Insufficient memory')
C       Increase the dimension MRAM in include file RAM.INC
      END IF
      ISBOLD=0
      DOLD(1)=0.
      DOLD(2)=0.
      DOLD(3)=0.
      DOLD(K1)=1.
      DNEW(1)=0.
      DNEW(2)=0.
      DNEW(3)=0.
      DNEW(K1)=0.
      DTMP(1)=0.
      DTMP(2)=0.
      DTMP(3)=0.
      DTMP(K1)=0.
C
      DO 33 I3=0,N3-1
        X3=O3+D3*FLOAT(I3)
        DO 32 I2=0,N2-1
          X2=O2+D2*FLOAT(I2)
          DO 23 IW3=-1,1,2
            XNEW(K3)=X3+W3*FLOAT(IW3)
            XOLD(K3)=XNEW(K3)
            XTMP(K3)=XNEW(K3)
            DO 22 IW2=-1,1,2
              XNEW(K2)=X2+W2*FLOAT(IW2)
              XOLD(K2)=XNEW(K2)
              XTMP(K2)=XNEW(K2)
C
              DO 21 I1=0,NN1
                XNEW(K1)=O1+D1*FLOAT(I1)-W1
                XOLD(K1)=O1+D1*FLOAT(I1-1)+W1
                IF(IW2.EQ.-1.AND.IW3.EQ.-1) THEN
                  RAM(IA)=0.
                  RAM(IB)=0.
                  RAM(IC)=0.
                  RAM(ID)=0.
                  RAM(IE)=0.
                  IF(FVP.NE.' ') THEN
                    RAM(IF)=0.
                  END IF
                  IF(FVS.NE.' ') THEN
                    RAM(IG)=0.
                  END IF
                END IF
C
C               Elastic parameters A,B,C are harmonically averaged along
C               four grid legs close to the leg between gridpoints
C               XOLD and XNEW.  The resulting elastic parameters are the
C               arithmetic averages of the respective four harmonic
C               averages.
C
C               Initial values for averaging:
                A=0.
                B=0.
                C=0.
                H=0.
C
C               Evaluation of block indices and material parameters:
                CALL BLOCK(XNEW,0,0,ISRF,ISBNEW,ICBNEW)
                CALL ELPAR
     *               (ICBNEW,XNEW,LA,LB,LC,ANEW,BNEW,CNEW,DEN,QFC,VP,VS)
C
C               Density and P-wave quality factor at the gridpoint:
                RAM(ID)=RAM(ID)+0.125*DEN
                RAM(IE)=RAM(IE)+0.125*QFC
                IF(FVP.NE.' ') THEN
                  RAM(IF)=RAM(IF)+0.125*VP
                END IF
                IF(FVS.NE.' ') THEN
                  RAM(IG)=RAM(IG)+0.125*VS
                END IF
C
C               Harmonic average of elastic parameters:
                IF(I1.GT.0) THEN
                  ITER=0
   10             CONTINUE
                  IF(ISBNEW.NE.ISBOLD) THEN
                    ITER=ITER+1
                    IF(ITER.GT.100) THEN
C                MODFD-02
                      CALL ERROR
     *                         ('MODFD-02: More than 100 intersections')
C                     More than 100 points of intersection of the
C                     gridline element with the structural interfaces.
C                     Check the input data and then contact the author.
                    END IF
C
C                   Determining the interface between points XOLD, XNEW:
                    IY(4)=ISBOLD
                    IY(5)=ICBOLD
                    IY(6)=0
                    XTMP(K1)=XNEW(K1)
                    XNEW1=XNEW(K1)
                    XOLD1=XOLD(K1)
                    XTMP1=XTMP(K1)
                    CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERR,
     *                       XOLD1,XOLD1,XOLD,DOLD,XNEW1,XNEW,DNEW,
     *                             XTMP1,XTMP,DTMP,X1,X,D)
                    IF(IY(6).EQ.0) THEN
C                     No more structural interface:
                      GO TO 20
                    END IF
C                   X is the point of intersection with the interface.
C
C                   Contributions to the average:
                    DIST=X(K1)-XOLD(K1)
                    A=A+DIST*AOLD
                    B=B+DIST*BOLD
                    C=C+DIST*COLD
                    XOLD(K1)=X(K1)
                    CALL ELPAR
     *               (ICBOLD,XOLD,LA,LB,LC,AOLD,BOLD,COLD,DEN,QFC,VP,VS)
                    A=A+DIST*AOLD
                    B=B+DIST*BOLD
                    C=C+DIST*COLD
                    H=H+DIST
                    ISBOLD=IY(7)
                    ICBOLD=IY(8)
                    CALL ELPAR
     *               (ICBOLD,XOLD,LA,LB,LC,AOLD,BOLD,COLD,DEN,QFC,VP,VS)
C                   Point XOLD is now shifted to the point of
C                   intersection.
C
                    GO TO 10
                  END IF
   20             CONTINUE
C
C                 Last contributions to the average:
                  DIST=XNEW(K1)-XOLD(K1)
                  A=A+DIST*(AOLD+ANEW)
                  B=B+DIST*(BOLD+BNEW)
                  C=C+DIST*(COLD+CNEW)
                  H=H+DIST
C
C                 Storing average elastic parameters in the memory:
                  IF(LA) THEN
                    RAM(IA)=RAM(IA)+0.5*H/A
                  END IF
                  IF(LB) THEN
                    RAM(IB)=RAM(IB)+0.5*H/B
                  END IF
                  IF(LC) THEN
                    RAM(IC)=RAM(IC)+0.5*H/C
                  END IF
C
C                 Screen output:
                  NOUT=NOUT+1
                  IF(100*NOUT.GT.IOUT*NALL) THEN
                    IOUT=IOUT+1
                    WRITE(*,'(A,I3)') '+',IOUT
                  END IF
                  IA=IA+J1
                  IB=IB+J1
                  IC=IC+J1
                END IF
C
C               Proceeding to the initial point of the next interval:
                XOLD(K1)=O1+D1*FLOAT(I1)+W1
                LA=.TRUE.
                LB=.TRUE.
                LC=.TRUE.
C
C               Evaluation of block indices and material parameters:
                CALL BLOCK(XOLD,0,0,ISRF,ISBOLD,ICBOLD)
                CALL ELPAR
     *               (ICBOLD,XOLD,LA,LB,LC,AOLD,BOLD,COLD,DEN,QFC,VP,VS)
C
C               Density and P-wave quality factor at the gridpoint:
                RAM(ID)=RAM(ID)+0.125*DEN
                RAM(IE)=RAM(IE)+0.125*QFC
                IF(FVP.NE.' ') THEN
                  RAM(IF)=RAM(IF)+0.125*VP
                END IF
                IF(FVS.NE.' ') THEN
                  RAM(IG)=RAM(IG)+0.125*VS
                END IF
C
                ID=ID+1
                IE=IE+1
                IF=IF+1
                IG=IG+1
   21         CONTINUE
              IA=IA-J1*NN1
              IB=IB-J1*NN1
              IC=IC-J1*NN1
              ID=ID-N1
              IE=IE-N1
              IF=IF-N1
              IG=IG-N1
   22       CONTINUE
   23     CONTINUE
          IA=IA+J2
          IB=IB+J2
          IC=IC+J2
          ID=ID+N1
          IE=IE+N1
          IF=IF+N1
          IG=IG+N1
   32   CONTINUE
        IA=IA+J3
        IB=IB+J3
        IC=IC+J3
   33 CONTINUE
C
C     Form of the output files:
      CALL RSEP3T('FORM',FORM,'formatted')
C
C     Writing the output:
      IA=1
      IB=IA+N1N2N3
      IC=IB+N1N2N3
      ID=IC+N1N2N3
      IE=ID+N1*N2*N3
      IF=IE+N1*N2*N3
      IG=IF+N1*N2*N3
      IF(FDEN.NE.' ') THEN
        CALL WARRAY(LU,FDEN,FORM,.TRUE.,0.,.FALSE.,0.,
     *                                              N1*N2*N3,RAM(ID))
      END IF
      IF(FQFC.NE.' ') THEN
        CALL WARRAY(LU,FQFC,FORM,.TRUE.,0.,.FALSE.,0.,
     *                                              N1*N2*N3,RAM(IE))
      END IF
      IF(FVP .NE.' ') THEN
        CALL WARRAY(LU,FVP ,FORM,.TRUE.,0.,.FALSE.,0.,
     *                                              N1*N2*N3,RAM(IF))
      END IF
      IF(FVS .NE.' ') THEN
        CALL WARRAY(LU,FVS ,FORM,.TRUE.,0.,.FALSE.,0.,
     *                                              N1*N2*N3,RAM(IG))
      END IF
      IF(NN1.GT.1) THEN
        IF(FA.NE.' ') THEN
          CALL WARRAY(LU,FA,FORM,.TRUE.,0.,.FALSE.,0.,N1N2N3,RAM(IA))
        END IF
        IF(FB.NE.' ') THEN
          CALL WARRAY(LU,FB,FORM,.TRUE.,0.,.FALSE.,0.,N1N2N3,RAM(IB))
        END IF
        IF(FC.NE.' ') THEN
          CALL WARRAY(LU,FC,FORM,.TRUE.,0.,.FALSE.,0.,N1N2N3,RAM(IC))
        END IF
      END IF
      WRITE(*,'(A,I3,A)') '+',IOUT,
     * ' per cent evaluated                                            '
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE ELPAR(ICB,X,LA,LB,LC,AINV,BINV,CINV,DEN,QFC,VP,VS)
      INTEGER ICB
      REAL X(3),AINV,BINV,CINV,DEN,QFC,VP,VS
      LOGICAL LA,LB,LC
C
C Subroutine to evaluate reciprocal values of elastic parameters
C A=lambda+2*mu, B=mu and C=lambda at a given point.
C
C Input:
C     ICB...  Index of the complex block.
C     X...    Coordinates of a given point.
C     LA,LB,LC... Defined logical variables.
C
C Output:
C     LA...   Unchanged input value if A is positive at the given point,
C             otherwise switched to .FALSE.
C     LB...   Unchanged input value if B is positive at the given point,
C             otherwise switched to .FALSE.
C     LC...   Unchanged input value if C is positive at the given point,
C             otherwise switched to .FALSE.
C     AINV... AINV=1/A if A=lambda+2*mu is positive, otherwise AINV=0.
C     BINV... BINV=1/B if B=mu          is positive, otherwise BINV=0.
C     CINV... CINV=1/C if C=lambda      is positive, otherwise CINV=0.
C     DEN...  Density at the given point.
C     QFC...  P-wave quality factor at the given point if it is finite.
C             Otherwise, QFC=0.
C     VP...   P-wave velocity at the given point.
C     VS...   S-wave velocity at the given point.
C
C Date: 1998, September 14
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:
      REAL UP(10),US(10),QP,QS,VD(10),QL
C
C.......................................................................
C
      IF(ICB.EQ.0) THEN
C       Free space:
        AINV=0.
        BINV=0.
        CINV=0.
        QL =0.
        DEN=0.
        VP =0.
        VS =0.
      ELSE
C       Material complex block:
        CALL PARM2(ICB,X,UP,US,DEN,QP,QS)
        CALL VELOC(ICB,UP,US,QP,QS,VP,VS,VD,QL)
        AINV=VP*VP*DEN
        BINV=VS*VS*DEN
        CINV=AINV-2.*BINV
      END IF
      IF(AINV.GT.0.) THEN
        AINV=1./AINV
      ELSE
        LA=.FALSE.
      END IF
      IF(BINV.GT.0.) THEN
        BINV=1./BINV
      ELSE
        LB=.FALSE.
      END IF
      IF(CINV.GT.0.) THEN
        CINV=1./CINV
      ELSE
        LC=.FALSE.
      END IF
      IF(QL.GT.0.) THEN
        QFC=1./QL
      ELSE
        QFC=0.
      END IF
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 'model.for'
C     model.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parm.for'
C     parm.for
      INCLUDE 'val.for'
C     val.for
      INCLUDE 'fit.for'
C     fit.for
      INCLUDE 'means.for'
C     means.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C