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