C 
C Subroutine file 'model.for' to specify a seismic model.
C
C Date: 1999, May 18
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following external procedures:
C     MODEL1..Subroutine designed to read the input data for the
C             description of the model and to store them in the memory.
C             Subroutine MODEL1 reads the input data (1)-(8) listed
C             below itself, then calls subroutine SRFC1 (included in the
C             subroutine file 'srfc.for') to read the input data (9) for
C             smooth surfaces, and finally calls subroutine PARM1
C             (included in the subroutine file 'parm.for') to read the
C             input data (10) for the parameters of the medium.
C             MODEL1
C     NSRFC...Integer function returning the number of surfaces covering
C             structural interfaces.
C             NSRFC
C     BLOCK...Subroutine designed to determine the mutual position of a
C             point and a simple and a complex block.  Just calls more
C             general subroutine BLOCKS which does the job.
C             BLOCK
C     BLOCKS..Subroutine designed to determine the mutual position of a
C             point and a simple and a complex block.  A more general
C             version of subroutine BLOCK, designed especially to work
C             with curves situated at structural interfaces.
C             BLOCKS
C     ISIDE...Auxiliary function to subroutine BLOCKS.
C             ISIDE
C     INTERF..Auxiliary subroutine to subroutine BLOCKS.
C             INTERF
C     SEPAR...Subroutine determining the surface separating two given
C             simple blocks.
C             SEPAR
C     VELOC...Subroutine transforming the values of a medium parameters
C             into velocities and loss factors.
C             VELOC
C     FPOWER..Subroutine evaluating the value and, possibly, the three
C             first and six second partial derivatives of a function,
C             if the value and the three first and six second partial
C             derivatives of the POWER-th power of the function are
C             known.  Particularly, this is an auxiliary subroutine
C             to the subroutine VELOC.
C             FPOWER
C
C Note:
C     The lines denoted by '*V' in the first two columns of file
C     'model.for' in subroutines VELOC and FPOWER are designed to
C     calculate the model variations with respect to the model
C     parameters.
C     File 'modelv.for', intended for the model inversion, is created
C     from 'model.for' by replacing each '*V' in the first two columns
C     by spaces using program 'clean.for'.  Subroutines VAR4 and VAR5
C     of file 'var.for' may then be called to handle the variations.
C
C.......................................................................
C
C                                                    
C Input data MODEL for the specification of the model:
C     The data are read in by the list directed input (free format).  In
C     the list of input data below, each numbered paragraph indicates
C     the beginning of a new input operation (new READ statement).  If
C     the first letter of the symbolic name of the input variable is
C     I-N, the corresponding value in input data must be of the type
C     INTEGER.  Otherwise (except TEXTM), the input parameter is of the
C     type REAL.
C (1) TEXTM
C     The string containing the name of the model.  Only the first 80
C     characters of the string are significant.
C (2) KOORS,NEXPV,NEXPQ,IVERT,/
C     KOORS...Specifies the type of the coordinate system:
C             KOORS.LE.0: Cartesian coordinates (default).
C             KOORS.EQ.1: polar spherical coordinates in radians,
C               (X1,X2,X3)=(colatitude,longitude,radius).
C             KOORS.GE.2: geographic spherical coordinates in radians,
C               (X1,X2,X3)=(longitude,latitude,radius).
C             If the coordinate system is right-handed (recommended),
C             all vectorial products are evaluated using the right-hand
C             rule, otherwise using the left-hand rule.
C             KOORS is passed to the subroutines of file 'metric.for'
C             by means of invocation of subroutine METR1 and presently
C             represents the only input data for the coordinate system.
C             Note that possible future additional data for the
C             coordinate system should be read by subroutine METR1 and
C             should be located between input data (2) and (3).
C             metric.for
C     NEXPV,NEXPQ... The default values are highly recommended!
C             Velocities powered to NEXPV and loss factors powered to
C             NEXPQ are reported by the subroutines evaluating isotropic
C             material parameters.
C             For example, unit values of NEXPV and NEXPQ indicate that
C             velocities and loss factors are the output parameters
C             of the subroutines evaluating isotropic material
C             parameters, indices equal -1 indicate reciprocal values of
C             these quantities, i.e.  slownesses and quality factors.
C             When using the basic submitted version of the subroutine
C             file 'parm.for', the default values of NEXPV=1, NEXPQ=1
C             are highly recommended.  Other values make sense only if a
C             user is submitting his own subroutines to evaluate
C             isotropic material parameters which, e.g., output the
C             slowness instead of the velocity.  In such a case,
C             switching NEXPV from 1 to -1 may avoid the modification
C             of the user's subroutines.
C     IVERT...Orientation of the vertical axis:
C             IVERT=0: unknown (default),
C             IVERT=+1: X1 vertical, pointing upwards,
C             IVERT=-1: X1 vertical, pointing downwards,
C             IVERT=+2: X2 vertical, pointing upwards,
C             IVERT=-2: X2 vertical, pointing downwards,
C             IVERT=+3: X3 vertical, pointing upwards,
C             IVERT=-3: X3 vertical, pointing downwards,
C             Has no influence on the calculations, and need not be
C             specified.  If it is non-zero, it may be considered for
C             plotting purposes.
C     /...    Obligatory slash at the end of line for future extensions.
C     Default: KOORS=0, NEXPV=1, NEXPQ=1, IVERT=0.
C (3) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX
C     Boundaries of the model.
C (4) NSRFC
C     Number of smooth surfaces in the model.  The surfaces are indexed
C     sequentially in any order by positive integers ISRFC from 1 to
C     NSRFC.
C     It is recommended to define only surfaces covering structural
C     interfaces (model surfaces) in this data set.  Auxiliary surfaces
C     related to particular source-receiver configurations, numerical
C     procedures, etc., should preferably be defined in other data sets.
C (5) NSB
C     Number of simple blocks in the model, defined in this data set.
C     The defined blocks are indexed sequentially by positive integers
C     ISB from 1 to NSB in the same order as they are specified in data
C     (6).  Intersecting simple blocks are allowed but not recommended.
C     All material simple blocks in the model must be defined.
C     Free-space blocks need not be defined in this data set but it is
C     recommended to define the free-space simple blocks in order to
C     speed up the determination of a free space during computations.
C     Defined free-space blocks also enhance the possibility to check
C     for the model consistency.  If the free-space simple blocks are
C     not, for some reason, defined here, they may be defined in the
C     additional data file designed just for the consistency check by
C     program modchk.for
C (6) NSB input operations (READ statements):
C     For each simple block with index ISB, the indices of the surfaces
C     forming the set F(+) and the indices of the surfaces forming the
C     set F(-).  The indices of surfaces from F(+) must be positive, the
C     indices of surfaces from F(-) must be indicated by negative signs.
C     The indices may be specified in an arbitrary order and must be
C     terminated by a slash.
C (7) NCB
C     Number of material complex blocks in the model.  The material
C     complex blocks are indexed sequentially by positive integers ICB
C     from 1 to NCB.  The free-space complex block is not indexed and
C     consists of all simple blocks not contained in the material
C     complex blocks.  Space covered by no simple block is also deemed
C     to be a free space.
C (8) NCB input operations (READ statements):
C     For each material complex block, the indices of material simple
C     blocks forming the complex block.  The indices may be specified
C     in an arbitrary order and must be terminated by a slash.
C     Each material simple block must appear exactly once within these
C     data lines.  Simple blocks defined by data (6) but not listed here
C     are deemed to be free-space simple blocks.
C (9) The data specifying NSRFC functions F(X1,X2,X3) describing the
C     smooth surfaces in the model. The data are read by subroutine
C     SRFC1.  For their description refer to subroutine SRFC1 (included
C     in the subroutine file 'srfc.for').
C     srfc.for: Input data
C (10) The data specifying the distribution of parameters of the model
C     in all NCB material complex blocks.  The data are read by
C     subroutine PARM1.  For their description refer to subroutine
C     PARM1 (included in the subroutine file 'parm.for').
C     parm.for: Input data
C For an example refer to the sample input data for the model.
C Example of data set MODEL
C
C.......................................................................
C
C Storage in the memory:
C     The input data (1) to (8) describing the structure of the model
C     are stored in common blocks /MODELT/ and /MODELC/ defined in the
C     include file 'model.inc'.
C     model.inc
C
C=======================================================================
C
C     
C
      SUBROUTINE MODEL1(LUN)
      INTEGER LUN
C
C Subroutine MODEL1 reads the input data (1)-(8) for the description
C of the model and stores them in common blocks /MODELT/ and /MODELC/.
C Then it calls subroutine SRFC1 (included in the subroutine file
C 'srfc.for') to read the input data (9) for smooth surfaces and to
C store them in the memory.  Finally, it calls subroutine PARM1
C (included in the subroutine file 'parm.for') to read the input data
C (10) for the parameters of the medium and to store them in the memory.
C
C Input:
C     LUN...  Logical unit number of the external input device
C             containing the input data.
C The input parameter is not altered.
C
C No output.
C
C Common blocks /MODELT/ and /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C All the storage locations of the common blocks are defined in this
C subroutine.
C
C Subroutines and external functions required:
      EXTERNAL METR1,SRFC1,PARM1
C     METR1...File 'metric.for'.
C     SRFC1 and subsequent routines... File 'srfc.for' and subsequent
C             files (like 'val.for' and 'fit.for').
C     PARM1 and subsequent routines... File 'parm.for' and subsequent
C             files (like 'val.for' and 'fit.for').
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,J,L
C
      READ(LUN,*) TEXTM
      I=0
      NEXPV=1
      NEXPQ=1
      IVERT=0
      READ(LUN,*) I,NEXPV,NEXPQ,IVERT
      CALL METR1(I)
      READ(LUN,*) BOUNDM
      READ(LUN,*) NSRFCS
C
C     Simple blocks:
C     Number of simple blocks
      READ(LUN,*) NSB
C     Initializing memory for indices of surfaces bounding simple blocks
      L=NSB+1
      DO 11 I=L,MSB
        KSB(I)=0
   11 CONTINUE
C     Reading indices of surfaces bounding simple blocks:
      DO 14 J=1,NSB
        READ(LUN,*) (KSB(I),I=L,MSB)
        DO 12 I=L,MSB
          IF(IABS(KSB(I)).GT.NSRFCS) THEN
C           311
            CALL ERROR('311 in MODEL: Block bounded by wrong interface')
C           Index of the surface bounding the simple block is greater
C           than the specified number of surfaces.
          ELSE IF(KSB(I).EQ.0) THEN
            KSB(J)=I-1
            L=I
            GO TO 13
          END IF
   12   CONTINUE
        GO TO 99
   13   CONTINUE
   14 CONTINUE
C
C     Complex blocks:
C     Number of complex blocks
      READ(LUN,*) NCB
C     Initializing memory for indices of simple blocks forming c. blocks
      L=NCB+1
      DO 21 I=L,MCB
        KCB(I)=0
   21 CONTINUE
C     Reading indices of simple blocks forming complex blocks
      DO 24 J=1,NCB
        READ(LUN,*) (KCB(I),I=L,MCB)
        DO 22 I=L,MCB
          IF(KCB(I).LT.0.OR.KCB(I).GT.NSB) THEN
C           312
            CALL ERROR
     *            ('312 in MODEL: C. block composed of wrong s. block.')
C           Complex block composed of wrong simple block:
C           Index of a simple block composing the complex block is
C           greater than the specified number of simple blocks.
          ELSE IF(KCB(I).EQ.0) THEN
            KCB(J)=I-1
            L=I
            GO TO 23
          END IF
   22   CONTINUE
        GO TO 99
   23   CONTINUE
   24 CONTINUE
C
C     Smooth surfaces:
      CALL SRFC1(LUN,NSRFCS)
C
C     Material parameters:
      CALL PARM1(LUN,NCB)
      RETURN
C
   99 CONTINUE
C       310
        CALL ERROR('310 in MODEL1: Insufficient memory in /MODELC/')
C       Insufficient memory for the input data in common block /MODELC/.
C       The dimensions MSB or MCB of arrays KSB or KCB, respectively,
C       must be enlarged.
C       Refer to include file model.inc
      END
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION NSRFC()
C
C Integer function NSRFC is designed to return the number of surfaces
C covering structural interfaces.
C
C No input.
C
C Output:
C     NSRFC...Number of surfaces covering structural interfaces.
C
C Common block /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1989, December 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     No auxiliary storage locations.
C
      NSRFC=NSRFCS
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE BLOCK(COOR,ISRF1,ISB1,ISRF2,ISB2,ICB2)
      REAL COOR(3)
      INTEGER ISRF1,ISB1,ISB2,ICB2,ISRF2
C
C This subroutine searches for the simple block and the complex block in
C which a given point is situated.  This routine may be also used to
C determine the index of a block touching a specified block at a given
C point situated on the boundary of the specified block (the situation
C which may occur when a ray impinges on a boundary of a block).
C Another function of the routine is to determine the index of the
C surface bounding a given simple block and separating it from the given
C point.
C
C Input:
C     COOR... Array containing coordinates X1, X2, X3 of a given point.
C     ISRF1...Controls the possible determination of the simple block
C             touching block ISB1 at the given point.
C             ISRF1.EQ.0:  The simple block touching ISB1 at the given
C               point need not be determined.  This is the most frequent
C               option, used, e.g., if the given point is not situated
C               at a structural interface or if the simple block
C               touching ISB1 is already known.
C             ISRF1.NE.0:  Index of the surface at which the given point
C               is situated and which separates simple block ISB1 from
C               another block to be determined and returned in ISB2.
C               Since the given point may be situated on either side of
C               the surface, surface IABS(ISRF1) is skipped in the list
C               of surfaces limiting simple blocks when determining the
C               position of the given point with respect to simple
C               blocks.  The sign of ISRF1 is ignored.
C               This option cannot be used together with ISB1.EQ.0.
C     ISB1... Index of the simple block which position with respect to
C             the given point will be checked.
C             ISB1.EQ.0:  No simple block given.  The simple block in
C               which the given point is situated will be determined
C               and returned in ISB2.  ISRF2 will be set to zero.
C               The determination starts with initial estimate ISBOLD,
C               continuing with ISBOLD+1, ISBOLD-1, ISBOLD+2, ISBOLD-2,
C               etc.  Here ISBOLD is the simple block returned in ISB2
C               during the last invocation of this subroutine with
C               ISB1.EQ.0 on input and ISB2.NE.0 on output.  ISBOLD=1
C               before such an invocation.
C             ISB1.GT.0:  The position of the given point with
C               respect to simple block ISB1 will be examined.
C               If the point is situated in ISB1, the output value
C                 of ISRF2 is set to 0.
C               If the point is situated in ISB1 and ISRF1.EQ.0,
C                 the output value of ISB2 is set to ISB1.
C               If the point is situated in ISB1 and ISRF1.NE.0,
C                 the simple block containing the given point,
C                 separated from block ISB1 by surface ISRF1 and not
C                 separated from block ISB1 by another surface will be
C                 determined and returned in ISB2.
C               If the point is not situated in ISB1, the surface
C                 separating the point from ISB1 is output in ISRF2.
C                 In addition, the simple block containing the given
C                 point, separated from block ISB1 by surface ISRF2 and
C                 not separated from block ISB1 by another surface will
C                 be determined and returned in ISB2.
C               The determination starts with initial estimate ISB1,
C               continuing with ISB1+1, ISB1-1, ISB1+2, ISB1-2, etc.
C None of the input parameters are altered.
C
C Output:
C     ISRF2...For the given point not situated inside block ISB1 or
C             at its boundary, ISRF2 has the meaning of the index of
C             one of the surfaces bounding simple block ISB1 and
C             separating the given point from simple block ISB1,
C             supplemented by a sign '+' or '-' for the simple block
C             ISB1 situated at the positive or negative side of the
C             surface, respectively.  If the point is separated from
C             simple block ISB1 by several surfaces, the first one of
C             the surfaces, for which ISB2 (see description of ISB2
C             below) will be positive, is preferred.  Roughly speaking,
C             surface ISRF2 at which simple block ISB1 touches simple
C             block containing the given point is preferred.
C             ISRF2=0 if ISB1=0 or if the given point is situated inside
C             block ISB1 or at its boundary.
C             If ISRF1.NE.0, the position with respect to surface
C             IABS(ISRF1) is not checked.
C     ISB2... Index of the simple block containing the given point.
C             If ISB1.EQ.0:
C               ISB2 is the index of the simple block containing the
C               given point.
C               The first simple block of ISBOLD, ISBOLD+1, ISBOLD-1,
C               ISBOLD+2, ISBOLD-2, etc., containing the given point,
C               is returned in ISB2.  Here ISBOLD is the simple block
C               returned in ISB2 during the previous invocation of this
C               subroutine with ISB1.EQ.0.  ISBOLD=1 during such the
C               first invocation.
C             If ISB1.NE.0:
C               If ISRF1.EQ.0 and ISRF2.EQ.0, ISB2 must not be
C                 separated from simple block ISB1 by any surface.
C                 A surface "separates" two simple blocks if both blocks
C                 are bounded by the surface and are situated at its
C                 opposite sides.
C               If ISRF1.NE.0 and ISRF2.EQ.0, ISB2 must be separated
C                 from simple block ISB1 by surface IABS(ISRF1) and
C                 must not be separated by another surface.
C               If ISRF2.NE.0, ISB2 must be separated from simple block
C                 ISB1 by surface ISRF2 and must not be separated by
C                 another surface.
C               The first simple block of ISB1, ISB1+1, ISB1-1,
C               ISB1+2, ISB1-2, etc., having the above properties, is
C               returned in ISB2.
C             If there is no simple block of the above properties,
C             ISB2=0.
C     ICB2... Index of the complex block in which simple block ISB2 is
C             situated.  ICB2=0 if ISB2=0 or if simple block ISB2 is not
C             situated in any material complex block.
C
C Examples:
C     For each point of a set of discrete points, we wish to determine
C             simple block ISB and complex block ICB in which the point
C             is situated.
C     (a)     We thus use
C               CALL BLOCK(COOR,0,0,I,ISB,ICB)
C     When tracing a curve (e.g., a ray), we wish rather to find each
C             crossed interface than only to determine the simple blocks
C             at discrete points along the curve.
C     (b)     If the previous point along the curve is situated in
C             simple block ISB, we use
C             CALL BLOCK(COOR,0,ISB,ISRF1,ISBNEW,ICBNEW)
C             and get ISRF1.EQ.0 if no surface limiting simple block ISB
C             has been crossed.
C     (c)     If ISRF1.NE.0, IABS(ISRF1) is the surface which has been
C             crossed.  We should determine the point of intersection of
C             the curve with surface IABS(ISRF1).  Then we use
C             CALL BLOCK(COOR,ISRF1,ISB,ISRF2,ISB2,ICB2)
C             to find simple block ISB2 and complex block ICB2 at the
C             other side of surface ISRF1.  If ISRF2.EQ.0, the blocks
C             are successfully determined.  If ISRF2.NE.0, surface
C             IABS(ISRF2) has also been crossed and we have to repeat
C             step (c) with ISRF1=ISRF2.
C     When tracing a curve along surface ISRFA, we call subroutine
C             BLOCKS instead of BLOCK, with analogous arguments:
C     (b')    ISRF(1)=0
C             ISRF(2)=ISRFA
C             CALL BLOCKS(COOR,2,ISRF,ISB,ISRF1,ISBNEW,ICBNEW)
C     (c')    ISRF(1)=ISRF1
C             ISRF(2)=ISRFA
C             CALL BLOCKS(COOR,2,ISRF,ISB,ISRF2,ISB2,ICB2)
C     Analogously, when tracing the curve of intersection of surfaces
C             ISRFA and ISRFB, we call subroutine BLOCKS in the
C             following way:
C     (b")    ISRF(1)=0
C             ISRF(2)=ISRFA
C             ISRF(3)=ISRFB
C             CALL BLOCKS(COOR,3,ISRF,ISB,ISRF1,ISBNEW,ICBNEW)
C     (c")    ISRF(1)=ISRF1
C             ISRF(2)=ISRFA
C             ISRF(3)=ISRFB
C             CALL BLOCKS(COOR,3,ISRF,ISB,ISRF2,ISB2,ICB2)
C
C Subroutines and external functions required:
      EXTERNAL BLOCKS
C     BLOCKS... This file.
C
C Date: 1999, March 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      INTEGER ISRF(1)
C
      ISRF(1)=ISRF1
      CALL BLOCKS(COOR,1,ISRF,ISB1,ISRF2,ISB2,ICB2)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE BLOCKS(COOR,NSRF1,ISRF,ISB1,ISRF2,ISB2,ICB2)
      REAL COOR(3)
      INTEGER NSRF1,ISRF(NSRF1),ISB1,ISRF2,ISB2,ICB2
C
C This subroutine is a generalization of subroutine BLOCK to follow
C points along structural interfaces and edges.  Subroutine BLOCKS
C searches for the simple block and the complex block in which a given
C point is situated.  The given point may be situated in the vicinity
C of one or more structural interfaces.  In such case, blocks at the
C proper side of the interface(s) will be searched for.  This routine
C may also be used to determine the index of a block touching a
C specified block at a given point situated on the boundary of the
C specified block (the situation which may occur when a ray impinges on
C a boundary of a block).  Another function of the routine is to
C determine the index of the surface bounding a given simple block and
C separating it from the given point.
C
C Input:
C     COOR... Array containing coordinates X1, X2, X3 of a given point.
C     NSRF1...Number of surface indices specified in array ISRF.
C             There must be NSRF1.GE.1.  For ISB1.EQ.0., there must be
C             NSRF1.EQ.1 and ISRF(1).EQ.0.
C     ISRF... Array containing the indices of the surfaces at which the
C             given point is situated.  Since the given point may be
C             situated on either side of the surfaces, the surfaces
C             listed in array ISRF are skipped in the list of surfaces
C             limiting simple blocks when determining the position of
C             the given point with respect to simple blocks.
C             The signs of the indices in array ISRF are ignored.
C             ISRF(1):  Controls the possible determination of the
C               simple block touching block ISB1 at the given point.
C               ISRF(1).EQ.0:  The simple block touching ISB1 at the
C                 given point need not be determined.  This is the most
C                 frequent option, used, e.g., if the given point is not
C                 situated at a structural interface or if the simple
C                 block touching ISB1 is already known.
C               ISRF(1).NE.0:  Index of the surface at which the given
C                 point is situated and which separates simple block
C                 ISB1 from another block to be determined.  Refer to
C                 the description of output parameter ISB2.
C                 This option cannot be used together with ISB1.EQ.0.
C             ISRF(2) to ISRF(NSRF1):  Indices of the surfaces at
C               which the given point is situated, other than ISRF(1).
C               These indices are designed to specify the surface along
C               which a curve is traced, or the surfaces which curve of
C               intersection is traced.
C     ISB1... Index of the simple block which position with respect to
C             the given point will be checked.
C             ISB1.EQ.0:  No simple block given.  The simple block in
C               which the given point is situated will be determined
C               and returned in ISB2.  ISRF2 will be set to zero.
C               The determination starts with initial estimate ISBOLD,
C               continuing with ISBOLD+1, ISBOLD-1, ISBOLD+2, ISBOLD-2,
C               etc.  Here ISBOLD is the simple block returned in ISB2
C               during the last invocation of this subroutine with
C               ISB1.EQ.0 and ISB2.NE.0.  ISBOLD=1 before such an
C               invocation.
C             ISB1.GT.0:  The position of the given point with
C               respect to simple block ISB1 will be examined.
C               If the point is situated in ISB1, the output value
C                 of ISRF2 is set to 0.
C               If the point is situated in ISB1 and ISRF(1).EQ.0,
C                 the output value of ISB2 is set to ISB1.
C               If the point is situated in ISB1 and ISRF(1).NE.0,
C                 the simple block containing the given point,
C                 separated from block ISB1 by surface ISRF(1) and not
C                 separated from block ISB1 by another surface will be
C                 determined and returned in ISB2.
C               If the point is not situated in ISB1, the surface
C                 separating the point from ISB1 is output in ISRF2.
C                 In addition, the simple block containing the given
C                 point, separated from block ISB1 by surface ISRF2 and
C                 not separated from block ISB1 by another surface will
C                 be determined and returned in ISB2.
C               The determination starts with initial estimate ISB1,
C               continuing with ISB1+1, ISB1-1, ISB1+2, ISB1-2, etc.
C None of the input parameters are altered.
C
C Output:
C     ISRF2...For the given point not situated inside block ISB1 or
C             at its boundary, ISRF2 has the meaning of the index of
C             one of the surfaces bounding simple block ISB1 and
C             separating the given point from simple block ISB1,
C             supplemented by a sign '+' or '-' for the simple block
C             ISB1 situated at the positive or negative side of the
C             surface, respectively.  If the point is separated from
C             simple block ISB1 by several surfaces, the first one with
C             positive ISB2 (see below) is preferred.
C             ISRF2=0 if ISB1=0 or if the given point is situated inside
C             block ISB1 or at its boundary.
C             Note that surfaces ISRF(1) to ISRF(NSRF1) are skipped
C             when determining the position of the given point with
C             respect to simple blocks.
C     ISB2... Index of the simple block containing the given point.
C             If ISB1.EQ.0:
C               ISB2 is the index of the simple block containing the
C               given point.
C               The first simple block of ISBOLD, ISBOLD+1, ISBOLD-1,
C               ISBOLD+2, ISBOLD-2, etc., containing the given point,
C               is returned in ISB2.  Here ISBOLD is the simple block
C               returned in ISB2 during the previous invocation of this
C               subroutine with ISB1.EQ.0.  ISBOLD=1 during such the
C               first invocation.
C             If ISB1.NE.0:
C               If ISRF(1).EQ.0 and ISRF2.EQ.0, ISB2 must not be
C                 separated from simple block ISB1 by any surface.
C                 A surface "separates" two simple blocks if both blocks
C                 are bounded by the surface and are situated at its
C                 opposite sides.
C               If ISRF(1).NE.0 and ISRF2.EQ.0, ISB2 must be separated
C                 from simple block ISB1 by surface IABS(ISRF(1)) and
C                 must not be separated by another surface.
C               If ISRF2.NE.0, ISB2 must be separated from simple block
C                 ISB1 by surface ISRF2 and must not be separated by
C                 another surface.
C               The first simple block of ISB1, ISB1+1, ISB1-1,
C               ISB1+2, ISB1-2, etc., having the above properties, is
C               returned in ISB2.
C             If there is no simple block of the above properties,
C             ISB2=0.
C     ICB2... Index of the complex block in which simple block ISB2 is
C             situated.  ICB2=0 if ISB2=0 or if simple block ISB2 is not
C             situated in any material complex block.
C
C Common block /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL ISIDE,INTERF,SRFC2
      INTEGER ISIDE
C     ISIDE,INTERF... This file.
C     SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C             files (like 'val.for' and 'fit.for').
C
C Date: 1999, March 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER MSEPAR
      PARAMETER (MSEPAR=20)
      INTEGER NSEPAR,KSEPAR(MSEPAR),JSEPAR,ISEPAR,ISIDE1,I,J,K,II(1)
      INTEGER ISBOLD,ISB0,ISRF1
      SAVE ISBOLD
      DATA ISBOLD/1/
C
C.......................................................................
C
      ISRF1=ISRF(1)
C
C     Checking input values:
      IF(ISRF1.LT.-NSRFCS.OR.ISRF1.GT.NSRFCS) THEN
C       313
        CALL ERROR('313 in BLOCK: Wrong index of surface')
C       Absolute value of the input parameter ISRFC1 (index of the
C       surface) is greater than the number NSRFC of the surfaces
C       covering structural interfaces.
      END IF
      IF(ISB1.LT.0.OR.ISB1.GT.NSB) THEN
C       314
        CALL ERROR('314 in BLOCK: Wrong index of simple block')
C       Parameter ISB1 (index of the simple block) is either
C       negative or greater than the number NSB of simple blocks.
      END IF
      IF(ISB1.EQ.0) THEN
        IF(NSRF1.NE.1.OR.ISRF1.NE.0) THEN
C         315
          CALL ERROR('315 in BLOCK: No simple block specified')
C         If no simple block ISB1 is specified, NSRF1 must be 1 and
C         ISRF1=ISRF(1) must be 0.
        END IF
      END IF
C
C     Initial simple block ISB0:
      IF(ISB1.EQ.0) THEN
        ISB0=ISBOLD
      ELSE
        ISB0=ISB1
      END IF
C
C     Position of the given point with respect to simple block ISB0:
      ISRF2=0
      CALL INTERF(COOR,NSRF1,ISRF,ISB0,MSEPAR,NSEPAR,KSEPAR)
      IF(NSEPAR.EQ.0) THEN
C       The point is inside simple block ISB0:
        IF(ISRF1.EQ.0) THEN
          ISB2=ISB0
          GO TO 90
        ELSE
          NSEPAR=1
          KSEPAR(1)=ISRF1
        END IF
      ELSE
        IF(ISB1.EQ.0) THEN
          NSEPAR=1
        ELSE
          ISRF2=KSEPAR(1)
        END IF
      END IF
C
C     Search for the simple block in which the given point is situated:
      DO 20 JSEPAR=1,NSEPAR
        IF(ISB1.NE.0) THEN
          ISEPAR=IABS(KSEPAR(JSEPAR))
          ISIDE1=-ISIDE(ISEPAR,ISB1)
        END IF
C       Loop over simple blocks
        DO 19 J=1,MAX0(ISB0-1,NSB-ISB0)
          DO 18 ISB2=ISB0+J,ISB0-J,-2*J
            IF(ISB2.GT.0.AND.ISB2.LE.NSB) THEN
C             Selecting simple block ISB2 according to ISB1
              IF(ISB1.NE.0) THEN
C               Loop for surfaces bounding block ISB1
                DO 11 I=KSB(ISB1-1)+1,KSB(ISB1)
C                 Skipping simple blocks separated from ISB1 by another
C                 surface than ISEPAR
                  K=KSB(I)
                  IF(IABS(K).NE.ISEPAR) THEN
                    IF(ISIDE(K,ISB1).EQ.-ISIDE(K,ISB2)) THEN
                      GO TO 17
                    END IF
                  END IF
   11           CONTINUE
C               Skipping simple blocks not separated from ISB1 by ISEPAR
                IF(ISIDE(ISEPAR,ISB2).NE.ISIDE1) THEN
                  GO TO 17
                END IF
              END IF
C             Determining the position of the given point with respect
C             to the simple block
              CALL INTERF(COOR,NSRF1,ISRF,ISB2,1,I,II)
              IF(I.EQ.0) THEN
                IF(ISB1.EQ.0) THEN
                  ISBOLD=ISB2
                ELSE
                  ISRF2=KSEPAR(JSEPAR)
                END IF
                GO TO 90
              END IF
   17         CONTINUE
            END IF
   18     CONTINUE
   19   CONTINUE
   20 CONTINUE
C     No simple block has been found:
      ISB2=0
      ICB2=0
      RETURN
C
C     Determination of the complex block:
   90 CONTINUE
      DO 92 J=1,NCB
        DO 91 I=KCB(J-1)+1,KCB(J)
          IF(KCB(I).EQ.ISB2) THEN
            ICB2=J
            RETURN
          END IF
   91   CONTINUE
   92 CONTINUE
C     No complex block:
      ICB2=0
      RETURN
      END
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION ISIDE(ISRF,ISB)
      INTEGER ISRF,ISB
C
C This is an auxiliary function to the subroutine BLOCKS.
C This function determines the mutual position of a surface and a simple
C block.
C
C Input:
C     ISRF... Index of the surface.  The sign is ignored.
C     ISB...  Index of the simple block.
C None of the input parameters are altered.
C
C Output:
C     ISIDE...ISIDE=-1: The simple block is bounded by the surface and
C               is situated on its negative side.
C             ISIDE= 1: The simple block is bounded by the surface and
C               is situated on its positive side.
C             ISIDE= 2: The simple block is not bounded by the surface.
C
C Common block /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1989, December 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER IS,LS,MS
C
      LS=KSB(ISB-1)+1
      MS=KSB(ISB)
C
C     Loop for surfaces bounding simple block ISB:
      DO 1 IS=LS,MS
        IF(IABS(KSB(IS)).EQ.IABS(ISRF)) THEN
          ISIDE=ISIGN(1,KSB(IS))
          RETURN
        END IF
    1 CONTINUE
C
      ISIDE=2
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE INTERF(COOR,NSRF1,ISRF1,ISB,MSRF2,NSRF2,ISRF2)
      REAL COOR(3)
      INTEGER NSRF1,ISRF1(NSRF1),ISB,MSRF2,NSRF2,ISRF2(MSRF2)
C
C This is an auxiliary subroutine to the subroutine BLOCKS.
C This subroutine determines the position of a given point with respect
C to a given simple block.
C
C Input:
C     COOR... Array containing coordinates X1, X2, X3 of a given point.
C     NSRF1...Number of surface indices specified in array ISRF1.
C             Must be NSRF1.GE.1.
C     ISRF1...Array containing the indices of the surfaces at which the
C             given point is situated.  Since the given point may be
C             situated on either side of the surfaces, the surfaces
C             listed in array ISRF1 are skipped in the list of surfaces
C             limiting the simple block when determining the position of
C             the given point with respect to the simple block.
C             The signs of the indices in array ISRF1 are ignored.
C             Zero indices are allowed and have no neaning.
C     ISB...  Index of the given simple block.
C     MSRF2...Dimension of array ISRF2.
C None of the input parameters are altered.
C
C Output:
C     NSRF2...Number of surfaces separating the given point from simple
C             block ISB.
C             NSRF2=0  if the given point is situated inside simple
C             block ISB.
C     ISRF2...Indices of surfaces separating the given point from simple
C             block ISB, supplemented by sign '+' or '-' for simple
C             block ISB situated at the positive or negative side of the
C             surface, respectively.
C             The first MSRF2 such surfaces from the list of surfaces
C             limiting the simple block are reported.
C
C Common block /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL SRFC2
C     SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C              files (like 'val.for' and 'fit.for').
C
C Date: 1999, March 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER IS,JS,KS,I
      REAL F(10)
C     F...    Auxiliary array to contain the value and partial
C             derivatives F, F1, F3, F11, F12, F22, F13, F23, F33 of the
C             function describing surfaces at the given point.
C
C     Speed-up storage locations:
      INTEGER MOLD
      PARAMETER (MOLD=20)
      INTEGER NOLD,IOLD(MOLD)
      REAL COLD(3)
      SAVE COLD,NOLD,IOLD
      DATA COLD/3*-999999./
C
C.......................................................................
C
C     Already examined surfaces:
      IF(COOR(1).NE.COLD(1).OR.
     *   COOR(2).NE.COLD(2).OR.COOR(3).NE.COLD(3)) THEN
        NOLD=0
      END IF
C
C     Loop for surfaces bounding simple block ISB:
      NSRF2=0
      DO 9 IS=KSB(ISB-1)+1,KSB(ISB)
        KS=KSB(IS)
        JS=IABS(KS)
C       Skipping surfaces of array ISRF1
        DO 1 I=1,NSRF1
          IF(JS.EQ.IABS(ISRF1(I))) THEN
            GO TO 8
          END IF
    1   CONTINUE
C       Surfaces already examined
        DO 2 I=1,NOLD
          IF(JS.EQ.IABS(IOLD(I))) THEN
            IF(IOLD(I)*KS.LT.0) THEN
              NSRF2=NSRF2+1
              ISRF2(NSRF2)=KS
              IF(NSRF2.GE.MSRF2) THEN
                RETURN
              END IF
            END IF
            GO TO 8
          END IF
    2   CONTINUE
C       The surface has to be examined
        CALL SRFC2(JS,COOR,F)
        IF(NOLD.LT.MOLD) THEN
          IF(F(1).LT.0.) THEN
            NOLD=NOLD+1
            IOLD(NOLD)=-JS
          ELSE IF(F(1).GT.0.) THEN
            NOLD=NOLD+1
            IOLD(NOLD)= JS
          END IF
        END IF
        IF(F(1)*FLOAT(KS).LT.0.) THEN
          NSRF2=NSRF2+1
          ISRF2(NSRF2)=KS
          IF(NSRF2.GE.MSRF2) THEN
            RETURN
          END IF
        END IF
    8   CONTINUE
    9 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SEPAR(ISB1,ISB2,NSRF,ISRF)
      INTEGER ISB1,ISB2,NSRF,ISRF
C
C Subroutine determining the surface separating given simple blocks,
C i.e., the surface limiting the simple blocks, which are situated one
C at the positive side of the surface and the other at the negative side
C of the surface.
C
C Input:
C     ISB1,ISB2... Indices of given simple blocks.
C None of the input parameters are altered.
C
C Output:
C     NSRF... Number of surfaces separating the simple blocks.
C     ISRF... Index of a surface separating the simple blocks,
C             supplemented by sign minus if simple block ISB1 is
C             situated at its negative side.
C             If NSRF.EQ.0:  ISRF contains unchanged input value.
C             If NSRF.GT.1:  ISRF is one of the surfaces separating the
C               simple blocks.
C
C Common block /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1999, May 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I1,I2,KSBI2
C
      NSRF=0
C
C     Loop for surfaces bounding simple block ISB2
      DO 2 I2=KSB(ISB2-1)+1,KSB(ISB2)
        KSBI2=-KSB(I2)
C       Loop for surfaces bounding simple block ISB1
        DO 1 I1=KSB(ISB1-1)+1,KSB(ISB1)
          IF(KSB(I1).EQ.KSBI2) THEN
            NSRF=NSRF+1
            IF(NSRF.EQ.1) THEN
              ISRF=KSB(I1)
            END IF
          END IF
    1   CONTINUE
    2 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE VELOC(IWAVE,UP,US,QP,QS,VP,VS,VD,QL)
      INTEGER IWAVE
      REAL UP(10),US(10),QP,QS,VP,VS,VD(10),QL
C
C This subroutine transforms the values of parameters of the medium into
C velocities and loss factors.
C
C Input:
C     IWAVE...Type of wave.
C             IWAVE.GE.0: P wave,
C             IWAVE.LT.0: S wave.
C     UP,US...Powers of P and S wave velocities and their first and
C             second partial derivatives (the exponent of the powers is
C             NEXPV, see 'Input data for the model'), in order U, U1,
C             U2, U3, U11, U12, U22, U13, U23, U33.
C     QP,QS...Powers of the loss factors of P and S waves (the exponent
C             of the powers is NEXPQ, see 'Input data for the model').
C None of the input parameters are altered.
C
C Output:
C     VP,VS...P and S wave velocities.
C     VD...   Velocity and its first and second partial derivatives
C             ordered as UP, US, corresponding to the wave specified by
C             IWAVE, in order V, V1, V2, V3, V11, V12, V22, V13, V23,
C             V33.
C     QL...   Loss factor corresponding to the wave specified by IWAVE.
C
C Common block /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL FPOWER
C     FPOWER...This file.
C
C Date: 1992, December 31
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage location:
      REAL POWER,AUX1(1),AUX2(1)
C
      POWER=FLOAT(NEXPV)
      IF(IWAVE.GE.0) THEN
        CALL FPOWER(10,UP,POWER,VD)
        CALL VAR5(1,1)
        VP=VD(1)
        CALL FPOWER(1,US,POWER,AUX2)
        VS=AUX2(1)
        AUX1(1)=QP
      ELSE
        CALL FPOWER(1,UP,POWER,AUX2)
        VP=AUX2(1)
        CALL FPOWER(10,US,POWER,VD)
        CALL VAR5(2,2)
        VS=VD(1)
        AUX1(1)=QS
      END IF
      CALL FPOWER(1,AUX1,FLOAT(NEXPQ),AUX2)
      QL=AUX2(1)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE FPOWER(N,FINP,POWER,FOUT)
      INTEGER N
      REAL FINP(N),POWER,FOUT(N)
C
C This subroutine evaluates the value and, possibly, the three first and
C six second partial derivatives of a function if the value and the
C three first and six second partial derivatives of the POWER-th power
C of the function are known.
C
C Input:
C     N...    For N=1: only the function value is evaluated. The
C               derivatives are ignored.
C             For N=4: the value and the three first partial derivatives
C               are evaluated.
C             For N=10: the value and the three first and six second
C               partial derivatives are evaluated.
C     FINP... Array containing the value, the first and second partial
C             derivatives of the POWER-th power of the function to be
C             evaluated, in the order F, F1, F2, F3, F11, F12, F22, F13,
C             F23, F33.  For N=1, only the function value is required.
C     POWER...The specified function is equal to the POWER-th power of
C             the corresponding physical quantity.
C             POWER=0: Zero output array FOUT is generated.
C None of the input parameters are altered (except FINP if this
C parameter and FOUT are identical in the calling sequence).
C
C Output:
C     FOUT... Array containing the value, the first and second partial
C             derivatives of the evaluated function, in the order F, F1,
C             F2, F3, F11, F12, F22, F13, F23, F33.  This parameter may
C             coincide with FINP, in which case FINP is destroyed on
C             output.  Note that this coincidence is an exception to
C             ANSI standard of FORTRAN 77.
C
C No subroutines and external functions required.
C
C Date: 1999, January 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL F,AUX1,AUX2
C
      IF(POWER.EQ.0.) THEN
        DO 1 I=1,N
          FOUT(I)=0.
    1   CONTINUE
      ELSE IF(0.999.LT.POWER.AND.POWER.LT.1.001) THEN
        DO 2 I=1,N
          FOUT(I)=FINP(I)
    2   CONTINUE
        CALL VAR4(0,1.)
      ELSE IF(FINP(1).LT.0.) THEN
C       317
        CALL ERROR('317 in FPOWER: Negative material parameter')
C       Nonunit power of a material parameter is not allowed to be
C       negative.  The negative value may be caused by oscillatory
C       character of interpolated positive values.
      ELSE IF(FINP(1).EQ.0.) THEN
        IF(POWER.LT.0.) THEN
C         318
          CALL ERROR('318 in FPOWER: Zero inverse material parameter')
C         Negative power of a material parameter cannot be zero.
        ELSE
          FOUT(1)=0.
          DO 3 I=2,N
            IF(FINP(I).NE.0.) THEN
C             319
              CALL ERROR('319 in FPOWER: Zero material parameter')
C             Nonunit power of zero material parameter is not allowed to
C             have nonzero derivatives.
            END IF
            FOUT(I)=0.
    3     CONTINUE
        END IF
      ELSE
        IF(-1.001.LT.POWER.AND.POWER.LT.-0.999) THEN
          F=1./FINP(1)
        ELSE
          F=FINP(1)**(1./POWER)
        END IF
        FOUT(1)=F
        IF(N.GT.1) THEN
          AUX1= F/(FINP(1)*POWER)
          AUX2= (POWER-1.)/F
          FOUT(2)=AUX1*FINP(2)
          FOUT(3)=AUX1*FINP(3)
          FOUT(4)=AUX1*FINP(4)
          IF(N.GT.4) THEN
            FOUT(5)=AUX1*FINP(5)-AUX2*FOUT(2)*FOUT(2)
            FOUT(6)=AUX1*FINP(6)-AUX2*FOUT(2)*FOUT(3)
            FOUT(7)=AUX1*FINP(7)-AUX2*FOUT(3)*FOUT(3)
            FOUT(8)=AUX1*FINP(8)-AUX2*FOUT(2)*FOUT(4)
            FOUT(9)=AUX1*FINP(9)-AUX2*FOUT(3)*FOUT(4)
            FOUT(10)=AUX1*FINP(10)-AUX2*FOUT(4)*FOUT(4)
          END IF
          CALL VAR4(0,AUX1)
          CALL VAR4(2,-AUX2*FOUT(2))
          CALL VAR4(3,-AUX2*FOUT(3))
          CALL VAR4(4,-AUX2*FOUT(4))
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C