C C ****************************************************************** C C P R O G R A M S E I S 8 8 C C ****************************************************************** C C NUMERICAL MODELLING OF SEISMIC WAVE FIELDS C IN 2-D LATERALLY VARYING LAYERED STRUCTURES BY THE RAY METHOD C C C 1) 2) C AUTHORS: VLASTISLAV CERVENY IVAN PSENCIK C ******* C C ADDRESSES: 1) INSTITUTE OF GEOPHYSICS C ********* CHARLES UNIVERSITY C KE KARLOVU 3 C 121 16 PRAHA 2 C CZECHOSLOVAKIA C C 2) GEOPHYSICAL INSTITUTE C CZECHOSL. ACAD. SCI. C BOCNI II C 141 31 PRAHA 4 C CZECHOSLOVAKIA C C ****************************************************************** C C DESCRIPTION OF THE PROGRAM. C *************************** C C PROGRAM SEIS88 IS DESIGNED FOR THE COMPUTATION OF RAYS OF SEIS- C MIC WAVES WHICH ARRIVE AT A SYSTEM OF RECEIVERS DISTRIBUTED REGU- C LARLY OR IRREGULARLY ALONG THE EARTH'S SURFACE, ALONG AN INTERFACE, C OR ALONG A VERTICAL PROFILE. THE GENERATION OF WAVES IS SEMIAUTOMA- C TIC. AT RECEIVERS, CORRESPONDING TRAVEL TIMES ARE AUTOMATICALLY C DETERMINED. OPTIONALLY, AMPLITUDES AND PHASE SHIFTS MAY BE EVALU- C ATED. EFFECTS OF A SLIGHT ABSORPTION MAY BE ALSO CONSIDERED. ALL C THESE QUANTITIES ARE STORED AND MAY BE OPTIONALLY PLOTTED OR USED C FOR THE COMPUTATION OF SYNTHETIC SEISMOGRAMS. C THE MODEL IS 2-D, LATERALLY INHOMOGENEOUS, WITH CURVED INTER- C FACES. INTERFACES ARE SPECIFIED BY POINTS READ FROM THE INPUT DA- C TA. THEY ARE APPROXIMATED BY CUBIC SPLINE INTERPOLATION. EACH C INTERFACE CROSSES THE WHOLE MODEL FROM ITS LEFT BORDER TO THE C RIGHT BORDER. THE INTERFACES MAY HAVE CORNER POINTS AND MAY BE C FICTITIOUS IN CERTAIN PARTS. VARIOUS INTERFACES MAY ALSO PARTIAL- C LY COINCIDE. THUS MODELS WITH VANISHING LAYERS, BLOCK STRUCTURES, C FRACTURES, ISOLATED BODIES MAY BE HANDLED BY THE PROGRAM. WITHIN C INDIVIDUAL LAYERS, THE VELOCITY MAY VARY BOTH VERTICALLY AND HO- C RIZONTALLY. C THE SOURCE MAY BE SITUATED AT ANY POINT OF THE MEDIUM, WITH C EXCEPTION OF LAYERS OF ZERO THICKNESS, I.E., IT CANNOT BE SITU- C ATED BETWEEN TWO COINCIDING INTERFACES. C ALL THE DIRECT AND PRIMARY REFLECTED WAVES P AND S, INCLUDING C THE CONVERTED WAVES AT THE POINT OF REFLECTION, CAN BE GENERATED C AUTOMATICALLY. MULTIPLE REFLECTIONS OF ARBITRARY TYPE ARE C OPTIONALLY GENERATED MANUALLY (BY INPUT DATA). THE REFRACTED C WAVES ARE CONSIDERED AS SPECIAL CASES OF REFLECTED WAVES WITH C COMPOUND RAY ELEMENTS, SEE A MORE DETAILED DESCRIPTION LATER. C THE DETERMINATION OF RAYS WHICH ARRIVE AT SPECIFIED RECEIVERS C SITUATED EITHER ON SURFACE OR ON AN INTERFACE, OR ON A VERTICAL C PROFILE - TWO-POINT RAY TRACING - IS PERFORMED BY THE MODIFIED C SHOOTING METHOD. FOR ITERATIONS TO A CONSIDERED RECEIVER, METHOD C OF HALVING OF INTERVALS, REGULA FALSI METHOD, OR THE COMBINATION C OF THESE METHODS MAY BE USED. THE ITERATIONS TO A RECEIVER CONTI- C NUE UNLESS A RAY ARRIVING TO A 'REPS' VICINITY (FOR REPS SEE C INPUT DATA CARD NO.11) OF THE RECEIVER IS FOUND. THE ARRIVAL TIME C AT THE RECEIVER IS THEN OBTAINED BY A LINEAR INTERPOLATION FROM C ARRIVAL TIMES OF RAYS TERMINATING CLOSEST TO THE RAY AND SITUATED C TO BOTH SIDES FROM THE RECEIVER. THE AMPLITUDE CORRESPONDING TO C THE RAY CLOSEST TO THE RECEIVER IS ATRIBUTED TO THE RECEIVER. C A SPECIAL CARE IS DEVOTED TO CERTAIN SINGULAR SITUATIONS IN C THE RAY FIELD. SPECIAL BRANCHES OF THE PROGRAM SEEK RAYS OF RE- C FRACTED WAVES, WHICH ARE CLOSE TO THE RAYS OF HEAD WAVES, RAYS C OF REFLECTED AND REFRACTED WAVES IN VICINITIES OF SHADOW BOUN- C DARIES, ETC. IN CASE THAT A RAY WITH IND.NE.3 (I.E. A RAY WHICH C DOES NOT COME TO THE EARTH'S SURFACE - SEE THE DESCRIPTION OF C THE PARAMETER IND) IS OBTAINED DURING THE SHOOTING, THE COR- C RESPONDING BOUNDARY RAY SEPARATING ILLUMINATED REGION (WITH C IND.EQ.3) FROM THE REGION WITH IND.NE.3, IS SOUGHT. AS SOON AS C THE BOUNDARY RAY IS FOUND THE ITERATIONS TO THE RECEIVERS IN C THE ILLUMINATED REGION ARE PERFORMED IN THE ABOVE DESCRIBED WAY. C THE PROGRAM WORKS IN A SINGLE PRECISION, SO THAT THE DETERMINA- C TION OF ALL WAVES IS IMPOSSIBLE IN CERTAIN SITUATIONS. C THE AMPLITUDES OF P, SV AND SH WAVES ARE EVALUATED BY STANDARD C RAY FORMULAE. NO MODIFICATIONS ARE APPLIED IN SINGULAR REGIONS. C THE PURE HEAD WAVES GENERALLY DO NOT EXIST IN LATERALLY INHOMO- C GENEOUS MEDIA WITH CURVED INTERFACES AND ARE NOT CONSIDERED IN C SEIS88. THEY MIGHT BE, OF COURSE,SIMULATED BY SLIGHTLY REFRACTED C WAVES. C GEOMETRICAL SPREADING IS DETERMINED BY SOLVING A SYSTEM OF TWO C LINEAR ORDINARY DIFFERENTIAL EQUATIONS OF FIRST ORDER (SO CALLED C DYNAMIC RAY TRACING SYSTEM) BY A MODIFIED EULER'S METHOD. C SLIGHT ABSORPTION MAY BE CONSIDERED. IN THIS CASE, A GLOBAL C ABSORPTION FACTOR (QUANTITY OBTAINED AS A TIME INTEGRAL OF 1/Q C ALONG A RAY, WHERE Q IS A QUALITY FACTOR) IS ALSO EVALUATED. C THE ABSORPTION IS CONSIDERED SLIGHT IF PRODUCT OF THE GLOBAL C ABSORPTION FACTOR WITH A CONSIDERED FREQUENCY IS APPROXIMATELY C EQUAL OR LESS THAN UNIT. QUALITY FACTOR DISTRIBUTION IN THE C MODEL IS RELATED TO THE DISTRIBUTION OF CORRESPONDING VELOCITY, C SEE BELOW. C C PROGRAM SEIS88 IS SUPPLEMENTED BY PROGRAMS FOR PREPARATION OF C DATA FOR VELOCITY MODEL (PROGRAM SMOOTH), TO COMPUTE SYNTHETIC C SEISMOGRAMS (PROGRAM SYNTPL), AND TO PLOT RAYS, TRAVEL-TIME AND C AMPLITUDE-DISTANCE CURVES (PROGRAM RAYPLOT), SYNTHETIC SEIS- C MOGRAMS (PROGRAM SEISPLOT) OR PARTICLE MOTION DIAGRAMS (PROGRAM C POLARPLOT). THE PLOTTING PROGRAMS REQUIRE CALCOMP ROUTINES PLOTS, C PLOT, NUMBER AND SYMBOL. C C IN THE PRESENT VERSION OF THE PROGRAM SEIS88, THERE ARE THREE C POSSIBILITIES OF THE APPROXIMATION OF THE VELOCITY DISTRIBUTION C INSIDE INDIVIDUAL LAYERS: 1) BICUBIC SPLINE INTERPOLATION (IDEN- C TICAL TO THAT USED IN SEIS81), 2) LINEAR INTERPOLATION BETWEEN C ISOVELOCITY INTERFACES, 3) PIECE-WISE BILINEAR INTERPOLATION. C APPROXIMATIONS 1) AND 3) ARE CONVENIENT FOR MODELLING SPECIFIC C FEATURES WITHIN LAYERS. THE APPROXIMATION 2) IS CONVENIENT C ESPECIALLY WHEN THERE IS NO DETAILED INFORMATION ABOUT THE C VELOCITY DISTRIBUTION WITHIN A LAYER. NOTE THAT THE APPROXIMA- C TION 2) NEEDS A MINIMUM OF INPUT DATA. ALSO NOTE THAT THIS C APPROXIMATION IS CONVENIENT FOR MODELLING SECOND ORDER INTER- C FACES (WITH CONTINUOUS VELOCITY, BUT DISCONTINUOUS GRADIENT) C AND FOR MODELLING OF MEDIA CONTAINING HOMOGENEOUS LAYERS. C THE APPROXIMATION 3) IS CONVENIENT FOR KINEMATIC MODELLING. C IT SHOULD NOT BE USED FOR THE COMPUTATION OF AMPLITUDES. IT C IS EASY TO WRITE A NEW ROUTINE MODEL AND CORRESPONDING SUB- C ROUTINES FOR ANY TYPE OF VELOCITY APPROXIMATION WITHIN LAYERS. C THE ONLY CONDITIONS ARE THAT THE ROUTINE VELOC MUST BE SPE- C CIFIED AS DESCRIBED BELOW AND IT MUS CONTAIN COMMON/AUXI/. C C THE USERS ARE KINDLY ASKED TO REFER TO PROGRAM SEIS88 IN C ANY CASE THEY PRESENT SOME RESULTS OBTAINED WITH ITS HELP. C C C C I N P U T D A T A C ******************* C C INPUT DATA FOR THE MODEL C ************************ C THE 2-D MODEL IS SPECIFIED IN A RIGHT=HANDED CARTESIAN COORDI- C NATE SYSTEM X,Y,Z. THE X- AND Y-AXES ARE HORIZONTAL, THE Z-AXIS C IS VERTICAL, POSITIVE DOWNWARDS. THE MODEL IS SITUATED IN THE C (X,Z) PLANE SO THAT X COORDINATE INCREASES FROM THE LEFT TO THE C RIGHT. THE MODEL IS BOUNDED BY TWO VERTICAL BOUNDARIES ON ITS C LEFT-HAND AND RIGHT-HAND SIDE AND BY THE FIRST AND THE LAST C INTERFACE AT THE TOP AND THE BOTTOM OF THE MODEL. THE TOP INTER- C FACE CORRESPONDS TO THE EARTH'S SURFACE. C C 1)ONE CARD C LU1,LU2,MPRINT,MTEXT,LU,MM FORMAT(3I2,17A4,2X,2I2) C LU1... THE NUMBER OF THE FILE IN WHICH RAY DIAGRAMS C TRAVEL TIMES AND AMPLITUDES ARE STORED. C LU2... THE NUMBER OF THE FILE IN WHICH DATA FOR CON- C STRUCTION OF SYNTHETIC SEISMOGRAMS ARE STORED. C NOTE: THE INPUT DATA LU1, LU2 HAVE ONLY FORMAL C MEANING HERE. THESE DATA ARE ACTUALLY SPECIFIED C IN INPUT DATA CARD NO.9. HERE THEY ARE USED ONLY C FOR CONSISTENCY OF INPUT DATA FOR SEIS88 WITH C THOSE FOR SEIS83. C MPRINT... CONTROLS THE PRINTING OF THE DESCRIPTION OF THE C MODEL ON THE LINE PRINTER. C MPRINT=0... ONLY PRINT OF INPUT DATA FOR THE MO- C DEL, NO OTHER DESCRIPTION OF THE MODEL IS PRINTED. C MPRINT=1... ONLY A SIMPLE PRINTER PLOT OF THE C VELOCITY DISTRIBUTION IN THE MODEL. C MPRINT=2... MORE DETAILED DESCRIPTION OF THE C MODEL. C MTEXT... ARBITRARY ALPHANUMERIC TEXT DESCRIBING THE MODEL. C THIS TEXT WILL APPEAR IN THE PLOTS OF SYNTHETIC C SEISMOGRAMS. C LU... THE NUMBER OF THE FILE, FROM WHICH THE INFORMA- C TION ABOUT THE VELOCITY DISTRIBUTION (DATA SUB 4) C CAN BE READ IN. IF LU=0, THE INFORMATION IS READ C FROM CARDS OR FROM TERMINAL. C MM... INDICATION OF THE TYPE OF THE ROUTINE MODEL USED C FOR THE APPROXIMATION OF VELOCITY DISTRIBUTION. C MM=0... BICUBIC SPLINE INTERPOLATION. C MM=1... LINEAR VERTICAL INTERPOLATION BETWEEN C ISOVELOCITY INTERFACES. C MM=2... PIECE-WISE BILINEAR INTERPOLATION C NOTE: IF THE VALUE OF MM DOES NOT CORRESPOND TO C THE USED ROUTINE MODEL, THE COMPUTATION TERMINATES. C C 2) ONE CARD, DESCRIPTION OF INTERFACES. C NINT,NPNT(1),NPNT(2),...,NPNT(NINT) FORMAT(16I5) C NINT... NUMBER OF INTERFACES,INCLUDING THE TOP AND THE C BOTTOM BOUNDARIES OF THE MODEL (NINT.GE.2, C NINT.LE.20). THE TOP BOUNDARY OF THE MODEL C (FREE SURFACE) AND THE BOTTOM BOUNDARY MAY C ALSO BE CURVED. C NOTE THAT THE NUMBER OF LAYERS IN THE MODEL IS C NINT-1. C NPNT(I)...NUMBER OF POINTS SPECIFYING THE I-TH INTERFACE. C (NPNT(I).GE.2,NPNT(I).LE.30). C NOTE: IN CASE OF MM=1 (SEE CARD NO.1), ALL INTERFA- C CES REPRESENT ISOLINES OF VELOCITY. THE INTERFACES C MAY BE OF THE FIRST ORDER WHEN THE VELOCITIES IMME- C DIATELY ABOVE AND BELOW INTERFACE DIFFER, OR THEY C MAY BE OF THE SECOND ORDER WHEN BOTH VELOCITIES ARE C EQUAL. C C 3) NINT SYSTEMS OF CARDS, DESCRIPTIONS OF INDIVIDUAL INTERFACES. C FOR EACH INTERFACE ONE SYSTEM OF CARDS. INTERFACES ARE TAKEN C FROM THE TOP TO THE BOTTOM. C EACH INTERFACE CROSSES THE WHOLE MODEL, FROM THE LEFT VERTI- C CAL BOUNDARY TO THE RIGHT VERTICAL BOUNDARY. C THE FIRST (TOP) INTERFACE REPRESENTS THE UPPER BOUNDARY OF THE C MODEL (THE EARTH'S SURFACE). C THE LAST (BOTTOM) INTERFACE REPRESENTS THE LOWER BOUNDARY OF C THE MODEL. C THE SYSTEM OF CARDS FOR THE I-TH INTERFACE CONSISTS OF NPNT(I) C TRIPLETS OF NUMBERS; EACH TRIPLET CORRESPONDS TO ONE POINT OF C THE INTERFACE (POINTS ARE TAKEN FROM THE LEFT TO THE RIGHT). C EACH TRIPLET CONTAINS THE FOLLOWING THREE NUMBERS: A) THE X C COORDINATE OF THE POINT, B) THE Z COORDINATE OF THE POINT, C C) A SPECIAL INDEX III, SEE BELOW. THE TRIPLETS ARE READ IN C INDEPENDENTLY FOR EACH INTERFACE, USING FORMAT(3(2F10.5,I5)). C THE FIRST POINTS OF ALL INTERFACES MUST LIE ON ONE VERTICAL C LINE, WHICH IS CALLED THE LEFT BOUNDARY OF THE MODEL. C THE LAST POINTS OF ALL INTERFACES MUST LIE ON ANOTHER VERTI- C CAL LINE, WHICH IS CALLED THE RIGHT BOUNDARY OF THE MODEL. C THE PART OF THE INTERFACE BETWEEN TWO SUCCESSIVE POINTS ON C THE INTERFACE IS CALLED THE ELEMENT OF THE INTERFACE. C C THE QUANTITY III CONTROLS THE CHARACTER OF THE INTERFACE IN C THE NEIGHBOURHOOD OF THE CORRESPONDING POINT: C III=0 THE INTERFACE ITSELF, AS WELL AS ITS FIRST AND C SECOND DERIVATIVES, ARE SMOOTH AT THE CORRESPON- C DING POINT. C III=-1 THE CORNER POINT IN THE INTERFACE (DISCONTINUOUS C FIRST DERIVATIVE). C III=-2 THE CORNER POINT OF THE INTERFACE. MOREOVER, THE C INTERFACE TO THE RIGHT SIDE OF THE POINT IS C FICTITIOUS. C III=K,K.GT.0 THE ELEMENT OF THE INTERFACE TO THE RIGHT C OF THE CORRESPONDING POINT COINCIDES WITH THE C K-TH ELEMENT OF THE ADJOINING UPPER INTERFACE, C STARTING FROM THE K-TH POINT OF THE UPPER INTER- C FACE. THE COORDINATES OF THE CORRESPONDING POINTS C AT BOTH INTERFACES MUST BE THE SAME. C C C 4) DESCRIPTION OF VELOCITY DISTRIBUTION IN INDIVIDUAL LAYERS. C NOTE THAT THESE INPUT DATA MAY DIFFER FOR DIFFERENT KINDS C OF APPROXIMATION OF THE VELOCITY DISTRIBUTION. C C *************************************************** C MM= BICUBIC SPLINE OR PIECE-WISE BILINEAR INTERPOLATION C 0,2 (MM=0 OR MM=2, SEE INPUT DATA CARD NO.1) C *************************************************** C C (NINT-1) SYSTEMS OF CARDS, DESCRIPTIONS OF VELOCITY DISTRIBU- C TION IN INDIVIDUAL LAYERS. LAYERS ARE TAKEN FROM THE TOP TO C THE BOTTOM. STARTING FROM THE FIRST LAYER. THE DATA CAN BE C READ IN EITHER FROM A CARD READER, TERMINAL, ETC., OR FROM C A DATA SET STORED IN A FORMATTED FORM IN THE FILE LU (SEE DA- C TA SUB 1). VELOCITY IS SPECIFIED AT GRID POINTS OF A RECTANGU- C LAR NETWORK. THE NETWORK MAY BE DIFFERENT IN DIFFERENT LAYERS. C THE NETWORK MUST ALWAYS COVER THE WHOLE CORRESPONDING LAYER. C NOTE THAT IN THIS CASE THE INPUT DATA ARE IDENTICAL TO THOSE C USED IN THE PROGRAM SEIS81. C THE SYSTEM OF CARDS FOR EACH LAYER CONSISTS OF THE FOLLOWING C CARDS: C C 4A) ONE CARD,NUMBER OF NETWORK LINES. C MX,MY FORMAT(2I5) C MX... NUMBER OF VERTICAL LINES IN THE NETWORK. C MY... NUMBER OF HORIZONTAL LINES IN THE NETWORK. C C 4B) ONE CARD (OR GROUP OF CARDS), X COORDINATES OF VERTICAL NET- C WORK LINES. C SX(1),SX(2),...,SX(MX) FORMAT(8F10.5) C SX(I)... X COORDINATE OF THE I-TH VERTICAL LINE. C SX(1) SPECIFIES THE LEFT-HAND BOUNDARY C OF THE MODEL (THE SAME FOR ALL LAYERS). C SX(MX) SPECIFIES THE RIGHT-HAND BOUNDARY C OF THE MODEL (THE SAME FOR ALL LAYERS). C C 4C) ONE CARD (OR GROUP OF CARDS), Z COORDINATES OF HORIZONTAL C NETWORK LINES. C SY(1),SY(2),...,SY(MY) FORMAT(8F10.5) C SY(I)... Z COORDINATE OF THE I-TH HORIZONTAL LINE. C SY(1) AND SY(MY) MUST BE CHOSEN IN SUCH A WAY C THAT THE NETWORK COVERS THE WHOLE LAYER. THEY C MAY BE COMPLETELY OUTSIDE THE LAYER. C C 4D) A GROUP OF CARDS, VELOCITY VALUES AT GRID POINTS ALONG THE C VERTICAL NETWORK LINES. FOR EACH VERTICAL NETWORK LINE,'MY' C VELOCITY VALUES ARE READ IN, INDEPENDENTLY FOR EACH LINE, C FORMAT(8F10.5). THE GRID POINTS ALONG THE VERTICAL NETWORK C LINE ARE TAKEN FROM THE TOP TO THE BOTTOM. THE VERTICAL NET- C WORK LINES ARE TAKEN FROM LEFT TO RIGHT. C C NOTES ON MAXIMUM PERMITTED DIMENSIONS: C THE SUM OF ALL VERTICAL LINES IN ALL LAYERS MAY NOT EXCEED C 350. C THE SUM OF ALL HORIZONTAL LINES IN ALL LAYERS MAY NOT EXCEED C 350. C THE TOTAL NUMBER OF ALL GRID POINTS IN ALL LAYERS (I.E., THE C TOTAL NUMBER OF ALL VELOCITY VALUES READ IN) MAY NOT EXCEED C 1000. C C ************************************************************* C C C *************************************************** C MM= LINEAR INTERPOLATION BETWEEN ISOVELOCITY INTERFACES C 1 (MM=1, SEE INPUT DATA CARD NO.1) C *************************************************** C C THESE DATA ARE ALWAYS READ IN FROM A CARD READER OR TERMINAL. C C C 4A) ONE CARD (OR A GROUP OF CARDS), VELOCITIES BELOW INTERFACES C V1(1),V1(2),....,V1(NINT-1) FORMAT(8F10.5) C V1(I)... VELOCITY IMMEDIATELY BELOW I-TH INTERFACE, C I.E., ALONG THE UPPER BOUNDARY OF THE I-TH C LAYER. C C 4B) ONE CARD (OR A GROUP OF CARDS), VELOCITIES ABOVE INTERFACES C V2(1),V2(2),...,V2(NINT-1) FORMAT(8F10.5) C V2(I)... VELOCITY IMMEDIATELY ABOVE (I+1)-TH INTERFACE, C I.E., ALONG THE LOWER BOUNDARY OF THE I-TH C LAYER. C C C 5) ONE CARD, SPECIFICATION OF DENSITIES, S VELOCITIES AND QUALITY C FACTORS. C NRO, NVS(1), NVS(2),...,NVS(NINT-1), NABS FORMAT(16I5)) C NRO... NRO=0... THE DENSITY AT ANY POINT OF THE MEDIUM C IS DETERMINED AUTOMATICALLY FROM THE RELATION C DENSITY=1.7+0.2*VP (WHERE VP IS THE P VELOCITY). C THE FOLLOWING CARD NO.6A IS NOT READ IN IN THIS C CASE. C NRO=1... THE DENSITY IN THE I-TH LAYER IS DETER- C MINED FROM THE RELATION C DENSITY=RHO1(I)+RHO2(I)*VP, C WHERE RHO1(I) AND RHO2(I) ARE READ IN IN THE INPUT C DATA CARD NO.6A. C FOR A LIQUID LAYER, DENSITY=1 AUTOMATICALLY. C NVS(I)... NVS(I)=0...THE VELOCITY VALUES READ IN SUB 4 C CORRESPOND TO P WAVE VELOCITY DISTRIBUTION C IN THE I-TH LAYER. S WAVE VELOCITIES IN THE I-TH C LAYER ARE DETERMINED FROM THE RELATION C VS=VP/PTOS(I). THE VALUES OF PTOS(I) ARE READ IN C IN THE INPUT CARD NO.7 C NVS(I)=1... THE VELOCITY VALUES READ IN SUB 4 C CORRESPOND TO THE S WAVE VELOCITY DISTRIBUTION IN C THE I-TH LAYER. P WAVE VELOCITIES IN THIS LAYER C ARE DETERMINED FROM THE RELATION C VP=VS*PTOS(I), WHERE PTOS(I) ARE READ IN IN THE C CARD NO.7. C NOTE THAT THE FIRST LAYER MAY BE ALSO LIQUID C (E.G. OCEAN). IN THIS CASE, ONLY NVS(1)=0 IS ALLO- C WED, NVS(1)=1 IS NOT PERMITTED. PUT PTOS(1).GE.100 C IN THIS CASE, THEN THE S WAVE VELOCITY IN THE C FIRST LAYER IS AUTOMATICALLY ZERO. C NABS... NABS=0... NO ABSORPTION IS CONSIDERED. C NABS=1... SLIGHT ABSORPTION IS CONSIDERED. C FOR NABS=1, THE DATA SUB 6B ARE READ IN. THE C QUALITY FACTORS QP AND QS ARE DETERMINED FROM C ONE OF THE RELATIONS GIVEN IN DESCRIPTION OF C INPUT DATA SUB 6B. NO SPECIFICATION OF THE TYPE C OF THE ABSORPTION (FUTTERMAN'S OR NON-CAUSAL) IS C MADE IN THIS PROGRAM. THE TYPE IS SPECIFIED LATER C IN PROGRAMS RAYPLOT OR SYNTPL, INCLUDED IN THIS C PACKAGE. HERE ONLY THE GLOBAL ABSORPTION FACTOR C IS COMPUTED. C C C 6A)ONE CARD, SPECIFICATION OF DENSITIES. C NOTE: THIS CARD IS INCLUDED ONLY WHEN NRO=1. C RHO1(1),RHO2(1),RHO1(2),RHO2(2),..,RHO1(NINT-1),RHO2(NINT-1) C FORMAT(8F10.5) C RHO1(I),RHO2(I)... PARAMETERS OF THE RELATION FOR THE C DETERMINATION OF THE DENSITY FROM P WAVE VELO- C CITIES IN THE I-TH LAYER, C DENSITY=RHO1(I)+RHO2(I)*VP. C C C 6B)(NINT-1) CARDS, SPECIFICATION OF QUALITY FACTORS. C NOTE: THESE CARDS ARE INCLUDED ONLY IF NABS.NE.0. C ONE CARD FOR EACH LAYER, FROM TOP TO THE BOTTOM. C EACH CARD HAS A FORM: C NQP,NQS,QP1,QP2,QP3,QS1,QS2,QS3 FORMAT(2I5,6F10.5) C NQP... NQP=0... THE QUALITY FACTOR QP IN THE CONSIDERED C LAYER IS DETERMINED FROM THE RELATION: C QP=QP1+QP2*VP+QP3*VP*VP, C WHERE VP IS THE P WAVE VELOCITY; C NQP=1... THE QUALITY FACTOR QP IN THE CONSIDERED C LAYER IS DETERMINED FROM THE RELATION: C QP=QP1+QP2*VP+QP3/VP. C NQS... THE SAME AS FOR QP; IT IS ONLY NECESSARY TO RE- C PLACE THE QUANTITIES VP,QP1,QP2,QP3 BY VS,QS1, C QS2,QS3, WHERE VS IS THE S WAVE VELOCITY. C QP1,QP2,QP3,QS1,QS2,QS3... COEFFICIENTS OF THE RELA- C TIONS FOR THE DETERMINATION OF QUALITY FACTORS C IN A CONSIDERED LAYER. C C C 7) ONE CARD, SPECIFICATION OF THE RATIOS BETWEEN S AND P VELO- C CITIES IN INDIVIDUAL LAYERS C PTOS(1),PTOS(2),...,PTOS(NINT-1) FORMAT(8F10.5) C PTOS(I)...RATIO OF THE P WAVE VELOCITY TO THE S WAVE VELO- C CITY IN THE I-TH LAYER,VS=VP/PTOS(I). C NOTE THAT PTOS(I) MUST BE ALWAYS GREATER THAN C SQRT(2.). FOR PTOS(I).LT.SQRT(2.), E.G. WHEN C PTOS(I) IS NOT SPECIFIED, PTOS(I)=1.732 AUTO- C MATICALLY. C FOR PTOS(I).GE.100., VS VELOCITY IN THE I-TH C LAYER IS ZERO AUTOMATICALLY (LIQUID LAYER). C C C 8) ONE CARD, CONTROLS THE PRINTER PLOT OF THE P VELOCITY DISTRI- C BUTION IN THE MODEL, WHICH IS PRINTED ON A LINE PRINTER WHEN C MPRINT.GE.1 (SEE CARD NO.1). THE PRINTER PLOT CONSISTS OF A C SYSTEM OF DIGITS FROM 0 TO 9 AND ASTERISKS. THE DIGITS ARE C DETERMINED FROM THE COMPUTED VELOCITY DISTRIBUTION USING C THE RELATION I=IFIX(10.*(V-VMIN)/(VMAX-VMIN)). OUTSIDE THE C RANGE VMIN, VMAX, THE ASTERISKS ARE PRINTED. THE CARD SHOULD C BE FORMALLY INCLUDED EVEN WHEN MPRINT=0. C VMIN,VMAX,BMIN,BMAX,BLEFT,BRIGHT FORMAT(8F10.5) C VMIN,VMAX... MINIMUM AND MAXIMUM P VELOCITY SHOWN IN THE C PRINTER PLOT C BMIN,BMAX... Z COORDINATES OF THE TOP AND THE BOTTOM HORI- C ZONTAL BOUNDARIES OF THE PRINTER PLOT. C WHEN BMIN AND BMAX ARE NOT SPECIFIED,Z-COORDINATES C OF THE FIRST POINTS OF THE TOP AND THE BOTTOM C INTERFACES ARE TAKEN AS BMIN AND BMAX. C BLEFT,BRIGHT... X COORDINATES OF THE VERTICAL BOUNDARIES C OF THE PRINTER PLOT. THEY MAY DIFFER FROM VERTI- C CAL BOUNDARIES OF THE MODEL. IF BLEFT AND BRIGHT C ARE NOT SPECIFIED, THEN THEY ARE TAKEN AS VERTICAL C BOUNDARIES OF THE MODEL. C C C C SPECIFICATION OF PARAMETERS OF COMPUTED SEISMIC BODY WAVES C C C ***************************************************************** C C 9) ONE CARD, VARIOUS SWITCHES. C ICONT,MEP,MOUT,MDIM,METHOD,IBP,IBS,IUP,IUS,IDP,IDS,IREAD,MREG, C ITMAX,NLAY,LU1,LU2,LTAPE,KSH,ILOC FORMAT(26I3) C ICONT... CONTROLS THE CONTINUATION OF THE COMPUTATION. C ICONT=0... TERMINATION OF COMPUTATION. THIS IS C THE LAST CARD IN THE INPUT DATA CARD SET. C ICONT=1... THE COMPUTATIONS CONTINUE, THE NEXT C CARD SHOULD BE CARD NO 10. C ICONT=-1... THE COMPUTATIONS CONTINUE FOR A NEW C MODEL OF MEDIUM.THE NEXT CARD SHOULD BE THE CARD C NO.1. C IN CASE OF ICONT=0 OR ICONT=-1, THE REMAINING C QUANTITIES IN THIS CARD MAY BE ARBITRARY. C MEP... ABS(MEP)... NUMBER OF RECEIVER POSITIONS. THE RE- C CEIVERS CAN BE SITUATED EITHER ALONG THE EARTH'S C SURFACE OR ALONG ONE OF THE INTERFACES, OR ALONG C A VERTICAL PROFILE. IN THE FORMER TWO CASES, THE C RECEIVERS ARE SPECIFIED BY A VARYING X COORDINATE, C Z COORDINATES FOLLOW AUTOMATICALLY FROM THE EQUA- C TION OF THE SURFACE OR INTERFACE. IN THE LATTER C CASE, THE RECEIVERS ARE SPECIFIED BY A VARYING Z C COORDINATE, X COORDINATES FOLLOW AUTOMATICALLY C FROM THE SPECIFICATION OF THE VERTICAL PROFILE. C THE SIGN OF MEP CONTROLS THE WAY IN WHICH THE C SYSTEM OF RECEIVER POSITIONS IS SPECIFIED (SEE C CARD NO.10). C IF(MEP.GT.0),RECEIVERS ARE DISTRIBUTED REGULARLY. C ONLY THE POSITIONS OF THE FIRST RECEIVER AND THE C STEP IN THE VARYING COORDINATE ARE READ IN, SEE C CARD NO.10. C IF (MEP.LT.0), THE RECEIVERS ARE DISTRIBUTED IRRE- C GULARLY ALONG THE PROFILE. THE VARYING COORDINATES C OF ALL RECEIVER POSITIONS ARE READ IN, SEE CARD C NO.10. C ABS(MEP).LE.100. C MOUT... CONTROLS THE PRINTING OF THE OUTPUT ON A LINE C PRINTER. C MOUT=0... ONLY INPUT DATA ARE REPRODUCED ON A C LINE PRINTER AND LISTING OF WAVE CODES IS PRINTED. C MOUT=1... ONLY THE TABLE OF THE TERMINATION POINTS C OF RAYS WHICH WERE SUCCESFULLY ITERATED AND TER- C MINATE AT SPECIFIED RECEIVERS AND OF BOUNDARY RAYS C MOUT=2... THE TABLE OF TERMINATION POINTS OF ALL C RAYS COMPUTED IN THE PROCESS OF ITERATIONS. C MDIM... CONTROLS ON WHICH LEVEL AND WITH WHICH TYPE OF C SOURCE THE PROGRAM SHOULD WORK. C MDIM=0... ONLY THE RAYS AND TRAVEL TIMES ARE COM- C PUTED, NO AMPLITUDES AND PHASE SHIFTS. C FOR MDIM.GT.0, THE AMPLITUDES AND PHASE SHIFTS C ARE COMPUTED. C THE SPREADING MAY BE OPTIONALLY COMPUTED IN TWO C DIFFERENT WAYS: C MDIM=2... LINE SOURCE APPROXIMATION C MDIM=3... POINT SOURCE APPROXIMATION. C FOR MDIM=1, THE GEOMETRICAL SPREADING IS NOT CON- C SIDERED IN THE EVALUATION OF RAY AMPLITUDES. C METHOD... TO FIND THE RAYS WHICH HAVE TERMINATION POINTS C AT A SPECIFIED SYSTEM OF RECEIVERS, THE SHOOTING C METHOD IS USED. TO PERFORM ITERATIONS, VARIOUS C METHODS CAN BE APPLIED: C METHOD=0... REGULA FALSI. C METHOD=1... COMBINATION OF REGULA FALSI METHOD C AND METHOD OF HALVING INTERVALS. C METHOD=2... HALVING OF INTERVALS. C METHOD=3... METHOD BASED ON PARAXIAL RAY APPRO- C XIMATION. C IBP,IBS,IUP,IUS,IDP,IDS... SWITCHES WHICH CONTROL THE AUTO- C MATIC GENERATION OF NUMERICAL CODES OF ELEMENTARY C WAVES. ONLY DIRECT WAVES AND PRIMARY REFLECTED C WAVES (POSSIBLY CONVERTED AT THE POINT OF REFLEC- C TION) CAN BE GENERATED AUTOMATICALLY. C IBP... CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF- C LECTED WAVES WHICH START AS A P WAVE FROM THE C SOURCE. THE REFLECTION TAKES PLACE BELOW THE C SOURCE. C IBP=0... NO PRIMARY REFLECTED WAVES ARE GENERATED. C IBP=1... PP PRIMARY REFLECTED WAVES ARE GENERATED. C IBP=2... PS PRIMARY REFLECTED WAVES ARE ALSO GENE- C RATED (TOGETHER WITH PP PRIMARY REFLECTED WAVES). C IBS... CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF- C LECTED WAVES WHICH START AS AN S WAVE FROM THE C SOURCE. THE REFLECTION TAKES PLACE BELOW THE SOUR- C CE. C IBS=0... NO PRIMARY REFLECTED WAVES ARE GENERATED. C IBS=1... SS PRIMARY REFLECTED WAVES ARE GENERATED. C IBS=2... SP PRIMARY REFLECTED WAVES ARE ALSO GENE- C RATED (TOGETHER WITH SS PRIMARY REFLECTED WAVES). C IUP... CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF- C LECTED WAVES WHICH START AS A P WAVE FROM THE C SOURCE. THE REFLECTION TAKES PLACE ABOVE THE SOUR- C CE. C IUP=0... NO PRIMARY REFLECTED WAVES ARE GENERATED. C IUP=1... PP PRIMARY REFLECTED WAVES ARE GENERATED. C IUP=2... PS PRIMARY REFLECTED WAVES ARE ALSO GENE- C RATED (TOGETHER WITH PP PRIMARY REFLECTED WAVES). C IUS... CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF- C LECTED WAVES WHICH START AS AN S WAVE FROM THE C SOURCE. THE REFLECTION TAKES PLACE ABOVE THE SOUR- C CE. C IUS=1... NO PRIMARY REFLECTED WAVES ARE GENERATED. C IUS=2... SS PRIMARY REFLECTED WAVES ARE GENERATED. C IUS=3... SP PRIMARY REFLECTED WAVES ARE ALSO GENE- C RATED (TOGETHER WITH SS PRIMARY REFLECTED WAVES). C IDP... CONTROLS THE AUTOMATIC GENERATION OF THE DIRECT C P WAVE. C IDP=0... NO DIRECT P WAVE IS GENERATED. C IDP=1... DIRECT P WAVE IS GENERATED. C IDS... CONTROLS THE AUTOMATIC GENERATION OF THE DIRECT C S WAVE. C IDS=0... NO DIRECT S WAVE IS GENERATED C IDS=1... DIRECT S WAVE IS GENERATED. C IREAD... CONTROLS THE MANUAL GENERATION OF NUMERICAL CODES C OF ELEMENTARY WAVES. C IREAD=0... NO NUMERICAL CODES ARE MANUALLY GENE- C RATED. C IREAD=1... NUMERICAL CODES OF CERTAIN ELEMENTARY C WAVES ARE MANUALLY GENERATED, SEE CARD NO. 13. C MREG... CONTROLS THE EVALUATION OF AMPLITUDES AT THE C TERMINATION POINT OF THE RAY. C MREG=0.. AT THE EARTH'S SURFACE, THE COEFFICIENTS C OF CONVERSION ARE APPLIED (THEY TAKE INTO ACCOUNT C THE EFFECTS OF THE EARTH'S SURFACE). C MREG=1... AT THE EARTH'S SURFACE, THE COEFFICIENTS C OF CONVERSION ARE NOT APPLIED, THE COMPONENTS OF C THE DISPLACEMENT VECTOR ARE CALCULATED AS IF WITHIN C AN INFINITE HALFSPACE. C NOTE: WHEN THE TERMINATION POINT OF THE RAY IS C SITUATED INSIDE THE MEDIUM (I.E. WHEN IND.NE.3), C THE HORIZONTAL AND VERTICAL COMPONENTS OF THE DIS- C PLACEMENT ARE ALWAYS COMPUTED AS IF WITHIN AN C INFINITE HALFSPACE (EVEN IF THE RECEIVER IS SITU- C ATED AT AN INTERFACE). C SIMILARLY, WHEN THE FIRST LAYER IS LIQUID (OCEAN), C THE COEFFICIENTS OF CONVERSION ARE NOT USED. C ITMAX... NUMBER OF ITERATIONS PERMITTED IN A 2-POINT RAY C TRACING. IF ITMAX=0, THEN ITMAX=20. C NLAY... CONTROLS THE GENERATION OF THE NUMERICAL CODE OF C PRIMARY REFLECTED WAVES FROM THE LAST INTERFACE. C NLAY=0... THE PRIMARY REFLECTED WAVES FROM THE C LAST INTERFACE ARE GENERATED (WHEN IBP.NE.0 OR C IBS.NE.0). C NLAY=1... THEY ARE NOT GENERATED. C LU1... THE NUMBER OF THE FILE IN WHICH RAY DIAGRAMS, C TRAVEL TIMES AND AMPLITUDES ARE STORED. C FOR LU1=0, THESE QUANTITIES ARE NOT STORED. C LU2. THE NUMBER OF THE FILE IN WHICH DATA FOR THE C CONSTRUCTION OF RAY SYNTHETIC SEISMOGRAMS ARE C STORED. C FOR LU2=0, THESE QUANTITIES ARE NOT STORED. C LTAPE... SWITCH CONTROLLING THE FORM, IN WHICH THE RESULTS C OF COMPUTATIONS ARE STORED IN THE FILE LU2. C LTAPE=0... THE DATA ARE STORED IN A FORMATTED C FORM. C LTAPE=1... THE DATA ARE STORED IN A BINARY FORM. C KSH... SWITCH CONTROLLING SELECTION OF THE TYPE OF COM- C PUTED WAVES. C KSH=0... ONLY P AND SV WAVES ARE CONSIDERED. NO C SH WAVES ARE COMPUTED. C KSH=1... ONLY SH WAVES ARE CONSIDERED. NO P, SV C WAVES ARE COMPUTED. ANY WAVE CONTAINING AN ELEMENT C ALONG WHICH THE WAVE PROPAGATES AS A P WAVE IS C AUTOMATICALLY EXCLUDED. C KSH=2... P, SV AND SH WAVES ARE CONSIDERED TOGE- C THER. FOR ANY PURELY SHEAR WAVE SPECIFIED BY THE C CODE, BOTH SV AND SH AMPLITUDES ARE COMPUTED. C ILOC... SWITCH SPECIFYING THE LOCATION OF RECEIVERS. C ILOC=0... THE RECEIVERS ARE LOCATED AT THE EARTH'S C SURFACE. C ILOC=1... THE RECEIVERS ARE LOCATED ON A VERTICAL C PROFILE. C ILOC.GT.100... THE RECEIVERS ARE LOCATED ON THE C (ILOC-100)TH INTERFACE. C C 10) ONE CARD OR GROUP OF CARDS: SPECIFICATION OF RECEIVER C POSITIONS. C C WHEN THE RECEIVERS ARE SITUATED ALONG THE EARTH'S SURFACE OR C AN INTERFACE (ILOC.NE.1), THEY ARE SPECIFIED BY THEIR X-CO- C ORDINATES. C WHEN THE RECEIVERS ARE SITUATED ALONG A VERTICAL PROFILE C (ILOC.EQ.1), THEY ARE SPECIFIED BY THEIR Z-COORDINATES. C C WHEN MEP.LT.0: (DENOTE NDST=-MEP) C DST(1),DST(2),...,DST(NDST),XVSP FORMAT(8F10.5) C THE RECEIVERS ARE DISTRIBUTED IRREGULARLY ALONG C THE EARTH'S SURFACE, AN INTERFACE, OR A VERTICAL C PROFILE. C DST(I) DENOTES THE VARYING COORDINATE OF THE I-TH C RECEIVER. DST(I) SHOULD FORM A MONOTONOUSLY IN- C CREASING (DECREASING) SERIES. C XVSP... THE X-COORDINATE OF THE VERTICAL PROFILE C FOR ILOC.EQ.1. OTHERWISE, DUMMY PARAMETER. C C WHEN MEP.GT.0: C RMIN,RSTEP,XVSP FORMAT(8F10.5) C THE RECEIVERS ARE DISTRIBUTED REGULARLY ALONG C THE EARTH'S SURFACE, AN INTERFACE, OR A VERTICAL C PROFILE. C THE VARYING COORDINATE OF THE I-TH RECEIVER IS C GIVEN BY THE FORMULA DST(I)=RMIN+RSTEP*(I-1). C XVSP... THE X-COORDINATE OF THE VERTICAL PROFILE C FOR ILOC.EQ.1. OTHERWISE, DUMMY PARAMETER. C C 11) ONE CARD, POSITION OF THE SOURCE,ETC. C XSOUR,ZSOUR,TSOUR,REPS,REPS1,BRD(1),BRD(2) FORMAT(8F10.5) C XSOUR... THE X-COORDINATE OF THE SOURCE C ZSOUR... THE Z-COORDINATE OF THE SOURCE C THE SOURCE MAY BE SITUATED AT ANY POINT OF THE C MODEL, INCLUDING ITS BOUNDARIES. THE SOURCE C SITUATED ON AN INTERFACE (WITH EXCEPTION OF THE C LAST INTERFACE) IS CONSIDERED AS A POINT OF THE C LAYER LYING UNDER THE MENTIONED INTERFACE. IN CASE C OF LAST INTERFACE, THE SOURCE IS CONSIDERED AS C A POINT OF THE LAST LAYER. THE SOURCE CANNOT BE C SITUATED IN THE LAYER OF ZERO THICKNESS (IN CASE C OF COINCIDENCE OF NEIGHBOURHOOD INTERFACES). C IF ABS(ZSOUR).LT.0.0001, THE SOURCE IS CONSIDERED C AS LYING ON THE EARTH SURFACE AND ITS REAL C Z COORDINATE IS DETERMINED AUTOMATICALLY (THIS IS C VERY IMPORTANT IN CASE WHEN THE EARTH'S SURFACE IS C CURVED AND ACTUAL Z-COORDINATE OF THE SOURCE IS NOT C A PRIORI KNOWN). C TSOUR... INITIAL (HYPOCENTRAL) TIME C REPS... RADIUS OF THE VICINITY OF A RECEIVER. IF A RAY C HAS ITS TERMINATION POINT WITHIN THIS VICINITY, C THE ITERATIONS (2-POINT RAY TRACING) TERMINATE. C IF REPS IS NOT SPECIFIED, THEN REPS=0.05. C REPS1... CONTROLS THE TERMINATION OF ITERATIONS TO BOUNDARY C RAYS> IF THE DISTANCE BETWEEN TERMINATION POINTS C OF TWO SUCCESSIVELY GENERATED RAYS IS LESS THAN C REPS1, THE ITERATIONS TERMINATE. IF REPS1 IS NOT C SPECIFIED, THEN REPS1=0.05*REPS. C BRD(1),BRD(2)... X-COORDINATES OF THE LEFT AND RIGHT VER- C TICAL BOUNDARIES OF THE PART OF THE MODEL - CALLED C FURTHER ON SUBMODEL - IN WHICH THE COMPUTATIONS C FOR THE GIVEN SOURCE WILL BE PERFORMED. THESE C BOUNDARIES MAY GENERALLY DIFFER FROM VERTICAL C BOUNDARIES OF THE WHOLE MODEL. THEY MUST BE, HOW- C EVER, SITUATED WITHIN IT. IF ANY OF THE BOUN- C DARIES OF THE SUBMODEL IS SPECIFIED OUTSIDE THE C WHOLE MODEL, THEN IT IS AUTOMATICALLY CORRECTED C AND TAKEN AS THE CORRESPONDING BOUNDARY OF THE C MODEL. C IF BRD(1) AND BRD(2) ARE NOT SPECIFIED, THEN THEY C ARE AUTOMATICALLY CONSIDERED AS VERTICAL BOUNDARIES C OF THE WHOLE MODEL. C C 12) ONE CARD, QUANTITIES THAT CONTROL THE BASIC SYSTEM OF INITIAL C ANGLES IN THE TWO-POINT RAY TRACING, ETC. C DT,AMIN1,ASTEP1,AMAX1,AMIN2,ASTEP2,AMAX2,AC FORMAT(8F10.5) C DT... TIME STEP IN THE INTEGRATION OF THE RAY-TRACING C SYSTEM. DT SHOULD BE GREATER THAN ZERO. C IF DT.LT.0.00001, THEN DT=1. C AMIN1,ASTEP1,AMAX1... DETERMINE THE SYSTEM OF INITIAL C ANGLES FOR PRIMARY REFLECTED WAVES (AND POSSIBLY C FOR OTHER MANUALLY GENERATED ELEMENTARY WAVES,THE C FIRST ELEMENT OF WHICH STRIKES THE INTERFACE SITU- C ATED BELOW THE SOURCE). SPECIFIED IN RADIANS. THE C ANGLES ARE MEASURED IN THE PLANE (X,Z) FROM THE PO- C SITIVE DIRECTION OF THE X-AXIS, POSITIVE CLOCKWISE. C AMIN2,ASTEP2,AMAX2... DETERMINE THE SYSTEM OF INITIAL ANGLES C FOR THE DIRECT WAVES (AND POSSIBLY FOR OTHER MA- C NUALLY GENERATED ELEMENTARY WAVES,THE FIRST ELE- C MENT OF WHICH STRIKES THE INTERFACE ABOVE THE C SOURCE). SPECIFIED IN RADIANS. THE ANGLES ARE MEA- C SURED IN THE PLANE (X,Z) FROM THE POSITIVE DIRECTION C OF THE X-AXIS< POSITIVE CLOCKWISE. C AC... REQUIRED ACCURACY OF INTEGRATION OF THE RAY TRA- C CING SYSTEM. RECOMENDED VALUES: 0.0001 - 0.001. C IF AC NOT SPECIFIED, THEN AC=0.0001. C C 13) ARBITRARY NUMBER OF CARDS, ENDING BY A BLANK CARD. C MANUAL GENERATION OF NUMERICAL CODES OF ELEMENTARY WAVES. C THE CARDS ARE INCLUDED ONLY WHEN IREAD.NE.0, SEE CARD NO.9 C FOR IREAD. THE LAST CARD OF THE SYSTEM IS A BLANK CARD. C KC,KREF,CODE(1),CODE(2),...,CODE(KREF) FORMAT(24I3) C KC,KREF,CODE(1),CODE(2),...,CODE(KREF) THE NUMERICAL CODE C OF THE WAVE. THE WHOLE RAY IS DIVIDED INTO ELE- C MENTS, EACH OF WHICH LIES BETWEEN TWO SUCCESSIVE C POINTS AT WHICH THE RAY STRIKES AN INTERFACE. C IF THE END POINTS OF AN ELEMENT LIE ON DIFFE- C RENT INTERFACES, IT IS CALLED A SIMPLE ELEMENT, C WHEN THE END POINTS LIE ON THE SAME INTERFACE, C IT IS CALLED A COMPOUND ELEMENT. ANY COMPOUND ELE- C MENT IS FORMALLY REGARDED AS TWO SIMPLE ELEMENTS. C THE PART OF A RAY BETWEEN THE SOURCE AND THE FIRST C POINT AT WHICH THE RAY STRIKES AN INTERFACE IS C CONSIDERED AS ONE ELEMENT OF THE RAY. IF THE SOUR- C CE IS SITUATED ON AN INTERFACE AND A RAY FROM THE C SOURCE HITS THE SAME INTERFACE, THE ELEMENT OF C THE RAY IS ALSO COMPOUND ELEMENT. C WHEN THE RAY CROSSES A LAYER OF ZERO THICKNESS C (COINCIDING INTERFACES), THE RAY ELEMENT WITHIN C THE LAYER IS FORMALLY CONSIDERED IN THE SAME WAY C AS IN THE CASE OF A LAYER OF A NON-ZERO THICKNESS. C KC... KC.NE.0... THE CODE CORRESPONDS TO A REFLECTED, OR C MULTIPLY REFLECTED, POSSIBLY CONVERTED WAVE. C KC=1... AFTER LEAVING THE SOURCE THE RAY SHOULD C STRIKE THE INTERFACE BELOW THE SOURCE. C KC=-1... AFTER LEAVING THE SOURCE THE RAY SHOULD C STRIKE THE INTERFACE ABOVE THE SOURCE. C KC=0... DIRECT WAVE (NO REFLECTIONS AND NO C CONVERSIONS AT INTERFACES, ONLY TRANSMISSIONS). C KREF... THE NUMBER OF SIMPLE ELEMENTS OF THE RAY. C IN CASE OF REFRACTED WAVE KREF=1. C CODE(I)...ABS(CODE(I))... THE NUMBER OF THE LAYER IN WHICH C THE I-TH SIMPLE ELEMENT OF THE RAY IS SITUATED. C CODE(I).GT.0 FOR A P-WAVE ELEMENT, C CODE(I).LT.0 FOR A S-WAVE ELEMENT. C COMPOUND ELEMENTS ARE REGARDED AS TWO SIMPLE C ELEMENTS AND ARE DENOTED BY TWO EQUAL NUMBERS. C NOTE: CODE(I) IS AN INTEGER. C C C ****************************************************************** C C C THE TERMINATION OF COMPUTATIONS. C ******************************** C C THE SYSTEM OF CARDS NO. 9 - 13 CAN BE REPEATED ANY NUMBER OF C TIMES TO PERFORM VARIOUS COMPUTATIONS FOR THE SAME MODEL. C EACH SYSTEM OF CARDS NO.9-13 PRODUCES ONE FILE WITH FORTRAN C REFERENCE NUMBER LU1 (IF LU1.NE.0) AND ANOTHER WITH REFERENCE C NUMBER LU2 (IF LU2.NE.0). THIS POSSIBILITY CAN BE USED ONLY IN C OPERATING SYSTEMS ALLOWING TO SPECIFY THE NUMBER OF A FILE c OUTSIDE THE SOURCE PROGRAM (E.G., IBM OPERATING SYSTEM). C TO TERMINATE THE COMPUTATIONS, PUT ICONT=0 IN THE CARD NO.9. C IF ICONT=-1, THE COMPUTATIONS CONTINUE, BUT FOR A NEW MODEL. C THE CARDS ARE THEN READ IN STARTING FROM THE CARD NO.1. C C ****************************************************************** C C C MOST IMPORTANT OUTPUT INFORMATION. C ********************************** C C C THE PARAMETER IND. C ****************** C THE PARAMETER IND SPECIFIES THE REASON FOR THE TERMINATION OF C THE COMPUTED RAY OR TERMINATION OF COMPUTATIONS AT ALL. PARAMETER C IND IS STORED IN COMMON/RAY/, SEE BELOW. C C IND=1 TERMINATION OF THE RAY AT THE LEFT BORDER OF THE SUBMODEL. C IND=2 TERMINATION OF THE RAY AT THE RIGHT BORDER OF THE SUBMODEL. C IND=3 TERMINATION OF THE RAY AT THE UPPER BORDER OF THE MODEL C (I.E., ON THE FREE SURFACE). C IND=4 TERMINATION OF THE RAY AT THE BOTTOM BORDER OF THE MODEL. C IND=5 IN THE RUNGE-KUTTA METHOD, THE BASIC TIME STEP DT IS C HALVED MORE THAN TEN TIMES, WITHOUT REACHING THE REQUIRED C ACCURACY AC, SEE CARD NO. 12. THERE IS PERHAPS SOMETHING C WRONG WITH YOUR VELOCITY DISTRIBUTION IN THE MODEL. C IND=6 IN CARD NO. 12, DT=0 (NOT ALLOWED). THE COMPUTATION C OF THE RAY DIAGRAM DOES NOT START. C IND=7 IN CARD NO. 12, DT.LT.0 (NOT ALLOWED). THE COMPUTATION C OF THE RAY DIAGRAM DOES NOT START. C IND=8 THE RAY STRIKES A DIFFERENT INTERFACE THAN INDICATED IN C THE NUMERICAL CODE OF THE WAVE. C IND=9 OVERCRITICAL INCIDENCE AT AN INTERFACE WHERE THE NUMERICAL C CODE REQUIRES A TRANSMISSION. C IND=12 TRAVEL TIME ALONG THE RAY IS GREATER THAN TMAX. IN THE C PROGRAM SEIS88, TMAX=10000 AUTOMATICALLY. C IND=13 THE NUMBER OF POINTS ALONG THE RAY IS LARGER THAN 400. C IND=14 THE POSITION OF THE SOURCE DOES NOT CORRESPOND TO THE C NUMERICAL CODE OF THE WAVE, I.E. ABS(CODE(1)) IS DIFFE- C RENT FROM THE NUMBER OF THE LAYER WITHIN WHICH THE SOURCE C IS SITUATED. THE SAME INDICATION APPEARS IF THE SOURCE IS C SITUATED IN THE LAYER OF ZERO THICKNESS. (FOR DETAILS, C SEE CARD NO.11). C IND=15 IN THE CASE WHERE THE COMPUTATION OF THE AMPLITUDES C IS REQUIRED (MDIM.NE.0, SEE CARD NO. 9 FOR MDIM), AND THE C NUMERICAL CODE OF THE WAVE REQUIRES THE COMPUTATION OF A C WAVE REFLECTED FROM THE BOTTOM BOUNDARY OF THE MODEL, THE C COMPUTATION OF THE RAY IS TERMINATED, AS THE VELOCITIES C AND THE DENSITY BELOW THIS BOUNDARY ARE NOT KNOWN. C IND=16 THE NUMERICAL CODE REQUIRES A COMPUTATION OF A WAVE C REFLECTED FROM A FICTITIOUS INTERFACE (III=-2). THE C AMPLITUDE OF SUCH A WAVE EQUALS ZERO, THEREFORE THE C COMPUTATION OF THE RAY IS TERMINATED. THE SAME INDICATION C APPEARS IF THE TERMINATION POINT OF THE RAY LIES ON C A FICTITIOUS PART OF AN INTERFACE. C IND=17 THE NUMERICAL CODE REQUIRES A REFLECTION FROM THE SECOND C INTERFACE OF TWO COINCIDING INTERFACES. THE COMPUTATION C OF THE RAY IS TERMINATED. THE SAME APPLIES TO OTHER IN- C TERFACES, WHEN THE NUMBER OF COINCIDING INTERFACES IS C LARGER THAN TWO. C IND=18 THE NUMERICAL CODE PRESCRIBES CONVERSION OF WAVE TRANS- C MITTED ON AN INTERFACE WHICH COINCIDES WITH ANOTHER INTER- C FACE. THE COMPUTATION OF THE RAY IS TERMINATED SINCE C KINEMATICALLY AND DYNAMICALLY IDENTICAL RAY CAN BE OBTAI- C NED WITH NUMERICAL CODE WHICH DOES NOT INCLUDE THE CORRES- C PONDING CONVERSION. C IND=19 THE ROUTINE MODEL USED FOR THE COMPUTATIONS DOES NOT COR- C RESPOND TO THAT SPECIFIED IN THE INPUT DATA, SEE PARAMETER C MM IN INPUT DATA CARD NO.1. THE COMPUTATION DOES NOT START. C IND=20 TWO NEIGHBOURHOOD INTERFACES IN THE MODEL INTERSECT EACH C OTHER. THE COMPUTATION TERMINATES. C IND=21 A RECTANGULAR NETWORK OF VELOCITY VALUES FOR ONE LAYER C DOES NOT COVER THE WHOLE LAYER. THE COMPUTATION TERMINATES. C IND=22 TERMINATION OF THE RAY AT THE VERTICAL PROFILE. C IND=50 THE SOURCE IS SITUATED OUTSIDE THE MODEL, THE COMPUTATION C OF THE RAY DIAGRAM DOES NOT START. C IND.GT.100 IF THE NUMERICAL CODE OF THE WAVE REQUIRES THE TER- C MINATION OF THE COMPUTATION AT SOME INNER INTERFACE AND C THE COMPUTATION OF THE RAY IS SUCCESSFULLY FINISHED, THEN C (IND-100) GIVES THE NUMBER OF THE INTERFACE AT WHICH THE C COMPUTATION OF THE RAY IS TERMINATED. C C C THE HISTORY OF THE RAY. C *********************** C C C ALL THE IMPORTANT INFORMATION ABOUT ANY COMPUTED RAY IS STORED IN C COMMON/RAY/AY(12,400),DS(11,50) KINT(50),MREG,N,IREF,IND,IND1. C THIS INFORMATION IS LOST AS SOON AS A NEW RAY BEGINS TO BE COM- C PUTED. A MAXIMUM NUMBER OF 400 POINTS ALONG THE RAY IS ALLOWED. C C THE MEANING OF THE ABOVE QUANTITIES IS THE FOLLOWING: C C N... NUMBER OF POINTS ALONG THE RAY. C IREF... NUMBER OF REFLECTIONS/TRANSMISSIONS, INCLUDING C REFLECTIONS/TRANSMISSIONS AT FICTITIOUS INTERFA- C CES. C IND... EXPLAINED ABOVE. C IND1... IABS(IND1)... THE NUMBER OF THE LAST INTERFACE C HIT BY THE RAY BEFORE THE COMPUTATION OF THE C RAY WAS TERMINATED (OR THE INTERFACE AT WHICH C THE COMPUTATION WAS TERMINATED). C IND1=0... THE COMPUTATION OF THE RAY TERMINATED C ON ONE OF THE VERTICAL BOUNDARIES OF THE MODEL. C THE RAY DID NOT CROSS ANY INTERFACE ON ITS WAY C FROM THE SOURCE TO THE TERMINATION POINT. C IND1.LT.0... THERE IS AT LEAST ONE COMPOUND ELE- C MENT ALONG THE RAY. C MREG... CONTROLS THE COMPUTATION OF AMPLITUDES AT THE TER- C MINATION POINT OF THE RAY,SEE INPUT CARD NO. 9. C C THE QUANTITIES AY(J,I) GIVE SOME INFORMATION ABOUT THE I-TH COM- C PUTED POINT OF THE RAY. I=1 CORRESPONDS TO THE SOURCE, I=N TO C THE TERMINATION POINT. C AY(1,I)...THE TIME CORRESPONDING TO THE POINT. C AY(2,I)...X-COORDINATE OF THE POINT. C AY(3,I)...Z-COORDINATE OF THE POINT. C AY(4,I)...X COMPONENT OF THE SLOWNESS VECTOR AT THE POINT. C AY(5,I)...Z COMPONENT OF THE SLOWNESS VECTOR AT THE POINT. C AY(6,I)...VELOCITY V AT THE POINT. C AY(7,I)...DV/DX C AY(8,I)...DV/DZ. C AY(9,1)...SECOND DERIVATIVE OF V WITH RESPECT TO X. C AY(10,I)..SECOND DERIVATIVE OF V WITH RESPECT TO X AND Z. C AY(11,I)..SECOND DERIVATIVE OF V WITH RESPECT TO Z. C AY(12,I)..QUALITY FACTOR AT THE POINT. C C C THE QUANTITIES DS(J,K) AND KINT(K) GIVE SOME INFORMATION ABOUT THE C K-TH POINT OF REFLECTION/TRANSMISSION. THE CASE K=1 CORRESPONDS C TO THE INTERFACE WHICH WAS HIT BY THE RAY FIRST. C C DS(1,K)...THE RADIUS OF THE CURVATURE OF THE INTERFACE. C DS(2,K)...X COMPONENT OF THE UNIT VECTOR TANGENT TO THE RAY C AT THE POINT OF INCIDENCE, IN LOCAL COORDINATE C SYSTEM CONNECTED WITH THE POINT OF INCIDENCE. C Z-AXIS OF THIS SYSTEM IS PERPENDICULAR TO THE C INTERFACE AND POINTS INTO THE MEDIUM IN WHICH C INCIDENT WAVE PROPAGATES. Y-AXIS IS PERPENDICULAR C TO THE PLANE OF THE MODEL, ITS POSITIVE DIRECTION C POINTS BEHIND THE (X,Z) PLANE. X-AXIS IS TANGENT C TO THE INTERFACE, ITS POSITIVE DIRECTION IS C CHOSEN SO THAT THE LOCAL COORDINATE SYSTEM IS C RIGHT-HANDED. C IN CASE OF TERMINATION POINT OF THE RAY ON ONE C OF THE VERTICAL BOUNDARIES OF THE MODEL, C DS(2,K)=AY(4,N). C DS(3,K)...Z COMPONENT OF THE UNIT VECTOR TANGENT TO THE RAY C AT THE POINT OF INCIDENCE, IN LOCAL COORDINATE C SYSTEM CONNECTED WITH THE POINT OF INCIDENCE C (FOR DETAILS SEE ABOVE). C IN CASE OF TERMINATION POINT OF THE RAY ON ONE C OF THE VERTICAL BOUNDARIES OF THE MODEL, C DS(3,K)=AY(5,N). C DS(4,K)...VP AT THE POINT OF INCIDENCE. C DS(5,K)...VS AT THE POINT OF INCIDENCE. C DS(6,K)...VP AT OTHER SIDE OF INTERFACE. C DS(7,K)...VS AT OTHER SIDE OF INTERFACE. C DS(8,K)...DENSITY AT THE POINT OF INCIDENCE. C DS(9,K)...DENSITY AT THE OTHER SIDE OF INTERFACE. C DS(10,K)..THE SAME AS FOR DS(2,K) BUT FOR GENERATED WAVE. C DS(11,K)..THE SAME AS FOR DS(3,K) BUT FOR GENERATED WAVE. C C KINT(K)...KINT(K).GT.0... GIVES THE SEQUENTIAL NUMBER OF C THE POINT OF THE RAY WHICH CORRESPONDS TO THE C K-TH POINT OF REFLECTION/TRANSMISSION. C KINT(K)=0... THE K-TH POINT OF REFLECTION/ C TRANSMISSION HAS A FORMAL MEANING, AS IT CORRES- C PONDS TO A COMPOUND ELEMENT OF THE RAY, WITH C BOTH END POINTS OF THE ELEMENT AT THE SAME INTER- C FACE. C KINT(K)=-1... THE K-TH POINT OF REFLECTION/ C TRANSMISSION CORRESPONDS TO A POINT ON THE INTER- C FACE WHICH COINCIDES WITH THE INTERFACE SITUATED C ABOVE IT. C C C DESCRIPTION OF INTERFACES. C ************************** C C C THE INFORMATION ABOUT THE INTERFACES IS STORED IN C COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20), C BRD(2),III(30,20),NPNT(20),NINT. C C NINT... SPECIFIED IN INPUT DATA CARD NO.2. C NPNT(I)...SPECIFIED IN INPUT DATA CARD NO.2. C BRD(1)... X-COORDINATE OF THE LEFT VERTICAL BOUNDARY C OF THE SUBMODEL. SPECIFIED IN INPUT CARD NO.11. C BRD(2)... X-COORDINATE OF THE RIGHT VERTICAL BOUNDARY C OF THE SUBMODEL. SPECIFIED IN INPUT CARD NO.11. C C BETWEEN TWO SUCCESSIVE POINTS OF THE INTERFACE, WITH THE X-COOR- C DINATE OF THE LEFT-HAND POINT DENOTED BY X1, THE CUBIC PARABOLA C APPROXIMATING THE INTERFACE HAS A FORM: C Z = A+B*(X-X1)+C*(X-X1)**2+D*(X-X1)**3 C C X1(K,I)...THE X-COORDINATE OF THE K-TH POINT OF THE I-TH C INTERFACE. C A1(K,I)...THE Z-COORDINATE OF THE K-TH POINT OF THE I-TH C INTERFACE (CORRESPONDS TO A IN THE ABOVE FORMULA). C B1(K,I),C1(K,I),D1(K,I)... THE COEFFICIENTS OF THE CUBIC C PARABOLA IN THE K-TH ELEMENT OF THE I-TH INTER- C FACE. THEY CORRESPOND TO B,C,D IN THE ABOVE FOR- C MULA. THE K-TH ELEMENT IS BETWEEN THE K-TH AND C (K+1)TH POINT. C III(K,I)..HAS THE SAME MEANING AS III READ FROM THE INPUT C CARD NO.3. C C C DESCRIPTION OF VELOCITY, DENSITY AND QUALITY FACTOR DISTRIBUTIONS. C ****************************************************************** C C C INFORMATION ABOUT THE VELOCITY, DENSITY AND QUALITY FACTOR C DISTRIBUTION IS STORED IN THE FOLLOWING COMMONS: C COMMON/MEDIM/V(1000),SX(350),SY(350),NX(20),NY(20),NVS(19), C PTOS(19) FOR MM=0,2, SEE INPUT CARD NO.1, C COMMON/MEDIM/V1(19),V2(19),NVS(19),PTOS(19) FOR MM=1, C COMMON/VCOEF/A02(1000),A20(1000),A22(1000) FOR MM=0,2, C COMMON/DEN/RHO1(19),RHO2(19),NRHO. C COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19), C NQP(19),NQS(19),NABS, C COMMON/SOUR/ROS,VPS,VSS. C C V... VELOCITY VALUES SPECIFIED IN INPUT DATA CARD NO.4D C FOR MM=0,2. C SX,SY... COORDINATES OF VERTICAL AND HORIZONTAL LINES IN THE C VELOCITY NETWORK, SPECIFIED IN INPUT DATA CARDS C NO.4B AND 4C FOR MM=0,2. C NX,NY... NUMBER OF VERTICAL AND HORIZONTAL LINES IN THE C VELOCITY NETWORK, SPECIFIED IN INPUT DATA CARDS C NO.4A (SEE MX,MY) FOR MM=0,2. C NVS... SEE A DETAILED DESCRIPTION IN DATA CARD NO.5. C PTOS... RATIOS OF P AND S VELOCITIES IN INDIVIDUAL LAYERS, C SEE DETAILS IN CARD NO. 7. C V1,V2... VELOCITY VALUES SPECIFIED IN INPUT DATA CARDS NO. C 4A AND 4B FOR MM=1. C A02,A20,A22.. BASIC COEFFICIENTS OF THE BICUBIC SPLINE OR BI- C LINEAR INTERPOLATION AT GRID POINTS SPECIFIED IN THE C SAME WAY AS FOR V (SEE DESCRIPTION UNDER THE INPUT C DATA CARD NO.4D FOR MM=0,2). A02 CORRESPONDS TO THE C FIRST (SECOND) DERIVATIVE OF VELOCITY WITH RESPECT C TO Z, A20 WITH RESPECT TO X AND A22 TO THE SECOND C (FOURTH) DERIVATIVE WITH RESPECT TO X AND Z IN CASE C OF BILINEAR (BICUBIC SPLINE) INTERPOLATION. C RHO1,RHO2... QUANTITIES THAT SPECIFY THE RELATION BETWEEN C THE DENSITIES AND P WAVE VELOCITIES IN INDIVI- C DUAL LAYERS, DENSITY=RHO1(I)+RHO2(I)*VP IN THE C I-TH LAYER, SEE CARD NO.6A. FOR NRO=0 (SEE CARD C NO.5 FOR NRO) RHO1(I)=1.7, RHO2(I)=0.2 IN ALL C LAYERS. FOR PTOS(1).GE.100 (SEE CARD NO. 7 FOR C PTOS(1)), RHO1(1)=1,RHO2(1)=0 (A LIQUID LAYER). C NRHO... FOR NRHO=0, THE UPPERMOST LAYER IS SOLID, FOR C NRHO=1, THE UPPERMOST LAYER IS LIQUID. IN THE C LATTER CASE, THE S VELOCITY=0 AND THE DENSITY=1 C IN THE LAYER. ALSO MREG=1 IN THIS CASE. C QP1,QP2,QP3,QS1,QS2,QS3... QUANTITIES WHICH SPECIFY THE RE- C LATION BETWEEN THE QUALITY FACTOR AND THE CORRES- C PONDING VELOCITY, SEE INPUT DATA CARD NO.6B. C NQP,NQS...SWITCHES WHICH DETERMINE THE TYPE OF THE RELATION C BETWEEN THE QUALITY FACTOR AND THE CORRESPONDING C VELOCITY, SEE INPUT DATA CARD NO.6B. C ROS,VPS,VSS... VALUES OF DENSITY, P-WAVE AND S-WAVE VELOCITY C AT THE SOURCE. C NOTE: FOR MM=0,2 THE INFORMATION ABOUT VELOCITY DISTRIBUTION C IN ONE LAYER IS STORED IN THE ABOVE VARIABLES IN THE SAME C WAY AS DESCRIBED IN INPUT DATA CARDS NO.4. THE DATA FOR INDI- C VIDUAL LAYERS ARE STORED SUCCESSIVELY FROM THE TOP TO THE C BOTTOM, ONE AFTER ANOTHER. C C DESCRIPTION OF CERTAIN AUXILIARY QUANTITIES C ******************************************* C C C AUXILIARY QUANTITIES WHICH ARE USED IN VARIOUS PARTS OF THE PRO- C GRAM ARE STORED IN COMMON /AUXI/ AND COMMON /AUXX/. C C COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT, C MPRINT,NTR C C INTR... NUMBER OF THE INTERFACE JUST STRUCK BY THE RAY. C INT1... NUMBER OF THE INTERFACE WHICH WAS HIT BEFORE INTR. C IOUT... INDEX USED WHEN THE POINT OF INCIDENCE OF A RAY C AT AN INTERFACE IS DETERMINED IN AN ITERATIVE WAY C (NOT USED IN THIS VERSION). C KRE... INDEX SPECIFYING WHAT TYPE OF WAVE IS CONSIDERED. C KRE=0... REFRACTED WAVE IS CONSIDERED. C KRE=KREF... (SEE CARD NO.14 FOR KREF), MULTIPLY C REFLECTED WAVE IS CONSIDERED. C IREFR... INDEX SPECIFYING IF THERE IS AT LEAST ONE COMPOUND C ELEMENT OF THE RAY (IREFR=1) OR NOT (IREFR=0). C LAY... NUMBER OF LAYER JUST PASSED THROUGH THE RAY. C ITYPE... INDEX SPECIFYING THE TYPE OF THE WAVE ALONG THE C INVESTIGATED ELEMENT OF THE RAY. ITYPE=1 FOR P C WAVE, ITYPE=-1 FOR S WAVE. C NDER... INDEX SPECIFYING WHETHER THE SECOND DERIVATIVES OF C VELOCITY SHOULD BE EVALUATED (NDER=1 FOR MDIM.GT.0) C OR NOT (NDER=0 FOR MDIM=0). C IPRINT...BLANK PARAMETER IN THIS VERSION. C MPRINT...CONTROLS THE PRINTING OF THE DESCRIPTION OF THE C MODEL ON THE LINE PRINTER, SEE INPUT DATA CARD NO.1. C NTR... SERVES FOR THE PURPOSE OF A MORE DETAILED DETER- C MINATION OF THE REASON OF THE TERMINATION OF THE RAY. C IT SHOWS THE PLACE IN THE ROUTINE OUTP WHERE THE COM- C PUTATION OF THE RAY WAS TERMINATED (SEE THE LISTING C OF THE ROUTINE OUTP). C C C COMMON/AUXX/MMX(20),MMY(20),MMXY(20),MAUX C C USED ONLY IN THE ROUTINE MODEL AND RELATED ROUTINES IN CASE C THAT MM=0 OR 2. C C MMX(I)...SPECIFIES THE NUMBER OF THE FIRST COMPONENT OF THE C VECTOR SX (SEE /MEDIM/) CORRESPONDING TO THE I-TH C LAYER. C MMY(I)...SPECIFIES THE NUMBER OF THE FIRST COMPONENT OF THE C VECTOR SY (SEE /MEDIM/) CORRESPONDING TO THE I-TH C LAYER. C MMXY(I)...SPECIFIES THE NUMBER OF THE FIRST COMPONENTS OF C THE VECTORS V (SEE /MEDIM/) AND AO2,A20,A22 (SEE C /VCOEF/) CORRESPONDING TO I-TH LAYER. C MAUX... INDEX SPECIFYING IF IT IS NECESSARY TO EVALUATE C THE COEFFICIENTS OF BICUBIC SPLINE OR BILINEAR C INTERPOLATION,OR WHETHER IT IS POSSIBLE TO USE THE C COEFFICIENT DETERMINED IN THE PREVIOUS STEP. C C C ****************************************************************** C C OUTPUT TABLES. C ************** C C ALL THE INPUT DATA ARE REPRODUCED ON THE LINE PRINTER. THE C OTHER OUTPUT IS CONTROLLED BY THE PARAMETERS MPRINT AND C MOUT, SEE INPUT CARDS NO.1 AND 9. C C C THE INDEX MPRINT: C **************** C THE INDEX MPRINT CONTROLS THE PRINTING OF THE DESCRIPTION OF THE C MODEL ON THE LINE PRINTER. C FOR MPRINT=0: C ************* C INPUT DATA FOR THE MODEL ARE REPRODUCED. C FOR MPRINT=1: C ************ C PRINTER PLOT OF THE VELOCITY DISTRIBUTION WITHIN THE PART OF THE C MODEL BOUNDED BY HORIZONTAL BOUNDARIES BMIN AND BMAX AND BY VER- C TICAL BOUNDARIES BLEFT AND BRIGHT, SEE DETAILS IN THE DESCRIPTION C OF THE INPUT CARD NO.8. THE PRINTERPLOT SHOWS APPROXIMATELY THE C POSITION OF INTERFACES AND/OR THE POSITION OF VELOCITY ISOLINES C IN THE MODEL. C FOR MPRINT=2: C ************ C THE SAME PRINTOUT AS IN THE CASE OF MPRINT=1, BUT WITH A DETAILED C DESCRIPTION OF INDIVIDUAL INTERFACES. FOR EACH INTERFACE, THE C FOLLOWING QUANTITIES ARE PRINTED: C A) COEFFICIENTS OF CUBIC PARABOLAS CORRESPONDING TO INDIVIDUAL C ELEMENTS OF THE INTERFACE (FROM THE LEFT TO THE RIGHT VERTICAL C BOUNDARY OF THE WHOLE MODEL). C B) COORDINATES OF 26 POINTS ALONG EACH INTERFACE, WITH REGULAR C STEP IN THE X-COORDINATE, FROM THE LEFT BOUNDARY 'BLEFT' TO THE C RIGHT BOUNDARY 'BRIGHT', SEE INPUT DATA CARD NO.8. C C C THE INDEX MOUT: C ************** C C THE INDEX MOUT CONTROLS THE PRINTING OF A DESCRIPTION OF THE WHOLE C PROCESS OF ITERATIONS (TWO-POINT RAY TRACING) AND OF ITS RESULTS C ON A LINE PRINTER. C FOR MOUT=0: C ********** C NO PRINT OF RESULTS. ONLY THE LISTING OF ELEMENTARY WAVES, FOR C WHICH THE COMPUTATIONS ARE PERFORMED, IS GIVEN. THE "EXTERNAL C CODE" (FULL CODE) CORRESPONDS TO THE NUMERICAL CODE OF THE WAVE C DESCRIBED IN INPUT DATA CARD NO.14. THE "INTERNAL CODE" IS THE C SUCCESSIVE NUMBER OF THE ELEMENTARY WAVES GENERATED BY THE ROUTINE C GENER. FOR MORE DETAILS, SEE DESCRIPTION OF THE ROUTINE GENER. C FOR MOUT=1 OR MOUT=2: C ******************** C IN ADDITION TO THE NUMERICAL CODES OF INDIVIDUAL ELEMENTARY WAVES, C CERTAIN RESULTS OF ITERATIONS ARE PRINTED AT THE LINE PRINTER. C THE PRINTING IS DIFFERENT FOR MOUT=1 (ONLY BASIC PRINT) AND C MOUT=2 (HEAVY PRINTING). C THE FOLLOWING QUANTITIES ARE PRINTED IN VARIOUS SITUATIONS: C IND, IND1... SEE EXPLANATION IN THE TABLE "THE PARAMETER C IND" AND IN "THE HISTORY OF RAY". C ITER... THE SUCCESSIVE NUMBER OF THE ITERATION. IN THE C PRINT OF THE FINAL RESULT OF ITERATIONS, IT C GIVES THE FULL NUMBER OF ITERATIONS. C II... THE SUCCESSIVE NUMBER OF THE RECEIVER POSITION C CORRESPONDING TO THE TERMINATION POINT OF THE C RAY. THE RECEIVERS ARE NUMBERED FROM THE LEFT C TO THE RIGHT. C INDEX... THE NUMBER OF SUCCESSFULLY ITERATED RAYS FOR A C GIVEN ELEMENTARY WAVE. C IRES... IRES=1... BOUNDARY RAY TERTMINATING ON THE RE- C CEIVER PROFILE WAS FOUND. C IRES=0... BOUNDARY RAY TERMINATING ON THE RE- C CEIVER PROFILE WAS NOT FOUND. C KMAH... KMAH INDEX, SEE DETAILS IN THE DESCRIPTION OF C THE ROUTINE JPAR. C X,Z... COORDINATES OF THE TERMINATION POINT OF THE RAY. C T... TRAVEL TIME AT THE TERMINATION POINT OF THE RAY. C TAST... GLOBAL ABSORPTION FACTOR AT THE TERMINATION POINT C OF THE RAY. C AA... PARAMETER OF THE RAY (INITIAL ANGLE) IN THE BASIC C SYSTEM OF INITIAL ANGLES IN THE RAY SHOOTING, SEE C DATA CARD NO.12. C PNEW... THE PARAMETER OF THE RAY (INITIAL ANGLE) IN ITE- C RATIONS. C DD... COORDINATE OF THE RECEIVER (X-COORDINATE FOR RE- C CEIVERS SITUATED ALONG THE EARTH'S SURFACE OR AN C INTERFACE, Z-COORDINATE FOR RECEIVERS SITUATED C ALONG A VERTICAL PROFILE). C AX,AY,AZ... AMPLITUDES OF RADIAL (HORIZONTAL ALONG THE PRO- C FILE, POSITIVE TO THE RIGHT), TRANSVERSE (HORI- c ZONTAL, PERPENDICULAR TO THE PROFILE) AND VERTI- C CAL (POSITIVE UPWARDS) COMPONENTS OF DISPLACEMENT C VECTOR AT THE TERMINATION POINT OF RAY (AX,AY,AZ C SPECIFIED IN A RIGHT-HANDED COORDINATE SYSTEM). C PHX,PHY,PHZ... PHASE SHIFTS CORRESPONDING TO THE ABOVE COMPO- C NENTS AT THE TERMINATION POINT OF THE RAY. C SPR... THE GEOMETRICAL SPREADING AT THE TERMINATION C POINT OF THE RAY. C C ALL TABLES ARE PRINTED WITHOUT HEADINGS, BUT CERTAIN LINES ARE C SUPPLEMENTED BY THE ADDITIONAL TEXT: "ITERATIONS", "SUCCESSFUL C ITERATION", "BOUNDARY RAY", "CRITICAL RAY". C GENERALLY, FOR EACH RAY ONE LINE IS PRINTED. THE EXCEPTION IS C THE SUCCESSFULLY ITERATED RAY, WHICH MAY BE SUPPLEMENTED BY SOME C ADDITIONAL INFORMATION ON AMPLITUDES (WHEN MDIM.GT.0,SEE LATER). C IN THE FOLLOWING, THE SYMBOL WW MEANS EITHER X (WHEN RECEIVERS C ARE SITUATED ALONG THE EARTH'S SURFACE OR ALONG AN INTERFACE), C OR Z (WHEN RECEIVERS ARE SITUATED ALONG A VERTICAL PROFILE). C C IN THE BASIC SYSTEM OF SHOOTING, THE FOLLOWING QUANTITIES C ARE PRINTED WHEN MOUT=2: C IND,IND1,W,T,AA. C IN ITERATIONS, FOR MOUT=2: C "ITERATIONS", IND,IND1,ITER,DD,W,T,PNEW. C IN SUCCESSFUL ITERATIONS, WHEN MOUT=2: C "SUCCESSFUL ITERATION",IND,IND1,ITER,II,INDEX,DD,W,T,PNEW, C AX,PHX,AY,PHY,AZ,PHZ,SPR,TAST,KMAH. C IN SUCCESSFUL ITERATIONS, WHEN MOUT=1 AND MDIM=0: C DD,X,Z,T,PNEW,IND,IND1,ITER,II. C IN SUCCESSFUL ITERATIONS, WHEN MOUT=1 AND MDIM.NE.0: C DD,X,Z,T,TAST,PNEW,SPR,AX,PHX,AY,PHY,AZ,PHZ,IND,IND1,ITER, C II,KMAH. C IN CASE OF BOUNDARY RAYS, ALL ITERATIONS ARE PRINTED WHEN MOUT=2: C "BOUNDARY RAY",IND,IND1,ITER,W,T,PNEW. C IN CASE OF BOUNDARY RAYS, ONLY THE FINAL ITERATION IS PRINTED WHEN C MOUT=1: C X,Z,T,PNEW,IND,IND1,IRES. C IN CASE OF CRITICAL RAYS, ALL ITERATIONS ARE PRINTED WHEN MOUT=2: C "CRITICAL RAY", IND,IND1,ITER,W,T,PNEW. C IN CASE OF CRITICAL RAYS, ONLY THE FINAL ITERATION IS PRINTED C WHEN MOUT=1: C X,Z,T,PNEW,IND,IND1,IRES. C C C C OUTPUT TO THE FILE LU1. C *********************** C C C IN THE CASE THAT LU1.NE.0, CERTAIN INFORMATION ABOUT THE MODEL, C AND THE MOST IMPORTANT RESULTS OF COMPUTATION ARE STORED IN THE C FILE LU1 FOR POSSIBLE FURTHER PROCESSING IN OTHER PROGRAMS. C THE PROGRAM RAYPLOT, INCLUDED IN THIS PACKAGE, MAY BE USED FOR C PLOTTING OF DATA STORED IN THE FILE LU1. THE DATA ARE STORED IN C LU1 IN UNFORMATTED FORM. C C WHEN LU1=0, NO DATA ARE STORED. C C THE DATA CONSIST OF THE FOLLOWING PARTS: C C A) INFORMATION ABOUT POSSIBLE TERMINATION OF COMPUTATION, ABOUT C THE SOURCE AND ABOUT INTERFACES: C 1) ICONT. C ICONT=0... INDICATION OF THE END OF DATA STORED IN THE FILE C LU1. THIS IS THE LAST QUANTITY STORED IN THE FILE LU1. C ICONT=1... INDICATION THAT A NEW DATA SET (2-8) FOR A NEW RAY C DIAGRAM MAY FOLLOW. C 2) NINT, (NPNT(I),I=1,NINT). C FOR THE MEANING OF THESE SYMBOLS SEE INPUT DATA CARD NO.2. C 3) FOR I=1,2,...,NINT: C (A1(J,I),B1(J,I),C1(J,I),D1(J,I),X1(J,I),III(J,I),J=1,NC), C WHERE NC=NPNT(I)-1. C FOR THE MEANING OF THESE SYMBOLS SEE INPUT DATA CARD NO.3 C AND THE DESCRIPTION OF INTERFACES. C C B) INFORMATION ABOUT THE SOURCE POSITION AND PARAMETERS OF THE C MEDIUM AT THE SOURCE: C 4) X0,Z0,ROS,VPS,VSS. C COORDINATES OF THE SOURCE, THE DENSITY AND P- AND S-WAVE VE- C LOCITY AT THE SOURCE. C C C) INFORMATION ABOUT RAYS: C FOR EACH SUCCESSFULLY ITERATED RAY WITH THE TERMINATION C POINT AT SOME RECEIVER AT THE EARTH'S SURFACE,THE QUAN- C TITIES SPECIFIED SUB 5 AND 6 ARE STORED. THE DATA SPECI- C FIED IN 5 AND 6 ARE STORED SUCCESSIVELY FOR ALL RAYS WITHIN C ONE RAY DIAGRAM, AN ARBITRARY NUMBER OF TIMES. AFTER THE C LAST RAY, THE PARAMETER N IS SET TO ZERO IN THE MAIN C PROGRAM. THEN THE QUANTITIES DESCRIBED SUB D ARE STORED. C 5) N,IND C IF(N.NE.0)... NUMBER OF POINTS ALONG THE RAY. C IF N=0... INDICATION OF THE END OF THE RAY DIAGRAM. C IND... THE INDEX SPECIFYING THE REASON OF THE TERMINATION C OF THE RAY. SEE DETAILED EXPLANATION ABOVE. WHEN THE RE- C CEIVERS ARE SITUATED ALONG THE EARTH'S SURFACE, THE RAYS C WITH IND=3 ARE STORED IN LU1. WHEN THE RECEIVERS ARE SI- C TUATED ALONG THE (IND-100)TH INTERFACE, THE RAYS WITH C IND.GT.100 ARE STORED IN LU1. WHEN THE RECEIVERS ARE SI- C TUATED ALONG A VERTICAL PROFILE, THE RAYS WITH IND=22 C ARE STORED IN LU1. C 6) (AY(2,I),AY(3,I),I=1,N). C FOR THE MEANING OF THESE SYMBOLS SEE THE DESCRIPTION OF C COMMON/RAY/IN THE HISTORY OF THE RAY. C C D) INFORMATION ABOUT TRAVEL-TIMES AND AMPLITUDES C 7) NS. C NUMBER OF RAYS WHICH TERMINATE AT THE SPECIFIED SYSTEM OF C RECEIVERS. C 8) (INDI(I),XCOOR(I),TIME(I),TAS(I),ANG(I),AMPX(I),AMPY(I), C AMPZ(I),PHAX(I),PHAY(I),PHAZ(I),I=1,NS). C THE MEANING OF THE ABOVE SYMBOLS IS AS FOLLOWS: C A SPECIAL PARAMETER IND WHICH SPECIFIES THE REASON OF THE TER- C MINATION OF THE COMPUTATION OF THE RAY, THE X-COORDINATE OF C THE TERMINATION POINT OF THE RAY, THE CORRESPONDING TRAVEL C TIME, GLOBAL ABSORPTION FACTOR T*, PARAMETER OF THE RAY (INI- C TIAL ANGLE), AMPLITUDES OF RADIAL, TRANSVERSE AND VERTICAL COM- C PONENTS AND CORRESPONDING PHASE SHIFTS. ALL THESE DATA FOLLOW C SUCCESSIVELY FOR ALL THE RAYS WHICH REACHED THE RECEIVER PRO- C FILE (I=1,2,....,NS). C C C C OUTPUT TO THE FILE LU2. C *********************** C C IN THE FILE LU2, THE DATA FOR COMPUTATION OF RAY SYNTHETIC SEIS- C MOGRAMS FOR INDIVIDUAL RECEIVER POSITIONS ARE STORED TOGETHER C WITH SOME OTHER RELEVANT INFORMATION. ALL THE DATA CORRESPOND C TO ONE SYSTEM OF CARDS NO.9-13. THE DATA MAY BE STORED IN A FOR- C MATTED FORM (LTAPE=0, SEE INPUT DATA CARD NO.9 FOR LTAPE) TO C ALLOW AN INSPECTION OF THE RESULTS OF COMPUTATION. TO INCREASE C THE SPEED OF COMPUTATIONS, IT IS POSSIBLE TO STORE THE DATA IN C A BINARY FORM (LTAPE=1). C TO PERFORM THE COMPUTATION OF SYNTHETIC SEISMOGRAMS AND THEIR C PLOTTING, THE PROGRAMS SYNTPL AND SEISPLOT INCLUDED IN THIS PAC- C KAGE MAY BE USED. C C THE DATA ARE STORED IN THE FOLLOWING SUCCESSION: C C 1) MTEXT FORMAT(17A4) C ARBITRARY ALPHANUMERIC TEXT DESCRIBING THE MODEL. C 2) NDST FORMAT(26I3) C NDST IS A NUMBER OF RECEIVER POSITIONS. C 3) XSOUR,ZSOUR,TSOUR,RSTEP,ROS,VPS,VSS FORMAT(8F10.5) C XSOUR,ZSOUR... X-COORDINATE AND Z-COORDINATE OF THE SOURCE C TSOUR... INITIAL (HYPOCENTRAL) TIME. C RSTEP... (AVERAGE) DIFFERENCE BETWEEN THE X-COORDINATES C OF TWO NEIGHBOURHOOD RECEIVER POSITIONS. C ROS... THE DENSITY AT THE SOURCE. C VPS,VSS... THE P- AND S-WAVE VELOCITY AT THE SOURCE. C 4) DST(I),I=1,NDST FORMAT(8F10.5) C COORDINATES OF RECEIVER POSITIONS (X-COORDINATES FOR RECEI- C VERS SITUATED ALONG THE EARTH'S SURFACE OR AN INTERFACE, Z- C COORDINATES FOR RECEIVERS SITUATED ALONG A VERTICAL PROFILE). C 5) NCD,II,T,AX,AY,AZ,PHX,PHY,PHZ,PNEW,TAST C FORMAT(2I3,F10.5,2E12.6,4F10.6) C NCD... NCODE=IABS(NCD), NCODE IS THE INTERNAL CODE OF THE C WAVE (NCODE=1,2,...), SEE DETAILED DESCRIPTION OF C THE INTERNAL CODE IN THE ROUTINE GENER. C NCD.GT.0... P WAVE. C NCD.LT.0... S WAVE. C II... SUCCESSIVE NUMBER OF THE RECEIVER POSITION. C T... ARRIVAL TIME. C AX... AMPLITUDE OF THE RADIAL COMPONENT. C AY... AMPLITUDE OF THE TRANSVERSE COMPONENT. C AZ... AMPLITUDE OF THE VERTICAL COMPONENT. C PHX... PHASE SHIFT OF THE RADIAL COMPONENT. C PHY... PHASE SHIFT OF THE TRANSVERSE COMPONENT. C PHZ... PHASE SHIFT OF THE VERTICAL COMPONENT. C PNEW... THE PARAMETER OF THE RAY (INITIAL ANGLE). C TAST... THE GLOBAL ABSORPTION FACTOR T*. C C THE DATA 5 ARE REPEATED AT LU2 MANY TIMES, DEPENDING ON THE C NUMBER OF ELEMENTARY WAVES AND ON THE NUMBER OF RECEIVER C POSITIONS. C C C ***************************************************************** C C A SHORT DESCRIPTION OF THE UNIVERSAL 2-D RAY TRACING ROUTINE C RAY2(X0,Z0,T0,FI0,X,Z,T,FI,TMAX,DT,AC) C ************************************** C C THE ROUTINE RAY2 COMPUTES THE RAY SPECIFIED BY INITIAL CONDI- C TIONS X0,Z0,T0,FI0. IT RETURNS THE QUANTITIES X,Z,T,FI AT THE C POINT OF TERMINATION OF THE RAY. THE COMPUTATION IS CONTROLLED C BY THE QUANTITIES TMAX,DT,AC. C THE MEANING OF THESE PARAMETERS IS AS FOLLOWS: C C X0,Z0... COORDINATES OF THE INITIAL POINT OF THE RAY. C T0... INITIAL TIME. C FI0... INITIAL ANGLE. C X,Z... COORDINATES OF THE POINT AT WHICH THE COMPUTATION C OF THE RAY IS TERMINATED. C T... CORRESPONDING TIME. C FI... CORRESPONDING ANGLE OF THE RAY. C DT,AC... THE SAME MEANING AS DESCRIBED IN THE INPUT C DATA CARD NO. 12. C TMAX... THE MAXIMUM PERMITTED TRAVEL TIME ALONG THE RAY. C IN THIS VERSION OF THE PROGRAM, TMAX=10000. C C THE HISTORY OF THE RAY IS STORED IN COMMON/RAY/, SEE ABOVE. THE C QUANTITY IND IN THIS COMMON SPECIFIES THE REASON FOR THE TERMINA- C TION OF THE RAY. C ROUTINES USED: RKGS,FCT,OUT,VELOC,ROOT. C THE INFORMATION ABOUT THE MODEL IS TRANSMITTED TO THE RAY TRACING C ROUTINE BY COMMON/INTRF/ AND COMMON/DEN/. ROUTINE VELOC IS USED C TO DETERMINE THE VELOCITY AND ITS DERIVATIVES AS WELL AS QUALITY C FACTOR AT ANY ARBITRARY POINT IN THE MEDIUM. C THE NUMERICAL CODE OF THE WAVE IS TRANSMITTED TO THE RAY TRACING C ROUTINE BY COMMON/COD/. C C C **************************************************************** C C C A SHORT DESCRIPTION OF THE UNIVERSAL TWO-POINT RAY TRACING C ROUTINE C TRAMP(XSOUR,ZSOUR,TSOUR,DT,AC,TMAX,ASTOP,ITMAX,MOUT,MDIM,MPSOUR, C MSSOUR,NCODE,LU1,LU2,METHOD,ISOUR,LTAPE,ITPR). C ********************************************** C C THE ROUTINE TRAMP IS DESIGNED FOR A TWO-POINT RAY TRACING. IT C COMPUTES THE RAYS WHICH ARRIVE AT A SYSTEM OF RECEIVERS DISTRI- C BUTED REGULARLY OR IRREGULARLY ALONG A PROFILE OF RECEIVERS. THE C RECEIVERS CAN BE SITUATED ON THE EARTH'S SURFACE, ALONG INTERFA- C CES, OR ALONG VERTICAL PROFILES. THE MODIFIED SHOOTING METHOD IS C USED. THE RAYS TERMINATING AT SPECIFIED RECEIVERS ARE SOUGHT BY C ITERATIVE USE OF INITIAL-VALUE RAY TRACING. EACH RAY IS SPECIFIED C BY ITS INITIAL ANGLE AT THE SOURCE. THE BASIC SYSTEM OF THE INI- C TIAL ANGLES OF RAYS AT THE SOURCE IS SPECIFIED IN THE INPUT DATA C CARD NO.12. TO BE SURE THAT NO RECEIVERS ARE OUTSIDE THE RANGE C OF INITIAL ANGLES, THE WHOLE RANGE OF ANGLES (E.G. FROM 3.1416 C TO -3.1416) SHOULD BE USED. TO SAVE THE COMPUTING TIME, THE RANGE C MAY BE SELECTED NARROWER. FROM THE COMPUTED TERMINATION POINTS OF C THESE RAYS, THE RAYS ARRIVING AT SPECIFIED RECEIVERS ARE ITERATI- C VELY DETERMINED. THE TRAVEL TIME CURVE MAY BE ARBITRARY, WITH C LOOPS, CAUSTIC, SHADOWS, ETC. C A SPECIAL CARE IS DEVOTED TO CERTAIN IRREGULAR SITUATIONS IN THE C RAY FIELD. THE ROUTINE IS AUTOMATICALLY LOOKING FOR ALL BOUNDARY C RAYS(I.E. THE RAYS CORRESPONDING TO THE BOUNDARIES OF SHADOW C ZONES AT THE EARTHS SURFACE, TO THE LEFT-HAND SIDE AND RIGHT-HAND C SIDE UPPER CORNERS OF THE MODEL, ETC.). A SPECIAL CARE IS ALSO C DEVOTED TO THE PROBLEM OF SLIGHTLY REFLECTED WAVES (WHICH HAVE C RAYS SIMILAR TO HEAD WAVES). AS THE ROUTINE USES ONLY A SINGLE C PRECISION, ALL THE RAYS OF THIS WAVE, ESPECIALLY THOSE RAYS C WHICH ARE VERY CLOSE TO THE CRITICAL RAY, CANNOT BE SOMETIMES C DETERMINED. C C THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS: C C XSOUR,ZSOUR,TSOUR,DT,AC,ASTOP... SEE INPUT CARD NO. 11 AND 12 C (WITH ASTOP=ASP). C TMAX... MAXIMUM TRAVEL TIME PERMITTED. IN SEIS88, TMAX=10000. C ITMAX,MOUT,MPSOUR,MSSOUR,METHOD,LU1,LU2,LTAPE... SEE INPUT C DATA CARD NO.9. C NCODE... INTERNAL WAVE CODE. SEE DETAILS IN THE DESCRIPTION C OF ROUTINE GENER. C ISOUR... NUMBER OF THE LAYER, IN WHICH THE SOURCE IS SITUATED. C ITPR... FOR BURIED RECEIVERS, THE NUMBER OF THE INTERFACE C ON WHICH THE RECEIVERS ARE SITUATED. C C THE INFORMATION ABOUT THE RECEIVER POSITIONS, ABOUT THE BASIC C SYSTEM OF INITIAL ANGLES FOR RAY SHOOTING, AND ABOUT THE REQUIRED C ACCURACY OF COMPUTATIONS (SEE REPS, INPUT CARD NO.11) IS TRANS- C FERRED TO TRAMP BY COMMON /DIST/. C ROUTINES USED: RAY2, AMPL,JPAR. C THE ROUTINE GENERALLY DOES NOT NEED ANY INFORMATION ON THE VELO- C CITY DISTRIBUTION IN THE MODEL AND ON THE HISTORY OF RAYS. IT C USES ONLY THE INFORMATION ON THE X-COORDINATE OF THE TERMINATION C POINT FOR RECEIVERS DISTRIBUTED ALONG THE EARTH'S SURFACE OR ALONG C AN INTERFACE, OR ON THE Z-COORDINATE OF THE TERMINATION POINT FOR C RECEIVERS DISTRIBUTED ALONG A VERTICAL PROFILE, AND THE INFORMA- C TION ON THE INITIAL ANGLES OF RAYS AND ON THE PARAMETER IND. IN C OTHER WORDS, THE ROUTINE CAN WORK WITH ANY OTHER INITIAL VALUE RAY C TRACING ROUTINE WHICH PROVIDES THESE QUANTITIES. IN FACT, IT HAS C BEEN SUCCESSFULLY APPLIED IN VARIOUS OTHER PROGRAMS, EVEN WHERE C THE RAY PARAMETER HAD A DIFERENT MEANING THAN THE INITIAL ANGLE. C IN THE ROUTINE TRAMP PARTS OF TWO FILES ARE GENERATED. THE FIRST C IS THE FILE LU1 (RAY DIAGRAMS TO SPECIFIED RECEIVERS, TRAVEL TIMES C AND AMPLITUDES), AND THE SECOND IS THE FILE LU2 (THE DATA FOR COM- C PUTATION OF RAY SYNTHETIC SEISMOGRAMS). THE CONTENT OF FILES IS C DESCRIBED ABOVE. LIKEWISE, THE OUTPUT ON THE LINE PRINTER IS DE- C SCRIBED ABOVE. C C ****************************************************************** C C C A SHORT DESCRIPTION OF THE UNIVERSAL ROUTINE FOR 2-D DYNAMIC C RAY TRACING ALONG A KNOWN RAY C JPAR(Q0,P0,Q,P,IPRINT,KMAH). C **************************** C C C WHEN THE RAY IS KNOWN AND ITS HISTORY IS STORED IN COMMON/RAY/, C 2-D DYNAMIC RAY TRACING CAN BE PERFORMED ALONG THE RAY BY C THE PROCEDURE JPAR. C THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS: C C Q0,P0... INITIAL VALUES FOR THE DYNAMIC RAY TRACING C SYSTEM. IN THE PRESENT VERSION Q0=0 AND P0=1./V0, C WHERE V0 IS THE CORRESPONDING WAVE VELOCITY AT THE C SOURCE. C Q,P... THE SOLUTIONS OF THE DYNAMIC RAY TRACING SYSTEM C AT THE TERMINATION POINT OF THE RAY. THE QUAN- C TITY Q REPRESENTS THE SPREADING FUNCTION FOR A C LINE SOURCE. THE SPREADING FUNCTION FOR THE POINT C SOURCE CAN BE EASILY OBTAINED BY WELL-KNOWN METHODS. C IPRINT.. CONTROLS THE PRINTING OF DYNAMIC RAY TRACING ALONG C THE RAY ON A LINE PRINTER. IN SEIS83, NO PRINT OF DY- C NAMIC RAY TRACING RESULTS IS PERFORMED, IPRINT=0. C KMAH... GIVES A NUMBER OF CHANGES OF SIGN OF THE QUANTITY C Q ALONG THE RAY (IN OTHER WORDS: IT SPECIFIES HOW C MANY TIMES THE RAY TOUCHED A CAUSTIC). IT IS USED C TO DETERMINE THE PHASE SHIFT CAUSED BY CAUSTICS. C C THE NUMERICAL CODE OF THE WAVE IS TRANSMITTED TO THE ROUTINE C JPAR BY COMMON/COD/. C ROUTINES USED: RK,FCTA. C NOTE THAT ALL THE INFORMATION ABOUT THE MODEL USED IN JPAR,RK C AND FCTA IS TRANSMITTED TO THESE ROUTINES ONLY BY COMMON/RAY/. C IN OTHER WORDS, THE DYNAMIC RAY TRACING PACKAGE CAN ALSO BE C USED WITH OTHER RAY TRACING ROUTINES, WHEN THEY GENERATE C COMMON/RAY/ AS DESCRIBED ABOVE. C C C **************************************************************** C C A SHORT DESCRIPTION OF THE UNIVERSAL ROUTINE FOR THE EVALUATION C OF RAY AMPLITUDES OF SEISMIC BODY WAVES ALONG A KNOWN RAY C AMPL(FJ,AX,AY,AZ,PHX,PHY,PHZ). C *********************** C C THE ROUTINE AMPL COMPUTES THE RAY AMPLITUDES AND PHASE SHIFTS C OF ARBITRARY MULTIPLY-REFLECTED (POSSIBLY CONVERTED) SEISMIC C BODY WAVES IN A 2-D LATERALLY INHOMOGENEOUS MEDIA WITH CURVED C INTERFACES. THE COMPUTATION IS PERFORMED ALONG A KNOWN RAY, AND C ALL THE INFORMATION ABOUT THE RAY IS TRANSFERRED TO AMPL BY C COMMON/RAY/. IT IS ALSO ASSUMED THAT THE GEOMETRICAL SPREADING C IS KNOWN (E.G., IT COULD BE DETERMINED INDEPENDENTLY BY THE PRO- C CEDURE JPAR, SEE ABOVE). THE GEOMETRICAL SPREADING MAY CORRES- C TO A POINT SOURCE OR A LINE SOURCE. C THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS: C C FJ... GEOMETRICAL SPREADING. C AX,AY,AZ... AMPLITUDES OF RADIAL, TRANSVERSE AND VER- C TICAL COMPONENTS. C PHX,PHY,PHZ...PHASE SHIFTS OF RADIAL, TRANSVERSE AND C VERTICAL COMPONENTS. C C THE NUMERICAL CODE OF THE WAVE IS TRANSFERRED TO THE ROUTINE C AMPL BY COMMON/COD/. C ROUTINES USED: COEF8. C NOTE THAT THE MODEL INFORMATION USED IN AMPL IS TRANSFERRED TO C THIS ROUTINE BY COMMON/RAY/ AND /SOUR/. IN OTHER WORDS, ROUTINE C AMPL CAN BE USED WITH ANY OTHER RAY TRACING ROUTINE, WHEN THE C ROUTINE GENERATES THESE COMMONS, AND WHEN THE GEOMETRICAL SPREA- C DING IS ALSO DETERMINED IN ADVANCE. C C C ***************************************************************** C C C A SHORT DESCRIPTION OF THE UNIVERSAL ROUTINE FOR DETERMINATION C OF REFLECTION/TRANSMISSION COEFFICIENTS C COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH). C *************************************************** C C THE ROUTINE COEF8 IS A STANDARD ROUTINE FOR THE COMPUTATION OF C P, SV AND SH REFLECTION/TRANSMISSION COEFFICIENTS OF PLANE WAVES C AT A PLANE INTERFACE BETWEEN TWO HOMOGENEOUS MEDIA OR AT A FREE C SURFACE OF A HOMOGENEOUS HALFSPACE. IN THE SECOND CASE, THE CON- C VERSION COEFFICIENTS AT RECEIVERS CAN ALSO BE OPTIONALLY COMPUTED. C CONVERSION COEFFICIENTS (COMPONENTS OF A CONVERSION VECTOR) ARE C DETERMINED WITH RESPECT TO A LOCAL CARTESIAN COORDINATE SYSTEM, C DEFINED AS FOLLOWS: POSITIVE X-AXIS POINTS ALONG THE PROFILE TO C THE RIGHT, POSITIVE Z-AXIS POINTS UP INTO THE FREE SPACE, THE C Y-AXIS IS CHOSEN SUCH AS TO MAKE THE SYSTEM RIGHT-HANDED (IT C POINTS BEHIND THE PROFILE). C C THE CODES OF INDIVIDUAL COEFFICIENTS ARE SPECIFIED BY THE C FOLLOWING NUMBERS C A/ INTERFACE BETWEEN TWO SOLID HALFSPACES C P1P1....1 P1SV1....2 P1P2....3 P1SV2....4 SH1SH1...9 C SV1P1...5 SV1SV1...6 SV1P2...7 SV1SV2...8 SH1SH2..10 C B/ FREE SURFACE (FOR RO2.LT.0.00001) C PP.....1 PX.....5 PX,PZ...X- AND Z- COMPONENTS OF THE C PSV....2 PZ.....6 COEF.OF CONVERSION,INCIDENT P WAVE. C SVP....3 SVX....7 SVX,SVZ...X- AND Z- COMPONENTS OF C SVSV...4 SVZ....8 COEF.OF CONVERSION,INCIDENT SV WAVE. C SHY...Y-COMPONENT OF SH COEF.OF CONV. C C I N P U T P A R A M E T E R S C P...RAY PARAMETER C VP1,VS1,RO1...PARAMETERS OF THE FIRST HALFSPACE C VP2,VS2,RO2...PARAMETERS OF SECOND HALFSPACE. FOR THE FREE C SURFACE TAKE RO2.LT.0.00001,EG.RO2=0., AND C ARBITRARY VALUES OF VP2 AND VS2 C NCODE...CODE OF THE COMPUTED COEFFICIENT C ND...=0 WHEN THE POSITIVE DIRECTION OF THE RAY C AND THE X-AXIS MAKE AN ACUTE ANGLE C =1 WHEN THE WAVE IMPINGES ON THE INTERFACE C AGAINST THE POSITIVE DIRECTION OF THE X-AXIS C C O U T P U T P A R A M E T E R S C RMOD,RPH...MODUL AND ARGUMENT OF THE COEFFICIENT C C N O T E S C 1/ POSITIVE P...IN THE DIRECTION OF PROPAGATION. C 2/ POSITIVE SV..TO THE LEFT FROM P. C 3/ TIME FACTOR OF INCIDENT WAVE ... EXP(-I*OMEGA*T) C 4/ FORMULAE ARE TAKEN FROM CERVENY ,MOLOTKOV, PSENCIK, RAY METHOD C IN SEISMOLOGY, PAGES 30-35. DUE TO THE ABOVE DEFINITION OF C THE LOCAL CARTESIAN COORDINATE SYSTEM, IN WHICH THE CONVERSION C COEFFICIENTS ARE EVALUATED, THE SIGNS OF CERTAIN COEFFICIENTS C ARE OPPOSITE. C C C *************************************************************** C C A SHORT DESCRIPTION OF THE ROUTINE FOR THE GENERATION OF NUMERICAL C CODES OF ELEMENTARY WAVES, C GENER(N,IDP,IDS,IBP,IBS,IUP,IUS,IREAD,IS,IR,NS,IFIN,NCODE). C ************************************************************ C C C THE ROUTINE GENER GENERATES AUTOMATICALLY THE NUMERICAL CODES C OF DIRECT AND PRIMARY REFLECTED P, S AND CONVERTED WAVES. THE C CONVERSIONS CAN TAKE PLACE ONLY AT THE REFLECTION POINT. IN C ADDITION, ANY OTHER MULTIPLY REFLECTED WAVES (POSSIBLY CONVER- C TED ARBITRARY NUMBER OF TIMES) CAN BE ADDED MANUALLY, USING C IREAD.NE.0. C THE SOURCE MAY BE SITUATED IN AN ARBITRARY LAYER (NO. IS), AND C SIMILARLY THE RECEIVER MAY BE ALSO SITUATED IN AN ARBITRARY C LAYER (NO.IR). C THE GENERATION STARTS WITH NAWX=0. THEN THE NUMERICAL CODE OF C THE FIRST WAVE IS GENERATED. AFTER THIS, NAWX MUST BE PUT =1 C IN THE CALLING PROGRAM, AND ROUTINE GENERATES SUCCESSIVELY NEW C CODES AT EACH CALLING. C THE SYSTEM OF GENERATED WAVES IS CONTROLLED BY THE PARAMETERS C IDP,IDS,IBP,IBS,IUP,IUS. FOR IDP,IDS,IBP,IBS, SEE EXPLANATIONS C IN THE DESCRIPTION OF THE INPUT DATA CARD NO.9. THE INDICES C IUP AND IUS CONTROL THE GENERATION OF PRIMARY REFLECTIONS WITH C THE REFLECTION POINTS ABOVE THE SOURCE AND THE RECEIVER. C THE QUANTITY NCODE IS THE INTERNAL WAVE CODE, WHICH IS INTRODUCED C IN ADDITION TO THE FULL NUMERICAL CODE (CALLED HERE ALSO EXTER- C NAL WAVE CODE).THE INTERNAL WAVE CODE IS JUST ONE INTEGER, WHICH C GIVES THE SUCCESSIVE NUMBER OF THE ELEMENTARY WAVE, AS IT WAS C GENERATED BY THE ROUTINE GENER. THE FIRST GENERATED WAVE HAS IN- C TERNAL CODE 1, ETC. THIS INTERNAL WAVE CODE HAS AN IMPORTANCE C IN THE COMPUTATION OF SYNTHETIC SEISMOGRAMS IN THE PROGRAM SYNTPL, C AS ELEMENTARY WAVES CAN BE SELECTED FOR THE CONSTRUCTION OF THESE C SEISMOGRAMS. C THE INDEX NLAYER GIVES A FULL NUMBER OF LAYERS OF THE MODEL FOR C WHICH THE NUMERICAL CODES OF ELEMENTARY WAVES ARE GENERATED. C IT MIGHT BE, OF COURSE, LESS THAN THE REAL NUMBER OF LAYERS IN C THE MODEL USED FOR THE COMPUTATION. FOR EXAMPLE, THE INDEX NLAY C (SEE INPUT CARD NO.9) CAN BE USED TO DECREASE ARTIFICIALLY THE C NUMBER OF LAYERS IN THE MODEL USED FOR GENERATION. C FOR IREAD=0, ONLY THE AUTOMATICALLY GENERATED DIRECT AND C PRIMARY REFLECTED WAVES ARE USED. FOR IREAD=1, ARBITRARY NUMBER C OF OTHER ELEMENTARY WAVES (MULTIPLE REFLECTIONS) CAN BE OPTIONALLY C ADDED, SEE THE DESCRIPTION OF MANUAL GENERATION IN INPUT CARD C NO.9 AND 13. C THE PARAMETER IFIN=0 WHEN THE ROUTINE RETURNS A NEW NUMERICAL CODE. C THE END OF GENERATION IS SHOWN BY IFIN=1. C THE REFRACTED WAVES ARE AUTOMATICALLY INCLUDED IN THE NUMERICAL C CODES CORRESPONDING TO REFLECTED WAVES, AS WAVES WITH COMPOUND C ELEMENTS, SEE DETAILS IN INPUT DATA CARD NO.13. C NOTE THAT IN THE LAYER CONTAINING THE SOURCE, THE WAVES, WHICH C STRIKE FIRST THE INTERFACE ABOVE THE SOURCE ARE COMPUTED SEPARA- C TELY FROM THOSE, WHICH STRIKE FIRST THE INTERFACE BELOW THE C SOURCE. C THE ROUTINE IS QUITE UNIVERSAL, IT CAN BE USED WITH ANY OTHER C PROGRAM WHICH WORKS WITH THE SAME NUMERICAL CODES OF ELEMENTARY C WAVES AND IN WHICH THE INTERFACES CROSS THE WHOLE MODEL, FROM C THE LEFT TO THE RIGHT. C C C **************************************************************** C C C A SHORT DESCRIPTION OF THE ROUTINE FOR THE APPROXIMATION OF C INTERFACES AND VELOCITY DISTRIBUTION WITHIN INDIVIDUAL LAYERS C MODEL(MTEXT,LU,MM). C ******************* C C *********************************************************** C BICUBIC SPLINE INTERPOLATION OF VELOCITIES SPECIFIED IN THE C GRID POINTS OF A RECTANGULAR NETWORK: MM=0 C *********************************************************** C C THE ROUTINE MODEL CONSISTS OF TWO PRINCIPAL PARTS: C 1) IN THE FIRST PART, THE INPUT DATA CONCERNING INTERFACES C AND VELOCITY DISTRIBUTION INSIDE INDIVIDUAL LAYERS ARE C READ IN. INTERFACES ARE INTERPOLATED BY CUBIC SPLINES C (ROUTINE SPLIN), THE VELOCITY DISTRIBUTION BY BICUBIC C SPLINES (ROUTINES BIAP AND SPLIN). THE COEFFICIENTS SO C OBTAINED ARE STORED IN COMMONS /INTRF/, /MEDIM/ AND C /VCOEF/. THE INPUT DATA CONCERNING THE DENSITY AND QUA- C LITY FACTOR DISTRIBUTION ARE STORED IN COMMON/DEN/ AND C COMMON/ABSR/. C C 2) IN THE SECOND PART, CERTAIN TABLES DESCRIBING THE MODEL ARE C OPTIONALLY PRINTED ON THE LINE PRINTER, SEE THE DETAILED C DESCRIPTION IN THE SECTION ON OUTPUT TABLES. C C THE MEANING OF THE INPUT PARAMETERS MTEXT, LU AND MM IS EXPLAI- C NED IN INPUT DATA CARD NO.1. NOTE THAT MM=0 FOR THIS ROUTINE. C ROUTINES USED: BIAP, SPLIN, VELOC. THE ROUTINE VELOC (SEE ITS C DETAILED DESCRIPTION FOLLOWING) IS USED ONLY IN THE SECOND C PART OF THE PROGRAM, FOR THE CONSTRUCTION OF THE PRINTER PLOT. C THE QUANTITIES STORED IN COMMON/AUXX/AND THE QUANTITIES LAY, C ITYPE,NDER FROM COMMON/AUXI/ ARE USED FOR COMMUNICATION BET- C WEEN THE ROUTINES MODEL, BIAP, SPLIN, VELOC, SEE THE DESCRIP- C TION OF COMMONS /AUXI/ AND /AUXX/. C C C **************************************************************** C C C A SHORT DESCRIPTION OF THE ROUTINE FOR THE DETERMINATION OF C VELOCITY AND ITS DERIVATIVES AT AN ARBITRARY POINT OF THE MEDIUM C VELOC(Y,VEL). C ************* C C C *********************************************************** C BICUBIC SPLINE INTERPOLATION OF VELOCITIES SPECIFIED IN THE C GRID POINTS OF A RECTANGULAR NETWORK: MM=0 C *********************************************************** C C IN THIS ROUTINE, THE SPECIFIED POINT AT WHICH THE VELOCITY C IS TO BE DETERMINED IS FIRST LOCALIZED IN THE GIVEN VELO- C CITY NETWORK. THEN THE VALUES OF VELOCITY AND OF ITS FIRST C PARTIAL DERIVATIVES ARE DETERMINED BY EVALUATING BICUBIC C SPLINE POLYNOMIALS. WHEN MDIM.GT.1, THE SECOND DERIVATIVES C ARE ALSO DETERMINED. C THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS: C C Y(1),Y(2)...X AND Z COORDINATES OF THE POINT AT WHICH THE C VELOCITY IS TO BE DETERMINED. C Y(3),Y(4)...VALUES OF THE COMPONENTS OF THE SLOWNESS VECTOR. C THEY HAVE ONLY A FORMAL MEANING IN THE ROUTINE. C VEL(1)... VELOCITY AT THE POINT Y(1),Y(2). C VEL(2),VEL(3)... FIRST PARTIAL DERIVATIVES VX AND VZ C VEL(4),VEL(5),VEL(6)... SECOND PARTIAL DERIVATIVES VXX, VXZ C AND VZZ. C C NOTE THAT FOR THE PROPER OPERATION OF ROUTINE VELOC, ALL THE C PARAMETERS OF THE COMMON/AUXX/ MUST BE SPECIFIED A PRIORI C (SEE THE DESCRIPTION OF COMMON/AUXX/). ALSO THE PARAMETERS C LAY,ITYPE AND NDER,(SEE THE DESCRIPTION OF COMMON/AUXI/), MUST C BE SPECIFIED BEFORE THE USE OF ROUTINE VELOC. C C ************************************************************** C C C A SHORT DESCRIPTION OF THE ROUTINE FOR THE APPROXIMATION OF C INTERFACES AND VELOCITY DISTRIBUTION WITHIN INDIVIDUAL LAYERS C MODEL(MTEXT,LU,MM). C ******************* C C *************************************************** C LINEAR VERTICAL INTERPOLATION OF VELOCITIES BETWEEN C ISOVELOCITY INTERFACES: MM=1 C *************************************************** C C THE ROUTINE MODEL CONSISTS OF TWO PRINCIPAL PARTS: C 1) IN THE FIRST PART, THE INPUT DATA CONCERNING INTERFACES C AND VELOCITY DISTRIBUTION INSIDE INDIVIDUAL LAYERS ARE C READ IN. INTERFACES ARE INTERPOLATED BY CUBIC SPLINES C (ROUTINE SPLIN). THE COEFFICIENTS THUS OBTAINED ARE STORED C IN COMMON /INTRF/. FOR EACH LAYER TWO VELOCITIES ARE SPECI- C FIED, ONE CONSTANT ALONG THE UPPER INTERFACE, THE OTHER C CONSTANT ALONG THE LOWER INTERFACE. THEY ARE STORED IN C COMMON/MEDIM/. THE INPUT DATA CONCERNING THE DENSITY C AND QUALITY FACTOR DISTRIBUTION ARE STORED IN COMMON/DEN/ C AND COMMON/ABSR/. C 2) IN THE SECOND PART, CERTAIN TABLES DESCRIBING THE MODEL ARE C OPTIONALLY PRINTED ON THE LINE PRINTER, SEE THE DETAILED C DESCRIPTION IN THE SECTION ON OUTPUT TABLES. C C THE MEANING OF THE INPUT PARAMETERS MTEXT, LU AND MM IS EXPLAI- C NED IN INPUT DATA CARD NO.1. LU, HOWEVER, HAS ONLY FORMAL MEA- C NING HERE. INTPUT DATA IN THIS ROUTINE ARE ALWAYS READ IN FROM C CARDS OR TERMINAL. NOTE THAT MM=1 FOR THIS ROUTINE. C ROUTINES USED: SPLIN, VELOC. THE ROUTINE VELOC (SEE ITS C DETAILED DESCRIPTION FOLLOWING) IS USED ONLY IN THE SECOND C PART OF THE PROGRAM, FOR THE CONSTRUCTION OF THE PRINTER PLOT. C THE QUANTITIES LAY, ITYPE,NDER FROM COMMON/AUXI/ ARE USED C FOR COMMUNICATION BETWEEN THE ROUTINES MODEL, SPLIN, AND VELOC, C SEE THE DESCRIPTION OF COMMON /AUXI/. C C C ************************************************************** C C C A SHORT DESCRIPTION OF THE ROUTINE FOR THE DETERMINATION OF C VELOCITY AND ITS DERIVATIVES AT AN ARBITRARY POINT OF THE C MEDIUM VELOC(Y,VEL). C ********************* C C C *************************************************** C LINEAR VERTICAL INTERPOLATION OF VELOCITIES BETWEEN C ISOVELOCITY INTERFACES: MM=1 C *************************************************** C C IN THIS ROUTINE, THE VALUES OF VELOCITY AND OF ITS FIRST C PARTIAL DERIVATIVES ARE DETERMINED BY INTERPOLATING LINEARLY C BETWEEN VELOCITY VALUES SPECIFIED ALONG THE UPPER AND LOWER C INTERFACE OF THE GIVEN LAYER. THE INTERPOLATION IS PERFORMED C ALONG A VERTICAL LINE PASSING THROUGH THE POINT IN WHICH C VELOCITY IS TO BE DETERMINED. WHEN MDIM.GT.1, THE SECOND C DERIVATIVES ARE ALSO DETERMINED. C THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS C FOLLOWS: C C Y(1),Y(2)...X AND Z COORDINATES OF THE POINT AT WHICH THE C VELOCITY IS TO BE DETERMINED. C Y(3),Y(4)...VALUES OF THE COMPONENTS OF THE SLOWNESS VECTOR. C THEY HAVE ONLY A FORMAL MEANING IN THE ROUTINE. C VEL(1)... VELOCITY AT THE POINT Y(1),Y(2). C VEL(2),VEL(3)... FIRST PARTIAL DERIVATIVES VX AND VZ C VEL(4),VEL(5),VEL(6)... SECOND PARTIAL DERIVATIVES VXX, VXZ C AND VZZ. C C NOTE THAT FOR THE PROPER OPERATION OF ROUTINE VELOC, THE C PARAMETERS LAY, ITYPE AND NDER,(SEE THE DESCRIPTION OF COM- C MON/AUXI/), MUST BE SPECIFIED BEFORE THE USE OF ROUTINE VELOC. C C ************************************************************** C C C C A SHORT DESCRIPTION OF THE ROUTINE FOR THE APPROXIMATION OF C INTERFACES AND VELOCITY DISTRIBUTION WITHIN INDIVIDUAL LAYERS C MODEL(MTEXT,LU,MM). C ******************* C C ********************************************************** C BILINEAR INTERPOLATION OF VELOCITIES SPECIFIED IN THE GRID C POINTS OF A RECTANGULAR NETWORK: MM=2 C ********************************************************** C C THE ROUTINE MODEL CONSISTS OF TWO PRINCIPAL PARTS: C 1) IN THE FIRST PART, THE INPUT DATA CONCERNING INTERFACES C AND VELOCITY DISTRIBUTION INSIDE INDIVIDUAL LAYERS ARE C READ IN. INTERFACES ARE INTERPOLATED BY CUBIC SPLINES C (ROUTINE SPLIN), THE VELOCITY DISTRIBUTION BY BILINEAR C INTERPOLATION (ROUTINE BILIN). THE COEFFICIENTS SO OBTAI- C NED ARE STORED IN COMMONS /INTRF/, /MEDIM/ AND /VCOEF/. C THE INPUT DATA CONCERNING THE DENSITY AND QUALITY FACTOR C DISTRIBUTION ARE STORED IN COMMON/DEN/ AND COMMON/ABSR/. C C C 2) IN THE SECOND PART, CERTAIN TABLES DESCRIBING THE MODEL ARE C OPTIONALLY PRINTED ON THE LINE PRINTER,SEE THE DETAILED C DESCRIPTION IN THE SECTION ON OUTPUT TABLES. C C THE MEANING OF THE INPUT PARAMETERS MTEXT, LU AND MM IS EXPLAI- C NED IN INPUT DATA CARD NO.1. NOTE THAT MM=2 FOR THIS ROUTINE. C ROUTINES USED: BIAP, SPLIN, VELOC. THE ROUTINE VELOC (SEE ITS C DETAILED DESCRIPTION FOLLOWING) IS USED ONLY IN THE SECOND C PART OF THE PROGRAM, FOR THE CONSTRUCTION OF THE PRINTER PLOT. C THE QUANTITIES STORED IN COMMON/AUXX/AND THE QUANTITIES LAY, C ITYPE,NDER FROM COMMON/AUXI/ ARE USED FOR COMMUNICATION BETWEEN C THE ROUTINES MODEL, BIAP, SPLIN, VELOC, SEE THE DESCRIPTION OF C COMMONS /AUXI/ AND /AUXX/. C C **************************************************************** C C C A SHORT DESCRIPTION OF THE ROUTINE FOR THE DETERMINATION OF C VELOCITY AND ITS DERIVATIVES AT AN ARBITRARY POINT OF THE MEDIUM C VELOC(Y,VEL). C ************* C C C ********************************************************** C BILINEAR INTERPOLATION OF VELOCITIES SPECIFIED IN THE GRID C POINTS OF A RECTANGULAR NETWORK: MM=2 C ********************************************************** C C IN THIS ROUTINE, THE SPECIFIED POINT AT WHICH THE VELOCITY C IS TO BE DETERMINED IS FIRST LOCALIZED IN THE GIVEN VELO- C CITY NETWORK. THEN THE VALUES OF VELOCITY AND OF ITS FIRST C PARTIAL DERIVATIVES ARE DETERMINED BY EVALUATING BILINEAR C POLYNOMIALS. WHEN MDIM.GT.1, THE SECOND DERIVATIVES ARE C ALSO DETERMINED. C THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS: C C Y(1),Y(2)...X AND Z COORDINATES OF THE POINT AT WHICH THE C VELOCITY IS TO BE DETERMINED. C Y(3),Y(4)...VALUES OF THE COMPONENTS OF THE SLOWNESS VECTOR. C THEY HAVE ONLY A FORMAL MEANING IN THE ROUTINE. C VEL(1)... VELOCITY AT THE POINT Y(1),Y(2). C VEL(2),VEL(3)... FIRST PARTIAL DERIVATIVES VX AND VZ C VEL(4),VEL(5),VEL(6)... SECOND PARTIAL DERIVATIVES VXX, VXZ C AND VZZ. C C NOTE THAT FOR THE PROPER OPERATION OF ROUTINE VELOC, ALL THE C PARAMETERS OF THE COMMON/AUXX/ MUST BE SPECIFIED A PRIORI C (SEE THE DESCRIPTION OF COMMON/AUXX/). ALSO THE PARAMETERS C LAY,ITYPE AND NDER,(SEE THE DESCRIPTION OF COMMON/AUXI/), MUST C BE SPECIFIED BEFORE THE USE OF ROUTINE VELOC. C C ************************************************************** C