C SUBROUTINE FILE 'PARM.FOR' FOR SPECIFICATION AND INTERPOLATION OF THE
C MATERIAL PARAMETERS OF THE MODEL IN RECTANGULAR GRIDS.
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
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     PARM2...SUBROUTINE EVALUATING THE MATERIAL PARAMETERS INCLUDING
C             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 SUBROUTINES PARM1 AND PARM2 SUPPORTING THE COMPLETE RAY TRACING
C ALGORITHM 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' 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 ATTENTION:
C     THE LINES DENOTED BY 'V' IN THE FIRST COLUMN SHOULD BE REMOVED
C     FROM THE SUBROUTINE PARM2 BEFORE COMPILATION.
C     IF THE MODEL VARIATIONS WITH RESPECT TO THE MODEL PARAMETERS WERE
C     TO BE CALCULATED, THESE LINES WOULD REMAIN WHILE REPLACING EACH
C     'V' IN THE FIRST COLUMN BY A SPACE.  IN SUCH A CASE, SUBROUTINES
C     VAR4 AND VAR5 OF THE FILE 'VAR.FOR' WOULD BE CALLED TO HANDLE THE
C     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 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 'VP '
C             'VS ', 'DEN', 'QP ', 'QS ', 'END'.
C     ICB...  INDEX OF THE COMPLEX BLOCK.
C (1B) SEVERAL TIMES 'INPUT DATA FOR ONE MATERIAL PARAMETER', SEE BELOW.
C     AT LEAST ONE OF THE VELOCITIES MUST BE SPECIFIED.  THE DEFAULT
C     VALUES OF PARAMETERS NOT SPECIFIED 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 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 SPECIFIED):
C             QP=0.0, QS=0.0.
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     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 (2) IVAR1,IVAR2,IVAR3,SIGMA
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 (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 DATE: 1992, DECEMBER 31
C CODED BY LUDEK KLIMES
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 (P AND S WAVE VELOCITIES, DENSITY, P AND S WAVE
C LOSS FACTORS) IN THE COMPLEX BLOCKS, 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 MAPPING OF POINTS
C IT IS NECESSARY TO CALL THE SUBROUTINE PARM2, WHICH ALSO RETURNS THE
C FIRST AND SECOND PARTIAL DERIVATIVES.  SUBROUTINE PARM1 MAY BE CALLED
C SEVERAL TIMES.  THE COMPLEX BLOCKS ARE INDEXED SUCCESIVELY, FOLLOWING
C THE COMPLEX BLOCKS DEFINED DURING THE PREVIOUS 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: 1989, DECEMBER 19
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      CHARACTER*3 TFUNCT(5)
      DATA TFUNCT/'VP ','VS ','DEN','QP ','QS '/
C
      CALL VAL1(LUN,2,NCB,5,TFUNCT)
      RETURN
      END
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 USEFULL 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 SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL FPOWER,VAL2
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 ERROR MESSAGES:
C     321...  NO VELOCITY IS DEFINED:
C             NEITHER P NOR S WAVE VELOCITY IN THE CURRENT COMPLEX BLOCK
C             IS DEFINED IN THE INPUT DATA.
C     322...  PROHIBITED MATERIAL PARAMETER:
C             P WAVE VELOCITY MUST BE POSITIVE, OTHER MATERIAL
C             PARAMETERS MUST BE NON-NEGATIVE.
C
C DATE: 1992, DECEMBER 31
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I
      REAL FAUX(10,5),POWER(5),POWER1,POWER2,POWER3,POWER4,POWER5,AUX(1)
      EQUIVALENCE (POWER(1),POWER1),(POWER(2),POWER2),(POWER(3),POWER3)
      EQUIVALENCE (POWER(4),POWER4),(POWER(5),POWER5)
C
C.......................................................................
C
      CALL VAL2(2,ICB,5,COOR,FAUX,POWER)
C
C     VELOCITIES:
      IF(POWER1.NE.0.) THEN
        IF(FAUX(1,1).LE.0.) GO TO 9
        CALL FPOWER(10,FAUX(1,1),POWER1,UP)
V       CALL VAR5(1,1)
        IF(POWER2.NE.0.) THEN
          IF(FAUX(1,2).LT.0.) GO TO 9
          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.) GO TO 9
          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
          PAUSE 'ERROR 321 IN PARM2: NO VELOCITY IS DEFINED'
          STOP
        END IF
      END IF
C
C     DENSITY:
      IF(POWER3.NE.0.) THEN
        IF(FAUX(1,3).LT.0.) GO TO 9
        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.) GO TO 9
        CALL FPOWER(1,FAUX(1,4),POWER4,AUX)
        QP=AUX(1)
        IF(POWER5.NE.0.) THEN
          IF(FAUX(1,5).LT.0.) GO TO 9
          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.) GO TO 9
          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)
      PAUSE 'ERROR 322 IN PARM2: PROHIBITED MATERIAL PARAMETER'
      STOP
      END
C
C=======================================================================
C
