C SUBROUTINE FILE 'MODEL.FOR' TO SPECIFY A SEISMIC MODEL. C C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK C C THIS FILE CONSISTS OF THE FOLLOWING EXTERNAL PROCEDURES: C MODELB..BLOCK DATA SUBROUTINE DEFINING COMMON BLOCKS MODELT AND C MODELC TO STORE THE DATA DESCRIBING THE MODEL. 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 FINALY 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 NSRFC...INTEGER FUNCTION RETURNING THE NUMBER OF SURFACES COVERING C STRUCTURAL INTERFACES. C BLOCK...SUBROUTINE DESIGNED TO DETERMINE THE MUTUAL POSITION OF A C POINT AND A SIMPLE AND A COMPLEX BLOCK. C ISIDE...AUXILIARY FUNCTION TO SUBROUTINE BLOCK. C INTERF..AUXILIARY SUBROUTINE TO SUBROUTINE BLOCK. C VELOC...SUBROUTINE TRANSFORMING THE VALUES OF A MEDIUM PARAMETERS C INTO VELOCITIES AND LOSS FACTORS. 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 C ATTENTION: C THE LINES DENOTED BY 'V' IN THE FIRST COLUMN SHOULD BE REMOVED C FROM THE SUBROUTINES VELOC AND FPOWER 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 INPUT DATA 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 C KOORS...SPECIFIES THE TYPE OF THE RIGHT-HANDED COORDINATE SYSTEM C (SEE THE SUBROUTINE FILE 'METRIC.FOR'): C KOORS.LE.0: CARTESIAN COORDINATES. C KOORS.EQ.1: POLAR SPHERICAL COORDINATES (X1,X2,X3)= C (COLATITUDE,LONGITUDE,RADIUS). C KOORS.GE.2: GEOGRAPHIC SPHERICAL COORDINATES (X1,X2,X3)= C (LONGITUDE,LATITUDE,RADIUS). C NEXPV,NEXPQ... SPECIFY EXPONENTS OF THE POWER OF VELOCITIES C (NEXPV) AND LOSS FACTORS (NEXPQ) IN INPUT DATA. FOR C EXAMPLE, UNIT VALUES OF NEXPV AND NEXPQ INDICATE THAT THE C PARAMETERS OF THE MEDIUM ARE VELOCITIES AND LOSS FACTORS, C INDICES EQUAL -1 INDICATE RECIPROCAL VALUES OF THESE C QUANTITIES, I.E. SLOWNESSES AND QUALITY FACTORS. C WHILE USING THE BASIC SUBMITTED VERSION OF THE SUBROUTINE C FILE 'PARM.FOR', THE SPECIFICATION OF NEXPV=1, NEXPQ=1 IS C HIGHLY RECOMMENDED. C DEFAULT: KOORS=0, NEXPV=1, NEXPQ=1. 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 (5) NSB C NUMBER OF MATERIAL SIMPLE BLOCKS IN THE MODEL. THE BLOCKS ARE C INDEXED SEQUENTIALLY IN ANY ORDER BY POSITIVE INTEGERS ISB FROM 1 C TO NSB. FREE-SPACE BLOCKS ARE NOT INDEXED. 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 COMPLEX BLOCK, THE INDICES OF SIMPLE BLOCKS FORMING THE C COMPLEX BLOCK. THE INDICES MAY BE SPECIFIED IN AN ARBITRARY ORDER C AND MUST BE TERMINATED BY A SLASH. 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 (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 FOR AN EXAMPLE REFER TO THE SAMPLE INPUT DATA FOR THE MODEL. 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 FOLLOWING SUBROUTINE: C ------------------------------------------------------------------ BLOCK DATA MODELB CHARACTER*80 TEXTM COMMON/MODELT/TEXTM SAVE /MODELT/ INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB SAVE /MODELC/ END C ------------------------------------------------------------------ C TEXTM...THE NAME OF THE MODEL. STRING OF 80 CHARACTERS. C NEXPV,NEXPQ... SPECIFY EXPONENTS OF THE POWER OF VELOCITIES C (NEXPV) AND Q-FACTORS (NEXPQ) IN INPUT DATA. FOR EXAMPLE, C UNIT VALUES OF NEXPV AND NEXPQ INDICATE THAT THE C PARAMETERS OF THE MEDIUM ARE VELOCITIES AND Q FACTORS, C INDICES EQUAL -1 INDICATE RECIPROCAL VALUES OF THESE C QUANTITIES, I.E. SLOWNESSES AND LOSS FACTORS. C BOUNDM..BOUNDARIES X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX OF THE C MODEL. C NSRFCS..NUMBER OF SMOOTH SURFACES IN THE MODEL. THE SURFACES C ARE INDEXED SEQUENTIALLY BY POSITIVE INTEGERS, FROM 1 TO C NSRFCS. NSRFCS IS THE STORAGE LOCATION FOR NSRFC. C NSB... NUMBER OF MATERIAL SIMPLE BLOCKS IN THE MODEL. THE BLOCKS C ARE INDEXED SEQUENTIALLY BY POSITIVE INTEGERS ISB FROM 1 C TO NSB. FREE-SPACE BLOCKS ARE NOT INDEXED. C KSB... CONTAINS THE INDICES OF THE SURFACES BOUNDING INDIVIDUAL C SIMPLE BLOCKS. KSB(ISB), FOR ISB = 1 TO NSB, SPECIFY THE C PARTITION OF ARRAY KSB(NSB+1:NSB+NS) AMONG THE SIMPLE C BLOCKS. HERE NS IS THE TOTAL NUMBER OF ALL OCCURENCIES OF C THE INDICES OF THE SURFACES BOUNDING ALL INDIVIDUAL SIMPLE C BLOCKS IN THE INPUT DATA. THE INDICES OF THE SURFACES C BOUNDING INDIVIDUAL SIMPLE BLOCKS ARE STORED FROM C KSB(NSB+1) TO KSB(NSB+NS). THE LOCATIONS KSB(NSB+NS+1:MSB) C ARE UNDEFINED. IT MUST BE NSB+NS.LT.MSB. THE INDICES OF C THE SURFACES BOUNDING THE SIMPLE BLOCK ISB ARE STORED IN C KSB(I1) TO KSB(I2), WITH C I1 = KSB(ISB-1)+1 , C I2 = KSB(ISB) , C WHERE KSB(ISB-1)=NSB FOR ISB=1. FOR EACH SIMPLE BLOCK C WITH INDEX ISB, THE INDICES OF THE SURFACES FORMING THE C SET F(+) ARE STORED WITH POSITIVE SIGNS, THE INDICES OF C SURFACES FROM F(-) WITH NEGATIVE SIGNS. FOR AN EXAMPLE C REFER TO THE SAMPLE INPUT DATA FOR THE MODEL. C NCB... NUMBER OF MATERIAL COMPLEX BLOCKS IN THE MODEL. THE BLOCKS C ARE INDEXED SEQUENTIALLY BY POSITIVE INTEGERS ICB FROM 1 C TO NCB. THE FREE-SPACE BLOCKS ARE NOT INDEXED. C KCB... CONTAINS THE INDICES OF THE SIMPLE BLOCKS FORMING C INDIVIDUAL COMPLEX BLOCKS. KCB(ICB), FOR ICB = 1 TO NCB, C SPECIFY THE PARTITION OF ARRAY KSB(NCB+1:NCB+NB) AMONG THE C COMPLEX BLOCKS. HERE NB IS THE TOTAL NUMBER OF ALL C OCCURENCIES OF THE INDICES OF THE SIMPLE BLOCKS FORMING C INDIVIDUAL COMPLEX BLOCKS IN THE INPUT DATA. THE INDICES C OF THE SIMPLE BLOCKS FORMING INDIVIDUAL COMPLEX BLOCKS ARE C STORED FROM KSB(NCB+1) TO KSB(NCB+NB). THE LOCATIONS C KSB(NCB+NB+1:MCB) ARE UNDEFINED. IT MUST BE NCB+NB.LT.MCB. C THE INDICES OF THE SIMPLE BLOCKS FORMING COMPLEX BLOCK ICB C ARE STORED IN KCB(I1) TO KCB(I2), WHERE C I1 = KCB(ICB-1)+1 , C I2 = KCB(ICB) . C HERE KCB(ICB-1)=NCB FOR ICB=1. FOR AN EXAMPLE REFER TO C THE SAMPLE INPUT DATA FOR THE MODEL. C COMMON BLOCK /MODELT/ IS INCLUDED IN EXTERNAL PROCEDURE MODEL1. C COMMON BLOCK /MODELC/ IS INCLUDED IN EXTERNAL PROCEDURES MODEL1, C BLOCK, ISIDE, INTERF, AND MAY BE INCLUDED IN ANY OTHER SUBROUTINE. C ALL THE INPUT DATA ARE STORED SEQUENTIALLY IN THE SAME ORDER AS C THEY WERE READ. THE ONLY EXCEPTION ARE LOCATIONS KSB(1) TO C KSB(NSB) AND KCB(1) TO KCB(NCB) WHICH ARE INSERTED WHEN THE INPUT C DATA ARE BEING READ . THE INDEX OF THE LAST ALLOCATED NUMERIC C STORAGE UNIT OF ARRAY KSB IS NAMED MSB. THE INDEX OF THE LAST C ALLOCATED NUMERIC STORAGE UNIT OF ARRAY KCB IS NAMED MCB. THE C VALUES OF MSB AND MCB ARE GIVEN BY THE SIXTH AND SEVENTH STATEMENT C OF THE BLOCK DATA SUBROUTINE MODELB. IF THE VALUE OF MSB OR MCB C IS CHANGED, IT MUST BE ADJUSTED IN ALL SUBROUTINES WHICH INCLUDE C THE COMMON BLOCK /MODELC/. C C DATE: 1992, NOVEMBER 30 C CODED BY LUDEK KLIMES 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. FINALY, IT CALLS SUBROUTINE PARM1 (INCLUDED C IN THE SUBROUTINE FILE 'PARM.FOR') TO READ THE INPUT DATA (10) FOR THE C 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 BLOCK /MODELT/: CHARACTER*80 TEXTM COMMON/MODELT/TEXTM C COMMON BLOCK /MODELC/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB C ALL THE STORAGE LOCATIONS OF THE COMMON BLOCKS ARE DEFINED IN THIS C SUBROUTINE. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL MODELB,METR1,SRFC1,PARM1 C MODELB..BLOCK DATA SUBROUTINE OF THIS FILE. 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 ERROR MESSAGES: C 310... INSUFFICIENT MEMORY IN /MODELC/: C INSUFFICIENT MEMORY FOR THE INPUT DATA IN COMMON BLOCK C /MODELC/. THE DIMENSIONS MSB OR MCB OF ARRAYS KSB OR KCB, C RESPECTIVELY, MUST BE ENLARGED. SEE THE BLOCK DATA C SUBROUTINE MODELB. C 311... BLOCK BOUNDED BY WRONG INTERFACE: C INDEX OF THE SURFACE BOUNDING THE SIMPLE BLOCK IS GREATER C THAN THE SPECIFIED NUMBER OF SURFACES. C 312... C. BLOCK COMPOSED OF WRONG S. BLOCK: C INDEX OF A SIMPLE BLOCK COMPOSING THE COMPLEX BLOCK IS C GREATER THAN THE SPECIFIED NUMBER OF SIMPLE BLOCKS. C C DATE: 1993, DECEMBER 18 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 READ(LUN,*) I,NEXPV,NEXPQ 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 PAUSE 'ERROR 311 IN MODEL: BLOCK BOUNDED BY WRONG INTERFACE' STOP 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 STOP * 'ERROR 312 IN MODEL: C. BLOCK COMPOSED OF WRONG S. BLOCK.' 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 PAUSE 'ERROR 310 IN MODEL1: INSUFFICIENT MEMORY IN /MODELC/' STOP END 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/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB 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 SUBROUTINE BLOCK(COOR,ISRF1,ISB1,ISRF2,ISB2,ICB2,F) 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 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 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/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB 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 ERROR MESSAGES: C 313... 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. C 314... 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. C C DATE: 1993, MAY 18 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 PAUSE 'ERROR 313 IN BLOCK: WRONG INDEX OF SURFACE' STOP END IF IF(ISB1.LT.0.OR.ISB1.GT.NSB) THEN PAUSE 'ERROR 314 IN BLOCK: WRONG INDEX OF SIMPLE BLOCK' STOP END IF C C POSITION OF THE GIVEN POINT WITH RESPECT TO THE GIVEN SIMPLE BLOCK C ISB1: IF(ISB1.EQ.0) THEN C NO SIMPLE BLOCK SPECIFIED: ISRF2=0 ELSE 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=ISIDE(ISRF1,ISB1) ELSE C THE POINT IS SITUATED AT THE BOUNDARY OF SIMPLE BLOCK ISB1: ISIDE1=-ISIDE(ISRF1,ISB1) 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 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 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/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB 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 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/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB 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 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/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB 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 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 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: 1992, NOVEMBER 30 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER I REAL F,AUX1,AUX2 C IF(0.999.LT.POWER.AND.POWER.LT.1.001) THEN DO 1 I=1,N FOUT(I)=FINP(I) 1 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