C 
C Subroutine file 'model.for' to specify a seismic model.
C
C Date: 1998, October 26
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.
C             BLOCK
C     ISIDE...Auxiliary function to subroutine BLOCK.
C             ISIDE
C     INTERF..Auxiliary subroutine to subroutine BLOCK.
C             INTERF
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, if
C             the value and the three first and six second partial
C             derivatives of the POWER-th power of the function are
C             are 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 (and usually are not) defined in
C     this data set.  Free-space blocks may be defined at user's
C     discretion in order to enhance the possibility to check for the
C     model consistency.  Free-space blocks may be defined either here
C     or in the additional data file designed just for the consistency
C     check by program MODCHK.
C     Program MODCHK
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 blocks are
C     indexed sequentially by positive integers ICB from 1 to NCB.  The
C     free-space blocks are not indexed.
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(INSIDE,COOR,ISRF1,ISB1,ISRF2,ISB2,ICB2,F)
      LOGICAL INSIDE
      REAL COOR(3)
      INTEGER ISRF1,ISB1,ISB2,ICB2,ISRF2
      REAL F(10)
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 any of
C the surfaces bounding a block and separating the block from the given
C point.
C
C Input:
C     INSIDE..INSIDE=.TRUE.:  the given point is checked for its
C               location inside simple block ISB1.
C             INSIDE=.FALSE.:  the given point is expected to be
C               situated outside simple block ISB1, and simple block
C               ISB2 (different from ISB1) in which the point is
C               situated is being searched for.  If ISRF1.NE.0, the
C               position with respect to surface ISRF1 is not checked,
C               and the simple blocks situated at the other side of
C               ISRF1 than simple block ISB1 are preferred when
C               searching for ISB2.
C     COOR... Array containing coordinates X1, X2, X3 of a given point.
C     ISRF1...ISRF1.NE.0: index of the surface at which the given point
C               is situated.  The sign is ignored.
C             ISRF1.EQ.0: the given point is not situated at any surface
C     ISB1... For ISRF1.NE.0: index of a simple block bounded by the
C               surface IABS(ISRF1) - search for the index of a
C               neighbouring simple block touching ISB1 at the given
C               point.
C             For ISRF1.EQ.0: index of an arbitrary simple block or
C               ISB1=0.
C None of the input parameters are altered.
C
C Output:
C     ISRF2...for the given point not situated inside the block ISB1 or
C             at its boundary, ISRF2 has the meaning of the index of one
C             of the surfaces bounding the simple block ISB1 and
C             separating the given point from the 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.
C             Otherwise zero.
C             If ISRF1.NE.0, the position with respect to surface ISRF1
C             is not checked, and if INSIDE=.FALSE., ISRF2=0 without any
C             checking.
C     ISB2... For ISRF1.NE.0, ISRF2.EQ.0: index of a simple block
C               neighbouring to the simple block isb1 and separated from
C               the simple block ISB1 by surface IABS(ISRF1) at the
C               given point.
C             For ISRF1.NE.0, ISRF2.NE.0: index of a simple block in
C               which the given point is situated.  In this case, the
C               given point may be situated either inside the simple
C               block ISB2 or in the vicinity of its boundary formed by
C               the surface IABS(ISRF1).  From the two possible simple
C               blocks touching the surface IABS(ISRF1) at the given
C               point, that being situated at the same side of the
C               surface ISRF1 as the simple block ISB1, is selected.
C             For ISRF1.EQ.0: index of the simple block in which the
C               given point is situated.
C             If there is no simple material block of the above
C             properties, ISB2=0.
C     ICB2... Index of the complex block in which the simple block ISB2
C             is situated.  ICB2=0 if the simple block ISB2 is not
C             situated in any complex block.  For ISB2=0: ICB2=0.
C     F...    For ISRF2.NE.0: array containing the value and partial
C               derivatives F, F1, F3, F11, F12, F22, F13, F23, F33 of
C               the function describing the surface IABS(ISRF2).
C             For ISRF2.EQ.0: undefined.
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: 1995, October 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,J,I1,I2,ISIDE1,ISIDE2,ISRF
      REAL FAUX(10)
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
C
C     Initialization:
      ISRF2=0
      IF(ISRF1.EQ.0.OR.ISB1.EQ.0) THEN
        ISIDE1=2
      ELSE
        ISIDE1=-ISIDE(ISRF1,ISB1)
      END IF
C
C     Position of the given point with respect to the given simple block
C     ISB1:
      IF(INSIDE) THEN
        IF(ISB1.NE.0) THEN
          CALL INTERF(COOR,ISRF1,ISB1,ISRF2,F)
          IF(ISRF1.EQ.0) THEN
            IF(ISRF2.EQ.0) THEN
C             The point is inside simple block ISB1:
              ISB2=ISB1
              GO TO 10
            END IF
          ELSE
            IF(ISRF2.NE.0) THEN
C             The point being situated at surface ISRF1 bounding simple
C             block ISB1 is not situated at the boundary of simple block
C             ISB1:
              ISIDE1=-ISIDE1
            END IF
          END IF
        END IF
      END IF
C
C     Search for the simple block in which the given point is
C     situated:
      I2=MAX0(ISB1-1,NSB-ISB1)
      DO 2 J=1,I2
        DO 1 I=1,-1,-2
          ISB2=ISB1+I*J
          IF(ISB2.GT.0.AND.ISB2.LE.NSB) THEN
            IF(ISRF1.NE.0) THEN
              ISIDE2=ISIDE(ISRF1,ISB2)
            END IF
            IF(ISRF1.EQ.0.OR.ISIDE1.EQ.ISIDE2) THEN
              CALL INTERF(COOR,ISRF1,ISB2,ISRF,FAUX)
              IF(ISRF.EQ.0) GO TO 10
            END IF
          END IF
    1   CONTINUE
    2 CONTINUE
      IF(.NOT.INSIDE) THEN
        DO 4 J=1,I2
          DO 3 I=1,-1,-2
            ISB2=ISB1+I*J
            IF(ISB2.GT.0.AND.ISB2.LE.NSB) THEN
              CALL INTERF(COOR,ISRF1,ISB2,ISRF,FAUX)
              IF(ISRF.EQ.0) GO TO 10
            END IF
    3     CONTINUE
    4   CONTINUE
      END IF
C     No simple block has been found:
      ISB2=0
      ICB2=0
      RETURN
C
C     Determination of the complex block:
   10 CONTINUE
      DO 12 J=1,NCB
        I1=KCB(J-1)+1
        I2=KCB(J)
        DO 11 I=I1,I2
          IF(KCB(I).EQ.ISB2) THEN
            ICB2=J
            RETURN
          END IF
   11   CONTINUE
   12 CONTINUE
C     No complex block:
      ICB2=0
C
      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 BLOCK.
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,ISRF1,ISB,ISRF2,F)
      REAL COOR(3),F(10)
      INTEGER ISRF1,ISB,ISRF2
C
C This is an auxiliary subroutine to the subroutine BLOCK.
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     ISRF1...ISRF1.NE.0: Index of the surface at which the given point
C               is situated.  The sign is ignored.
C             ISRF1.EQ.0: The given point is not situated at any
C               surface.
C     ISB...  Index of the given simple block.
C None of the input parameters are altered.
C
C Output:
C     ISRF2...Index of a surface separating the given point and the
C               simple block ISB, supplemented by a sign '+' or '-' for
C               the simple block ISB1 situated at the positive or
C               negative side of the surface, respectively.
C             ISRF2=0  if the given point is situated inside the simple
C               block ISB.
C     F...    For ISRF2.NE.0: Array containing the value and partial
C               derivatives F, F1, F3, F11, F12, F22, F13, F23, F33 of
C               the function describing the surface IABS(ISRF2) at the
C               given point.
C             For ISRF2.EQ.0: Undefined.
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: 1989, December 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER IS,JS,KS,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
        KS=KSB(IS)
        JS=IABS(KS)
        IF(JS.NE.IABS(ISRF1)) THEN
          CALL SRFC2(JS,COOR,F)
          IF(F(1)*FLOAT(KS).LT.0.) THEN
            ISRF2=KS
            RETURN
          END IF
        END IF
    1 CONTINUE
C
      ISRF2=0
      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)
*V      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)
*V      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: 1996, September 30
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
*V      CALL VAR4(0,1.)
      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
*V        CALL VAR4(0,AUX1)
*V        CALL VAR4(2,-AUX2*FOUT(2))
*V        CALL VAR4(3,-AUX2*FOUT(3))
*V        CALL VAR4(4,-AUX2*FOUT(4))
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C