C 
C Subroutine file 'parm.for' for specification and interpolation of the
C material parameters of the model in rectangular grids.
C
C Date: 1999, August 18
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     PARM1...Subroutine reading the input data for the material
C             parameters of the model.
C             PARM1
C     PARM2...Subroutine evaluating the isotropic material parameters
C             including their first and second derivatives.
C             The functions may be embedded:  the independent variable
C             of the function may be another material parameter of the
C             same complex block foregoing in the input data.
C             PARM2
C     ADD10...Auxiliary subroutine to PARM2 summing 3 arrays of
C             dimension 10.
C             ADD10
C     LIN10...Auxiliary subroutine to PARM2 evaluating the linear
C             combination of 3 arrays of dimension 10.
C             LIN10
C     PARM3...Subroutine evaluating the anisotropic material parameters
C             including their first and second derivatives.
C             PARM3
C     PARM4...Entry of subroutine PARM3 answering whether the model is
C             isotropic or anisotropic.
C             PARM4
C Subroutines PARM1, PARM2, and PARM3, supporting isotropic complete ray
C tracing algorithm, anisotropic ray tracing and other seismic modelling
C algorithms, only mediate the work of subroutines VAL1, VAL2 and FPOWER
C which must be appended.  In addition, subroutines CURVN1 (or its
C alternative CURVB1), CURV2D (or its alternative CURVBD), SURFB1,
C SURFBD, VAL3B1, VAL3BD, VGEN, TERMS, SNHCSH, TRIDEC, TRISOL, DSPLNZ,
C INTRVL from the subroutine package 'FITPACK' by Alan Kaylor Cline,
C Department of Computer Sciences, University of Texas at Austin, are
C used.  In the complete ray tracing, this software file 'parm.for' may
C be replaced by any user-defined package containing subroutines PARM1
C and PARM2 with the same number, type and meaning of their parameters
C as in this file.
C
C Note:
C     The lines denoted by '*V' in the first two columns of file
C     'parm.for' are designed to calculate the model variations with
C     respect to the model parameters.
C     File 'parmv.for', intended for the model inversion, is created
C     from 'parm.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 If model variations are taken into account:
C     Model variations are assumed to be stored while evaluating the
C     functions during the invocation of subroutine VAL of file
C     'val.for' and subsequent routines of file 'fit.for'.
C     The variations of P-wave velocity are assumed to be stored in
C     register 1 of the system VAR*, the variations of S-wave velocity
C     are assumed to be stored in register 2 of the system VAR*.
C     Variations of the density and loss (or quality) factors are not
C     considered, although they may be stored in other registers.
C     Subroutines VAR4 and VAR5 are called within the subroutine PARM2
C     in order to deal with the variations of P and S wave velocities.
C
C.......................................................................
C
C                                                    
C Input data (read in by subroutine PARM1):
C     These input data define the complex blocks.  They are read in by
C     subroutine PARM1.  The number NCB of the complex blocks to be
C     defined is an input argument of subroutine PARM1.  The data are
C     read in by the list directed input (free format).
C (1) NCB-times (i.e. once for each complex block) input data (1A)+(1B):
C (1A) TEXTG,ICB
C     Identification of the complex block.
C     TEXTG...Any string.  Its first 3 characters must differ from
C             'VP ', 'VS ', 'DEN', 'QP ', 'QS ',
C             'A11', 'A12', 'A22', 'A13', 'A23', 'A33', 'A14', 'A24',
C             'A34', 'A44', 'A15', 'A25', 'A35', 'A45', 'A55', 'A16',
C             'A26', 'A36', 'A46', 'A56', 'A66',
C             'B11', 'B12', 'B22', 'B13', 'B23', 'B33', 'B14', 'B24',
C             'B34', 'B44', 'B15', 'B25', 'B35', 'B45', 'B55', 'B16',
C             'B26', 'B36', 'B46', 'B56', 'B66',
C             'Q11', 'Q12', 'Q22', 'Q13', 'Q23', 'Q33', 'Q14', 'Q24',
C             'Q34', 'Q44', 'Q15', 'Q25', 'Q35', 'Q45', 'Q55', 'Q16',
C             'Q26', 'Q36', 'Q46', 'Q56', 'Q66', 'END'.
C     ICB...  Index of the complex block.
C (1B) Several times 'Input data for one material parameter', see below.
C     Isotropic complex block:
C       At least one of velocities 'VP ' and 'VS ' must be specified.
C       Unspecified isotropic elastic parameters ('VP ', 'VS ', 'QP ',
C       'QS ') take their default values.  Anisotropic elastic
C       parameters correspond to the isotropic medium.
C     Anisotropic complex block with given isotropic reference medium:
C       Isotropic complex block with one to all anisotropic elastic
C       parameters specified.  Unspecified anisotropic elastic
C       parameters default to the isotropic medium.
C     Anisotropic complex block:
C       At least 9 reduced anisotropic elastic parameters 'A11', 'A12',
C       'A22', 'A13', 'A23', 'A33', 'A44', 'A55', and 'A66' must be
C       specified in the anisotropic complex block.  Unspecified
C       anisotropic elastic parameters default to zeros.
C (2) TEXTE,AUX
C     End of data.
C     TEXTE...String, the first 3 characters of which must be upper-case
C             'END'.
C     AUX...  Any number or a slash.
C For an example refer to the sample input data for the model.
C
C Input data for one material parameter:
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 TEXTF), the input parameter is of the
C     type REAL.
C (1) TEXTF,POWER
C     Physical meaning of the function.
C     TEXTF...String identifying which material parameter the function
C             describes. Only the first 3 characters are significant.
C             The first 3 characters of the string must be:
C             'VP ' for P wave velocity,
C             'VS ' for S wave velocity,
C             'DEN' for density,
C             'QP ' for P wave loss factor,
C             'QS ' for S wave loss factor.
C             'A11', 'A12', 'A22', 'A13', 'A23', 'A33', 'A14', 'A24',
C             'A34', 'A44', 'A15', 'A25', 'A35', 'A45', 'A55', 'A16',
C             'A26', 'A36', 'A46', 'A56', or 'A66' for reduced (i.e.
C               divided by the density) anisotropic elastic parameters
C               (components of the real part of the symmetric 6*6
C               stiffness matrix divided by the density).
C             'Q11', 'Q12', 'Q22', 'Q13', 'Q23', 'Q33', 'Q14', 'Q24',
C             'Q34', 'Q44', 'Q15', 'Q25', 'Q35', 'Q45', 'Q55', 'Q16',
C             'Q26', 'Q36', 'Q46', 'Q56', or 'Q66' for reduced (i.e.
C               divided by the density) imaginary anisotropic elastic
C               parameters (components of the imaginary part of the
C               symmetric 6*6 stiffness matrix divided by the density).
C             All strings must be entered in uppercase.
C     POWER...The specified function is equal to the POWER-th power of
C             the material parameter.
C             Examples:
C               For P wave velocity  TEXTF='VP ' and  POWER=1.0,
C               for P wave slowness  TEXTF='VP ' and  POWER=-1.0,
C               for S wave quadratic slowness  TEXTF='VS ' and
C                 POWER=-2.0,
C               for P wave loss factor  TEXTF='QP ' and  POWER=1.0,
C               for S wave quality factor  TEXTF='QS ' and  POWER=-1.0.
C     Default values of isotropic elastic parameters not specified in
C     the input data (subroutine PARM2):
C             Isotropic complex block or anisotropic complex block with
C             given isotropic reference medium (at least one of VP and
C             VS given):
C               P wave velocity: VP=1.73205*VS,
C               S wave velocity: VS=0.57735*VP,
C             Anisotropic complex block (neither VP nor VS given):
C               Default isotropic parameters are expressed in terms of
C               anisotropic parameters to allow for application of the
C               isotropic code to anisotropic models.
C               Average isotropic medium minimizing the squared norm of
C               the difference between the Christoffel matrices,
C               averaged over propagation directions:
C               AP=A11+A22+A33,
C               AL=A12+A13+A23,
C               AS=A44+A55+A66,
C               VP=SQRT( (3.*AP+2.*AL+4.*AS)/15 ),
C               VS=SQRT( (   AP  - AL+3.*AS)/15 ),
C               This reference isotropic medium is designed especially
C               to serve as a basis for anisotropic perturbations. From
C               this point of view, the variations of the these default
C               isotropic material parameters with respect to the model
C               parameters, evaluated by means of subroutines VAR*, are
C               unnecessary.  No variations of the above defaults are
C               thus set in PARM2.
C             In any case:
C               Density: RHO=1,
C               P wave loss factor (if S wave loss factor is specified):
C                       QP=1.333333*QS*(US(1)/UP(1))**2,
C               S wave loss factor (if P wave loss factor is specified):
C                       QS=0.750000*QP*(UP(1)/US(1))**2,
C               P and S wave loss factors (if none of them is given):
C                       QP=0.0, QS=0.0.
C     Default values of anisotropic elastic parameters not specified in
C     the input data (subroutine PARM3):
C             Isotropic complex block or anisotropic complex block with
C             given isotropic reference medium (at least one of VP and
C             VS given):
C               Default anisotropic parameters are expressed in terms of
C               isotropic parameters for the sake of continuity with
C               isotropic models:
C               A11=A22=A33= VP*VP,
C               A12=A13=A23= VP*VP-2*VS*VS,
C               A44=A55=A66= VS*VS,
C               A14=A24=A34=A15=A25=A35=A45=A16=A26=A36=A46=A56= 0.
C             Anisotropic complex block (neither VP nor VS given):
C               Aij=0.
C             In any case:
C               RHO=1,
C               Qij=0 in this version.
C (2) IVAR1,IVAR2,IVAR3,SIGMA,POWERW,/
C     The form of the function.
C     IVAR1,IVAR2,IVAR3... Denote the form of the function. The function
C             must be of the form
C               F(X1,X2,X3) = W(A1,A2,A3)-B1-B2-B3   .
C             X1, X2, X3 are the general coordinates. Each of A1, A2,
C             A3, B1, B2, B3 must be either:  (a) one of general
C             coordinates X1, X2, X3, (b) another previously defined
C             function F(X1,X2,X3) of the same complex block,  or (c)
C             must be left out.  At most 3 of parameters A1-B3 may be of
C             kind (a) or (b).  Note that IVAR1 controls the type of A1
C             and B1, IVAR2 controls the type of A2 and B2, IVAR3
C             controls the type of A3 and B3.
C             For IVAR1.EQ.0: A1, B1 are empty (left out).
C             For IVAR1.EQ.1: A1=X1, B1 is empty.
C             For IVAR1.EQ.2: A1=X2, B1 is empty.
C             For IVAR1.EQ.3: A1=X3, B1 is empty.
C             For IVAR1.GE.4: A1=F(X1,X2,X3), where F(X1,X2,X3) is
C               another function of the same complex block defined in
C               the input data as the (IVAR1-3)-th function of the
C               complex block. B1 is empty.
C               Example:
C                 If density=1.7+0.2*VP then the interpolated function
C                 is W(A1,A2,A3)=1.7+0.2*A1 with the independent
C                 variable A1=VP(X1,X2,X3).  This is specified by
C                 IVAR1=4, IVAR2=0, IVAR3=0 if VP is the first read in
C                 parameter, by IVAR1=5, IVAR2=0, IVAR3=0 if VP is the
C                 second read in parameter, etc.
C                 The possible alternatives are W(A1,A2,A3)=1.7+0.2*A2
C                 with A2=VP(X1,X2,X3) specified by IVAR1=0,
C                 IVAR2=(4 or 5 or the like), IVAR3=0, and
C                 W(A1,A2,A3)=1.7+0.2*A3 with A3=VP(X1,X2,X3) specified
C                 by IVAR1=0, IVAR3=2, IVAR3=(4 or 5 or the like).
C             For IVAR1.EQ.-1: B1=X1, A1 is empty.
C             For IVAR1.EQ.-2: B1=X2, A1 is empty.
C             For IVAR1.EQ.-3: B1=X3, A1 is empty.
C             For IVAR1.GE.-4: B1=F(X1,X2,X3), where F(X1,X2,X3) is
C               another function of the same complex block defined in
C               the input data as the (-IVAR1-3)-th function of the
C               complex block. A1 is empty.
C             The meaning of the parameters IVAR2, IVAR3 is similar.
C             Examples:
C             IVAR1: IVAR2: IVAR3: the form of the function:
C               1      2      3     F(X1,X2,X3)=W(X1,X2,X3)
C               3      1      2     F(X1,X2,X3)=W(X3,X1,X2)
C               1      2      0     F(X1,X2,X3)=W(X1,X2)
C               5      0      0     F(X1,X2,X3)=W(F2(X1,X2,X3)), where
C             F2(X1,X2,X3) is the second material parameter of the
C             complex block defined in the input data. Function W is
C             interpolated by means of splines under tension.
C     SIGMA...Is the tension factor (its sign is ignored).  This value
C             indicates the curviness desired.  If ABS(SIGMA) is nearly
C             zero (e.g. 0.001), the resulting surface is approximately
C             the tensor product of cubic splines.  If ABS(SIGMA) is
C             large (e.g. 50.), the resulting surface is approximately
C             tri-linear.  If SIGMA equals zero, tensor products of
C             cubic splines result.  A recommended value for SIGMA is
C             approximately 1. In absolute value.
C     POWERW..Given grid values (7) correspond to the POWERW-th power of
C             interpolated function W.  The given grid values (7) are
C             thus raised to the (1/POWERW)-th power immediately after
C             reading and then interpolated.
C     /...    Obligatory slash at the end of line for future extensions.
C     Default: IVAR1=0, IVAR2=0, IVAR3=0, SIGMA=0, POWERW=1.
C (3) NX(1),...,NX(NVAR)
C     The numbers of grid coordinates for the interpolation.
C     This input is performed if at least one of IVAR1, IVAR2, IVAR3 is
C     positive.
C     Each of NX(1),...,NX(NVAR) corresponds to one positive value of
C     IVAR1, IVAR2, IVAR3 and specifies the number of grid coordinates
C     corresponding to that independent variable of function W, see (2).
C     The sign of NX(1),...,NX(NVAR) is ignored.  NVAR (.LE.3) is the
C     number of positive values of the above quantities IVAR1, IVAR2,
C     IVAR3, i.e. the number of independent variables of function W,
C     see (1).
C (4) X1(1),...,X1(NX(1))
C     The grid coordinates corresponding to the first independent
C     variable of function W, see (2).
C     This input is performed if NX(1) is specified, see (3), and is not
C     zero. The grid coordinates may be specified in any order.
C (5) X2(1),...,X2(NX(2))
C     The grid coordinates corresponding to the second independent
C     variable of function W, see (2).
C     This input is performed if NX(2) is specified, see (3), and is not
C     zero. The grid coordinates may be specified in any order.
C (6) X3(1),...,X3(NX(3))
C     The grid coordinates corresponding to the third independent
C     variable of function W, see (2).
C     This input is performed if NX(3) is specified, see (3), and is not
C     zero. The grid coordinates may be specified in any order.
C (7) (((W(I,J,K),I=1,MAX(NX(1),1)),J=1,MAX(NX(2),1)),K=1,MAX(NX(3),1))
C     the values of function W at grid points. Function value W(i,j,k)
C     corresponds to point (X1(I),X2(J),X3(K)).
C
C=======================================================================
C
C     
C
      SUBROUTINE PARM1(LUN,NCB)
      INTEGER LUN,NCB
C
C This subroutine reads the input data for the distributions of the
C material parameters:
C In the isotropic complex blocks:
C     P and S wave velocities, density, P and S wave loss factors
C in the anisotropic complex blocks:
C     21 real parts of the reduced (divided by the density) elastic
C     parameters, 21 corresponding imaginary parts, density.
C Reference isotropic P and S velocities and loss factors may also be
C specified together with anisotropic material parameters.
C When reading the data, subroutine PARM1 also determines the parameters
C necessary to compute an interpolatory function on a three dimensional
C rectangular grid, and stores them in the memory.  The function
C determined can be represented as a tensor product of splines under
C tension.  The functions may be embedded.  For actual interpolation of
C material parameters it is necessary to call subroutine PARM2 for
C isotropic model (or mean isotropic model corresponding to the
C anisotropic model), or subroutine PARM3 for anisotropic material
C parameters of the isotropic or anisotropic model.  Subroutines PARM2
C and PARM3 also return the first and second partial derivatives of
C propagation velocities or reduced elastic parameters.  Subroutine
C PARM1 may be called several times.  The complex blocks are indexed
C successively, following the complex blocks defined during the previous
C invocations.
C
C Input:
C     LUN...  Logical unit number of the external input device
C             containing the input data.
C     NCB...  Number of the material complex blocks for which the input
C             data are specified during the current invocation of PARM1.
C None of the input parameters are altered.
C
C No output.
C
C Subroutines and external functions required:
      EXTERNAL VAL1
C     VAL1, SORTV, READV... File 'val.for'.
C     CURVN1 or CURVB1 (alternatives), SURFB1, VAL3B1, SNHCSH, VGEN,
C              TERMS, TRIDEC, TRISOL... Subroutine package 'FITPACK'
C              (file 'fit.for').
C
C Date: 1995, December 17
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      CHARACTER*3 TFUNCT(47)
      DATA TFUNCT/'VP ','VS ','DEN','QP ','QS ',
     *'A11','A12','A22','A13','A23','A33','A14','A24','A34','A44',
     *'A15','A25','A35','A45','A55','A16','A26','A36','A46','A56','A66',
     *'Q11','Q12','Q22','Q13','Q23','Q33','Q14','Q24','Q34','Q44',
     *'Q15','Q25','Q35','Q45','Q55','Q16','Q26','Q36','Q46','Q56','Q66'/
C
      CALL VAL1(LUN,2,NCB,47,TFUNCT)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PARM2(ICB,COOR,UP,US,RHO,QP,QS)
      INTEGER ICB
      REAL COOR(3),UP(10),US(10),RHO,QP,QS
C
C This subroutine evaluates P and S wave velocities, density, and P and
C S wave loss factors at a given point.  The three first and six second
C partial derivatives of the velocities are also evaluated.  The
C specified functions are represented as a tensor product of splines
C under tension.  The parameters may be dependent either on the general
C coordinates or on the distribution of another parameter, e.g.
C VS=0.577*VP or RHO=1.7+0.2*VP, where VP, VS and RHO are P and S
C velocities and density.  The coefficients of these functions are
C prepared in subroutine PARM1, in which the input data concerning the
C distribution of individual parameters within each complex block are
C read in.  The default values of parameters not specified in the input
C data are:
C     P wave velocity: VP=1.73205*VS,
C     S wave velocity: VS=0.57735*VP,
C     Density: RHO=1.0,
C     P wave loss factor (if the S wave loss factor is specified):
C             QP=1.333333*QS*(US(1)/UP(1))**2,
C     S wave loss factor (if the P wave loss factor is specified):
C             QS=0.750000*QP*(UP(1)/US(1))**2,
C     P and S wave loss factors (if none of them is specified):
C             QP=0.0, QS=0.0.
C Attention: The above default values are reasonable only if arguments
C     UP and US of this subroutine are velocities, and QP and QS loss
C     factors.  In other words, these default settings are useful only
C     if NEXPV=1 and NEXPQ=1 in the input data set model, line (2).
C Note that at least one of the velocities must be specified in the
C input data. P wave velocity must be positive, other material
C parameters must be non-negative.
C
C Input:
C     ICB...  Index of a complex block.
C     COOR... Array containing coordinates X1, X2, X3 of the given
C             point.
C None of the input parameters are altered.
C
C Output:
C     UP,US...Powers of the P and S wave velocities (the exponent of the
C             power is NEXPV, see the input data for the model) and
C             their first and second partial derivatives in order U, U1,
C             U2, U3, U11, U12, U22, U13, U23, U33, at the given point.
C     RHO...  Density at the given point.
C     QP,QS...Powers of the P and S wave loss factors (the exponent of
C             the power is NEXPQ, see the input data for the model) at
C             the given point.
C
C Common block /MODELC/ (to use NEGPAR):
      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:
      INTEGER  LUWARN
      EXTERNAL LUWARN,ERROR,FPOWER,VAL2
C     LUWARN,ERROR...File 'error.for'.
C     FPOWER...File 'model.for'.
C     VAL2...  File 'val.for'.
C     CURV2D or CURVBD (alternatives), SURFBD, VAL3BD, SNHCSH, DSPLNZ,
C              INTRVL... Subroutine package 'FITPACK' (file 'fit.for').
C
C Date: 1999, August 16
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL FAUX(10,47),POWER(47),POWER1,POWER2,POWER3,POWER4,POWER5
      EQUIVALENCE (POWER(1),POWER1),(POWER(2),POWER2),(POWER(3),POWER3)
      EQUIVALENCE (POWER(4),POWER4),(POWER(5),POWER5)
      REAL AUX(50)
C
C     Constants:
      REAL C1,C2,C3,C4
      PARAMETER (C1=1./15.,C2=2./15.,C3=3./15.,C4=4./15.)
C
C.......................................................................
C
      CALL VAL2(2,IABS(ICB),47,COOR,FAUX,POWER)
C
C     Velocities:
      IF(POWER1.NE.0.) THEN
        IF(FAUX(1,1).LE.0.) THEN
          IF(NEGPAR.EQ.0) GO TO 9
        END IF
        CALL FPOWER(10,FAUX(1,1),POWER1,UP)
*V      CALL VAR5(1,1)
        IF(POWER2.NE.0.) THEN
          IF(FAUX(1,2).LT.0.) THEN
            IF(NEGPAR.EQ.0) GO TO 9
          END IF
          CALL FPOWER(10,FAUX(1,2),POWER2,US)
*V        CALL VAR5(2,2)
        ELSE
          DO 1 I=1,10
            US(I)=0.57735*UP(I)
    1     CONTINUE
*V        CALL VAR4(0,0.57735)
*V        CALL VAR5(2,1)
        END IF
      ELSE
        IF(POWER2.NE.0.) THEN
          IF(FAUX(1,2).LE.0.) THEN
            IF(NEGPAR.EQ.0) GO TO 9
          END IF
          CALL FPOWER(10,FAUX(1,2),POWER2,US)
*V        CALL VAR5(2,2)
          DO 2 I=1,10
            UP(I)=1.73205*US(I)
    2     CONTINUE
*V        CALL VAR4(0,1.73205)
*V        CALL VAR5(1,2)
        ELSE
          IF(POWER(06).NE.0..AND.POWER(07).NE.0..AND.
     *       POWER(08).NE.0..AND.POWER(09).NE.0..AND.
     *       POWER(10).NE.0..AND.POWER(11).NE.0..AND.
     *       POWER(15).NE.0..AND.POWER(20).NE.0..AND.
     *       POWER(26).NE.0.) THEN
C           Isotropic reference medium to the anisotropic material:
            CALL FPOWER(10,FAUX(1,06),.5*POWER(06),AUX(01))
            CALL FPOWER(10,FAUX(1,08),.5*POWER(08),AUX(31))
            CALL FPOWER(10,FAUX(1,11),.5*POWER(11),AUX(41))
            CALL ADD10(AUX(01),AUX(31),AUX(41))
            CALL FPOWER(10,FAUX(1,07),.5*POWER(07),AUX(11))
            CALL FPOWER(10,FAUX(1,09),.5*POWER(09),AUX(31))
            CALL FPOWER(10,FAUX(1,10),.5*POWER(10),AUX(41))
            CALL ADD10(AUX(11),AUX(31),AUX(41))
            CALL FPOWER(10,FAUX(1,15),.5*POWER(15),AUX(21))
            CALL FPOWER(10,FAUX(1,20),.5*POWER(20),AUX(31))
            CALL FPOWER(10,FAUX(1,26),.5*POWER(26),AUX(41))
            CALL ADD10(AUX(21),AUX(31),AUX(41))
            CALL LIN10(C3,AUX(01), C2,AUX(11),C4,AUX(21),AUX(31))
            CALL LIN10(C1,AUX(01),-C1,AUX(11),C3,AUX(21),AUX(41))
            CALL FPOWER(10,AUX(31),2.,UP)
            CALL FPOWER(10,AUX(41),2.,US)
          ELSE
C           321
            CALL ERROR('321 in PARM2: No velocity is defined')
C           Neither P nor S wave velocity in the current complex block
C           is defined in the input data.
          END IF
        END IF
      END IF
C
C     Density:
      IF(POWER3.NE.0.) THEN
        IF(FAUX(1,3).LT.0.) THEN
          IF(NEGPAR.EQ.0) GO TO 9
        END IF
        CALL FPOWER(1,FAUX(1,3),POWER3,AUX)
        RHO=AUX(1)
      ELSE
        RHO=1.
      END IF
C
C     Loss factors:
      IF(POWER4.NE.0.) THEN
        IF(FAUX(1,4).LT.0.) THEN
          IF(NEGPAR.EQ.0) GO TO 9
        END IF
        CALL FPOWER(1,FAUX(1,4),POWER4,AUX)
        QP=AUX(1)
        IF(POWER5.NE.0.) THEN
          IF(FAUX(1,5).LT.0.) THEN
            IF(NEGPAR.EQ.0) GO TO 9
          END IF
          CALL FPOWER(1,FAUX(1,5),POWER5,AUX)
          QS=AUX(1)
        ELSE
          IF(US(1).GT.0.) THEN
            QS=0.750000*QP*(UP(1)/US(1))**2
          ELSE
            QS=0.
          END IF
        END IF
      ELSE
        IF(POWER5.NE.0.) THEN
          IF(FAUX(1,5).LT.0.) THEN
            IF(NEGPAR.EQ.0) GO TO 9
          END IF
          CALL FPOWER(1,FAUX(1,5),POWER5,AUX)
          QS=AUX(1)
          QP=1.333333*QS*(US(1)/UP(1))**2
        ELSE
          QP=0.
          QS=0.
        END IF
      END IF
      RETURN
C
    9 CONTINUE
        WRITE(*,'('' X='',F9.3,'' Y='',
     *        F9.3,'' Z='',F9.3/''   VP='',F7.3,'' VS='',F7.3,'' RO='',
     *        F7.3,'' QP='',F7.3,'' QS='',F7.3)') COOR,(FAUX(1,I),I=1,5)
        IF(LUWARN(0).NE.0) THEN
          WRITE(LUWARN(0),'('' X='',F9.3,'' Y='',
     *        F9.3,'' Z='',F9.3/''   VP='',F7.3,'' VS='',F7.3,'' RO='',
     *        F7.3,'' QP='',F7.3,'' QS='',F7.3)') COOR,(FAUX(1,I),I=1,5)
        END IF
C       322
        CALL ERROR('322 in PARM2: Prohibited material parameter')
C       P wave velocity must be positive, other material parameters must
C       be non-negative.
      END
C
C-----------------------------------------------------------------------
C
C     
C
      SUBROUTINE ADD10(A,B,C)
      REAL A(10),B(10),C(10)
C
C Auxiliary subroutine to PARM2 summing 3 arrays of dimension 10.
C
C.......................................................................
C
      INTEGER I
C
      DO 10 I=1,10
        A(I)=A(I)+B(I)+C(I)
   10 CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C
C     
C
      SUBROUTINE LIN10(C1,A1,C2,A2,C3,A3,A)
      REAL C1,A1(10),C2,A2(10),C3,A3(10),A(10)
C
C Auxiliary subroutine to PARM2 evaluating linear combination of
C 3 arrays of dimension 10.
C
C.......................................................................
C
      INTEGER I
C
      DO 10 I=1,10
        A(I)=C1*A1(I)+C2*A2(I)+C3*A3(I)
   10 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PARM3(ICB,COOR,A,RHO,Q)
      INTEGER ICB
      REAL COOR(3),A(10,21),RHO,Q(21)
C
C This subroutine is redundant for isotropic seismic models and codes.
C
C This subroutine evaluates 21 real parts of the reduced (divided by the
C density) elastic parameters, 21 corresponding imaginary parts, and the
C density at a given point.  The three first and six second partial
C derivatives of the 21 real parts of the reduced elastic parameters are
C also evaluated.  The specified functions are represented as a tensor
C product of splines under tension.  The parameters may be dependent
C either on the general coordinates or on the distribution of another
C parameter.  The coefficients of these functions are prepared in
C subroutine PARM1, in which the input data concerning the distribution
C of individual parameters within each complex block are read in.
C Variations of real parts of the reduced elastic parameters Aij with
C respect to model parameters (evaluated by subroutines VAR*) are
C stored in registers 6 to 26.
C
C Input:
C     ICB...  Index of a complex block.
C     COOR... Array containing coordinates X1, X2, X3 of the given
C             point.
C None of the input parameters are altered.
C
C Output:
C     A...    Values, first and second partial derivatives of real
C             parts of 21 reduced (divided by the density) elastic
C             parameters.  The order of the value, first and second
C             partial derivatives of each parameter Aij is:
C             Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33.
C             The order of parameters (second array index) is:
C             A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45,
C             A55,A16,A26,A36,A46,A56,A66.
C     RHO...  Density at the given point.
C     Q...    Imaginary parts of 21 reduced elastic parameters at the
C             given point, ordered as
C             Q11,Q12,Q22,Q13,Q23,Q33,Q14,Q24,Q34,Q44,Q15,Q25,Q35,Q45,
C             Q55,Q16,Q26,Q36,Q46,Q56,Q66.
C
C-----------------------------------------------------------------------
C
C     
C
C     ENTRY PARM4(ISOFLG)
      INTEGER ISOFLG
C
C Entry of subroutine PARM3 answering whether the model is isotropic.
C
C No input.
C
C Output:
C     ISOFLG..ISOFLG=0: Anisotropic model.
C             ISOFLG=1: Isotropic model.
C
C Common block /MODELC/ (to use NEGPAR in PARM3 and BOUNDM in PARM4):
      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 ERROR,FPOWER,VAL2
C     ERROR...File 'error.for'.
C     FPOWER...File 'model.for'.
C     VAL2...  File 'val.for'.
C     CURV2D or CURVBD (alternatives), SURFBD, VAL3BD, SNHCSH, DSPLNZ,
C              INTRVL... Subroutine package 'FITPACK' (file 'fit.for').
C
C Date: 1999, August 16
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER KA(21),IA,I,I1,I2,IERR
      REAL FAUX(10,47),POWER(47),POWER3,AUX(10)
      EQUIVALENCE (POWER(3),POWER3)
C     Order of processing the reduced elastic parameters:
      DATA KA/1,3,6,21,15,10,2,4,5,20,19,14,16,17,18,11,12,13,7,8,9/
C
C.......................................................................
C
      CALL VAL2(2,IABS(ICB),47,COOR,FAUX,POWER)
C
C     Density:
      IF(POWER3.NE.0.) THEN
        IF(FAUX(1,3).LT.0.) THEN
          IERR=3
          IF(NEGPAR.EQ.0) GO TO 49
        END IF
        CALL FPOWER(1,FAUX(1,3),POWER3,AUX)
        RHO=AUX(1)
      ELSE
        RHO=1.
      END IF
C
C     Real-valued elastic parameters:
      DO 19 I2=1,21
        IA=KA(I2)
        I=IA+5
        IF(POWER(I).NE.0.) THEN
C         Element of the stiffness matrix is specified in the data:
          IERR=3
          IF(I2.LE.6) THEN
            IF(I2.LE.3) THEN
              IF(FAUX(1,I).LE.0.) THEN
                IERR=I
                IF(NEGPAR.EQ.0) GO TO 49
              END IF
            ELSE
              IF(FAUX(1,I).LT.0.) THEN
                IERR=I
                IF(NEGPAR.EQ.0) GO TO 49
              END IF
            END IF
          END IF
          CALL FPOWER(10,FAUX(1,I),POWER(I),A(1,IA))
*V        CALL VAR5(I,I)
        ELSE IF(I2.LE.3) THEN
C         Default diagonal stiffnesses A11,A22,A33:
          IF(POWER(1).NE.0.) THEN
C           Aij=VP*VP:
            IF(FAUX(1,1).LE.0.) THEN
              IERR=1
              IF(NEGPAR.EQ.0) GO TO 49
            END IF
            CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA))
*V          CALL VAR5(I,1)
          ELSE IF(POWER(2).NE.0.) THEN
C           Aij=3.*VS*VS:
            IF(FAUX(1,2).LT.0.) THEN
              IERR=2
              IF(NEGPAR.EQ.0) GO TO 49
            END IF
            CALL FPOWER(10,FAUX(1,2),0.5*POWER(2),A(1,IA))
*V          CALL VAR5(I,2)
            DO 11 I1=1,10
              A(I1,IA)=3.*A(I1,IA)
   11       CONTINUE
*V          CALL VAR4(0,3.)
*V          CALL VAR5(I,I)
          ELSE
C           323
            CALL ERROR('323 in PARM3: Undefined elastic parameter')
C           If neither isotropic velocity VP nor VS is specified in the
C           input data, A11, A22 and A33  must be specified in order to
C           use this subroutine.
          END IF
        ELSE IF(I2.LE.6) THEN
C         Default diagonal stiffnesses A66,A55,A44:
          IF(POWER(2).NE.0.) THEN
C           Aij=VS*VS:
            IF(FAUX(1,2).LT.0.) THEN
              IERR=2
              IF(NEGPAR.EQ.0) GO TO 49
            END IF
            CALL FPOWER(10,FAUX(1,2),0.5*POWER(2),A(1,IA))
*V          CALL VAR5(I,2)
          ELSE IF(POWER(1).NE.0.) THEN
C           Aij=.333333*VP*VP:
            IF(FAUX(1,1).LE.0.) THEN
              IERR=1
              IF(NEGPAR.EQ.0) GO TO 49
            END IF
            CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA))
*V          CALL VAR5(I,1)
            DO 12 I1=1,10
              A(I1,IA)=.333333*A(I1,IA)
   12       CONTINUE
*V          CALL VAR4(0,.333333)
*V          CALL VAR5(I,I)
          ELSE
C           324
            CALL ERROR('324 in PARM3: Undefined elastic parameter')
C           If neither isotropic velocity VP nor VS is specified in the
C           input data, A44, A55 and A66  must be specified in order to
C           use this subroutine.
          END IF
        ELSE IF(I2.LE.9) THEN
C         Default non-diagonal stiffnesses A12,A13,A23:
          IF(POWER(1).NE.0.) THEN
            IF(FAUX(1,1).LE.0.) THEN
              IERR=1
              IF(NEGPAR.EQ.0) GO TO 49
            END IF
            IF(POWER(2).NE.0.) THEN
C             Aij=VP*VP-2*VS*VS:
              IF(FAUX(1,2).LT.0.) THEN
                IERR=2
                IF(NEGPAR.EQ.0) GO TO 49
              END IF
              CALL FPOWER(10,FAUX(1,2),0.5*POWER(1),AUX)
*V            CALL VAR5(I,2)
*V            CALL VAR4(0,-2.)
*V            CALL VAR5(I,I)
              CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA))
*V            CALL VAR5(I,1)
              DO 13 I1=1,10
                A(I1,IA)=A(I1,IA)-2.*AUX(I1)
   13         CONTINUE
            ELSE
C             Aij=.333333*VP*VP:
              CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA))
*V            CALL VAR5(I,1)
              DO 14 I1=1,10
                A(I1,IA)=.333333*A(I1,IA)
   14         CONTINUE
*V            CALL VAR4(0,.333333)
*V            CALL VAR5(I,I)
            END IF
          ELSE
            IF(POWER(2).NE.0.) THEN
C             Aij=VS*VS:
              IF(FAUX(1,2).LT.0.) THEN
                IERR=2
                IF(NEGPAR.EQ.0) GO TO 49
              END IF
              CALL FPOWER(10,FAUX(1,2),0.5*POWER(2),A(1,IA))
*V            CALL VAR5(I,2)
            ELSE
C             325
              CALL ERROR('325 in PARM3: Undefined elastic parameter')
C             If neither isotropic velocity VP nor VS is specified in
C             the input data, A45, A46 and A56  must be specified in
C             order to use this subroutine.
            END IF
          END IF
        ELSE
C         Default non-diagonal stiffnesses
C         A56,A46,A45,A16,A26,A36,A15,A25,A35,A14,A24,A34:
          DO 18 I1=1,10
            A(I1,IA)=0.
   18     CONTINUE
        END IF
   19 CONTINUE
C
C     Imaginary parts of the elastic parameters:
      DO 29 I2=1,21
        I=I2+26
        IF(POWER(I).NE.0.) THEN
          IF(I2.LE.6) THEN
            IF(FAUX(1,I).LT.0.) THEN
              IERR=I
              IF(NEGPAR.EQ.0) GO TO 49
            END IF
          END IF
          CALL FPOWER(1,FAUX(1,I),POWER(I),Q(I2))
        ELSE
C         Default:
          Q(I2)=0.
        END IF
   29 CONTINUE
C
C     Check for positive semidefinitness:
C     *** not coded ***
C       Error 327: Prohibited elastic parameter:
C       Both real and imaginary parts of the 6*6 matrix of elastic
C       parameters (stiffness matrix) must be positively
C       semi-definite, and the 3*3 upper-left minor of its real
C       part must be positively definite.  The density must be
C       positive.
      RETURN
C
   49 CONTINUE
        WRITE(*,'('' X='',F9.3,'' Y='',F9.3,'' Z='',F9.3,
     *             '' FAUX(1,'',I2,'')='',F7.3)') COOR,IERR,FAUX(1,IERR)
C       326
        CALL ERROR('326 in PARM3: Prohibited elastic parameter')
C       Following parameters must be positive:
C         FAUX(1, 1)=VP **POWER( 1)
C         FAUX(1, 6)=A11**POWER( 6)
C         FAUX(1, 8)=A22**POWER( 8)
C         FAUX(1,11)=A33**POWER(11)
C       Following parameters must be non-negative:
C         FAUX(1, 2)=VS **POWER( 2)
C         FAUX(1, 3)=RHO**POWER( 3)
C         FAUX(1,15)=A44**POWER(15)
C         FAUX(1,20)=A55**POWER(20)
C         FAUX(1,26)=A66**POWER(26)
C         FAUX(1,27)=Q11**POWER(27)
C         FAUX(1,29)=Q22**POWER(29)
C         FAUX(1,32)=Q33**POWER(32)
C         FAUX(1,36)=Q44**POWER(36)
C         FAUX(1,41)=Q55**POWER(41)
C         FAUX(1,47)=Q66**POWER(47)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY PARM4(ISOFLG)
C
      ISOFLG=1
      AUX(1)=(BOUNDM(1)+BOUNDM(2))/2.
      AUX(2)=(BOUNDM(3)+BOUNDM(4))/2.
      AUX(3)=(BOUNDM(5)+BOUNDM(6))/2.
      DO 62 I2=1,NCB
        CALL VAL2(2,I2,47,AUX,FAUX,POWER)
        DO 61 I1=5,47
          IF(POWER(I1).NE.0.) THEN
            ISOFLG=0
          END IF
   61   CONTINUE
   62 CONTINUE
      RETURN
      END
C
C=======================================================================
C