******************************************************************

P R O G R A M   S E I S

******************************************************************

NUMERICAL MODELLING OF SEISMIC WAVE FIELDS
IN 2-D LATERALLY VARYING LAYERED STRUCTURES BY THE RAY METHOD


                           1)              2)
AUTHORS: VLASTISLAV CERVENY    IVAN PSENCIK
*******

ADDRESSES:     1) INSTITUTE OF GEOPHYSICS
*********         CHARLES UNIVERSITY
                  KE KARLOVU 3
                  121 16 PRAHA 2
                  CZECH REPUBLIC

               2) INSTITUTE OF GEOPHYSICS
                  ACAD. SCI. OF THE CR
                  BOCNI II
                  141 31 PRAHA 4
                  CZECH REPUBLIC

******************************************************************

DESCRIPTION OF THE PROGRAM.
***************************

  PROGRAM SEIS IS DESIGNED FOR THE COMPUTATION OF RAYS OF SEIS-
MIC WAVES WHICH ARRIVE AT A SYSTEM OF RECEIVERS DISTRIBUTED REGU-
LARLY OR IRREGULARLY ALONG THE EARTH'S SURFACE, ALONG AN INTERFACE,
OR ALONG A VERTICAL PROFILE. THE GENERATION OF WAVES IS SEMIAUTOMA-
TIC. AT RECEIVERS, CORRESPONDING TRAVEL TIMES ARE AUTOMATICALLY
DETERMINED. OPTIONALLY, AMPLITUDES AND PHASE SHIFTS MAY BE EVALU-
ATED. EFFECTS OF A SLIGHT ABSORPTION MAY BE ALSO CONSIDERED. ALL
THESE QUANTITIES ARE STORED AND MAY BE OPTIONALLY PLOTTED OR USED
FOR THE COMPUTATION OF SYNTHETIC SEISMOGRAMS.
  THE MODEL IS 2-D, LATERALLY INHOMOGENEOUS, WITH CURVED INTER-
FACES. INTERFACES ARE SPECIFIED BY POINTS READ FROM THE INPUT DA-
TA. THEY ARE APPROXIMATED BY CUBIC SPLINE INTERPOLATION. EACH
INTERFACE CROSSES THE WHOLE MODEL FROM ITS LEFT BORDER TO THE
RIGHT BORDER. THE INTERFACES MAY HAVE CORNER POINTS AND MAY BE
FICTITIOUS IN CERTAIN PARTS. VARIOUS INTERFACES MAY ALSO PARTIAL-
LY COINCIDE. THUS MODELS WITH VANISHING LAYERS, BLOCK STRUCTURES,
FRACTURES, ISOLATED BODIES MAY BE HANDLED BY THE PROGRAM. WITHIN
INDIVIDUAL LAYERS, THE VELOCITY MAY VARY BOTH VERTICALLY AND HO-
RIZONTALLY.
  THE SOURCE MAY BE SITUATED AT ANY POINT OF THE MEDIUM, WITH
EXCEPTION OF LAYERS OF ZERO THICKNESS, I.E., IT CANNOT BE SITU-
ATED BETWEEN TWO COINCIDING INTERFACES.
  ALL THE DIRECT AND PRIMARY REFLECTED WAVES P AND S, INCLUDING
THE CONVERTED WAVES AT THE POINT OF REFLECTION, CAN BE GENERATED
AUTOMATICALLY. MULTIPLE REFLECTIONS OF ARBITRARY TYPE ARE
OPTIONALLY GENERATED MANUALLY (BY INPUT DATA). THE REFRACTED
WAVES ARE CONSIDERED AS SPECIAL CASES OF REFLECTED WAVES WITH
COMPOUND RAY ELEMENTS, SEE A MORE DETAILED DESCRIPTION LATER.
  THE DETERMINATION OF RAYS WHICH ARRIVE AT SPECIFIED RECEIVERS
SITUATED EITHER ON SURFACE OR ON AN INTERFACE, OR ON A VERTICAL
PROFILE - TWO-POINT RAY TRACING - IS PERFORMED BY THE MODIFIED
SHOOTING METHOD. FOR ITERATIONS TO A CONSIDERED RECEIVER, METHOD
OF HALVING OF INTERVALS, REGULA FALSI METHOD, OR THE COMBINATION
OF THESE METHODS MAY BE USED. THE ITERATIONS TO A RECEIVER CONTI-
NUE UNLESS A RAY ARRIVING TO A 'REPS' VICINITY (FOR REPS SEE
INPUT DATA CARD NO.11) OF THE RECEIVER IS FOUND. THE ARRIVAL TIME
AT THE RECEIVER IS THEN OBTAINED BY A LINEAR INTERPOLATION FROM
ARRIVAL TIMES OF RAYS TERMINATING CLOSEST TO THE RAY AND SITUATED
TO BOTH SIDES FROM THE RECEIVER. THE AMPLITUDE CORRESPONDING TO
THE RAY CLOSEST TO THE RECEIVER IS ATRIBUTED TO THE RECEIVER.
  A SPECIAL CARE IS DEVOTED TO CERTAIN SINGULAR SITUATIONS IN
THE RAY FIELD. SPECIAL BRANCHES OF THE PROGRAM SEEK RAYS OF RE-
FRACTED WAVES, WHICH ARE CLOSE TO THE RAYS OF HEAD WAVES, RAYS
OF REFLECTED AND REFRACTED WAVES IN VICINITIES OF SHADOW BOUN-
DARIES, ETC. IN CASE THAT A RAY WITH IND.NE.3 (I.E. A RAY WHICH
DOES NOT COME TO THE EARTH'S SURFACE - SEE THE DESCRIPTION OF
THE PARAMETER IND) IS OBTAINED DURING THE SHOOTING, THE COR-
RESPONDING BOUNDARY RAY SEPARATING ILLUMINATED REGION (WITH
IND.EQ.3) FROM THE REGION WITH IND.NE.3, IS SOUGHT. AS SOON AS
THE BOUNDARY RAY IS FOUND THE ITERATIONS TO THE RECEIVERS IN
THE ILLUMINATED REGION ARE PERFORMED IN THE ABOVE DESCRIBED WAY.
THE PROGRAM WORKS IN A SINGLE PRECISION, SO THAT THE DETERMINA-
TION OF ALL WAVES IS IMPOSSIBLE IN CERTAIN SITUATIONS.
  THE AMPLITUDES OF P, SV AND SH WAVES ARE EVALUATED BY STANDARD
RAY FORMULAE. NO MODIFICATIONS ARE APPLIED IN SINGULAR REGIONS.
THE PURE HEAD WAVES GENERALLY DO NOT EXIST IN LATERALLY INHOMO-
GENEOUS MEDIA WITH CURVED INTERFACES AND ARE NOT CONSIDERED IN
SEIS. THEY MIGHT BE, OF COURSE,SIMULATED BY SLIGHTLY REFRACTED
WAVES.
  GEOMETRICAL SPREADING IS DETERMINED BY SOLVING A SYSTEM OF TWO
LINEAR ORDINARY DIFFERENTIAL EQUATIONS OF FIRST ORDER (SO CALLED
DYNAMIC RAY TRACING SYSTEM) BY A MODIFIED EULER'S METHOD.
  SLIGHT ABSORPTION MAY BE CONSIDERED. IN THIS CASE, A GLOBAL
ABSORPTION FACTOR (QUANTITY OBTAINED AS A TIME INTEGRAL OF 1/Q
ALONG A RAY, WHERE Q IS A QUALITY FACTOR) IS ALSO EVALUATED.
THE ABSORPTION IS CONSIDERED SLIGHT IF PRODUCT OF THE GLOBAL
ABSORPTION FACTOR WITH A CONSIDERED FREQUENCY IS APPROXIMATELY
EQUAL OR LESS THAN UNIT. QUALITY FACTOR DISTRIBUTION IN THE
MODEL IS RELATED TO THE DISTRIBUTION OF CORRESPONDING VELOCITY,
SEE BELOW.

PROGRAM SEIS IS SUPPLEMENTED BY PROGRAMS FOR PREPARATION OF
DATA FOR VELOCITY MODEL (PROGRAM SMOOTH), TO COMPUTE SYNTHETIC
SEISMOGRAMS (PROGRAM SYNTPL), AND TO PLOT RAYS, TRAVEL-TIME AND
AMPLITUDE-DISTANCE CURVES (PROGRAM RAYPLOT), SYNTHETIC SEIS-
MOGRAMS (PROGRAM SEISPLOT) OR PARTICLE MOTION DIAGRAMS (PROGRAM
POLARPLOT). THE PLOTTING PROGRAMS REQUIRE CALCOMP ROUTINES PLOTS,
PLOT, NUMBER AND SYMBOL.

IN THE PRESENT VERSION OF THE PROGRAM SEIS, THERE ARE THREE
POSSIBILITIES OF THE APPROXIMATION OF THE VELOCITY DISTRIBUTION
INSIDE INDIVIDUAL LAYERS: 1) BICUBIC SPLINE INTERPOLATION (IDEN-
TICAL TO THAT USED IN SEIS81), 2) LINEAR INTERPOLATION BETWEEN
ISOVELOCITY INTERFACES, 3) PIECE-WISE BILINEAR INTERPOLATION.
APPROXIMATIONS 1) AND 3) ARE CONVENIENT FOR MODELLING SPECIFIC
FEATURES WITHIN LAYERS. THE APPROXIMATION 2) IS CONVENIENT
ESPECIALLY WHEN THERE IS NO DETAILED INFORMATION ABOUT THE
VELOCITY DISTRIBUTION WITHIN A LAYER. NOTE THAT THE APPROXIMA-
TION 2) NEEDS A MINIMUM OF INPUT DATA. ALSO NOTE THAT THIS
APPROXIMATION IS CONVENIENT FOR MODELLING SECOND ORDER INTER-
FACES (WITH CONTINUOUS VELOCITY, BUT DISCONTINUOUS GRADIENT)
AND FOR MODELLING OF MEDIA CONTAINING HOMOGENEOUS LAYERS.
THE APPROXIMATION 3) IS CONVENIENT FOR KINEMATIC MODELLING.
IT SHOULD NOT BE USED FOR THE COMPUTATION OF AMPLITUDES. IT
IS EASY TO WRITE A NEW ROUTINE MODEL AND CORRESPONDING SUB-
ROUTINES FOR ANY TYPE OF VELOCITY APPROXIMATION WITHIN LAYERS.
THE ONLY CONDITIONS ARE THAT THE ROUTINE VELOC MUST BE SPE-
CIFIED AS DESCRIBED BELOW AND IT MUS CONTAIN COMMON/AUXI/.

  THE USERS ARE KINDLY ASKED TO REFER TO PROGRAM SEIS IN
ANY CASE THEY PRESENT SOME RESULTS OBTAINED WITH ITS HELP.



I N P U T   D A T A
*******************

INPUT DATA FOR THE MODEL
************************
THE 2-D MODEL IS SPECIFIED IN A RIGHT=HANDED CARTESIAN COORDI-
NATE SYSTEM X,Y,Z. THE X- AND Y-AXES ARE HORIZONTAL, THE Z-AXIS
IS VERTICAL, POSITIVE DOWNWARDS. THE MODEL IS SITUATED IN THE
(X,Z) PLANE SO THAT X COORDINATE INCREASES FROM THE LEFT TO THE
RIGHT. THE MODEL IS BOUNDED BY TWO VERTICAL BOUNDARIES ON ITS
LEFT-HAND AND  RIGHT-HAND SIDE AND BY THE FIRST AND THE LAST
INTERFACE AT THE TOP AND THE BOTTOM OF THE MODEL. THE TOP INTER-
FACE CORRESPONDS TO THE EARTH'S SURFACE.
THE INPUT DATA ARE READ IN FROM STANDARD INPUT BY LIST-DIRECTED
INPUT (FREE FORMAT):


1a)ONE CARD
  MTEXT
     MTEXT...  ARBITRARY ALPHANUMERIC TEXT DESCRIBING THE MODEL.
               THIS TEXT WILL APPEAR IN THE PLOTS OF SYNTHETIC
               SEISMOGRAMS.

1b)ONE CARD
  MPRINT,MM
     MPRINT... CONTROLS THE PRINTING OF THE DESCRIPTION OF THE
               MODEL ON THE LINE PRINTER.
               MPRINT=0... ONLY PRINT OF INPUT DATA FOR THE MO-
               DEL, NO OTHER DESCRIPTION OF THE MODEL IS PRINTED.
               MPRINT=1... ONLY A SIMPLE PRINTER PLOT OF THE
               VELOCITY DISTRIBUTION IN THE MODEL.
               MPRINT=2... MORE DETAILED DESCRIPTION OF THE
               MODEL.
     MM...     INDICATION OF THE TYPE OF THE ROUTINE MODEL USED
               FOR THE APPROXIMATION OF VELOCITY DISTRIBUTION.
               MM=0... BICUBIC SPLINE INTERPOLATION.
               MM=1... LINEAR VERTICAL INTERPOLATION BETWEEN
                       ISOVELOCITY INTERFACES.
               MM=2...  PIECE-WISE BILINEAR INTERPOLATION
               NOTE: IF THE VALUE OF MM DOES NOT CORRESPOND TO
               THE USED ROUTINE MODEL, THE COMPUTATION TERMINATES.

2) ONE CARD, DESCRIPTION OF INTERFACES.
   NINT,NPNT(1),NPNT(2),...,NPNT(NINT)
      NINT...   NUMBER OF INTERFACES,INCLUDING THE TOP AND THE
                BOTTOM BOUNDARIES OF THE MODEL (NINT.GE.2,
                NINT.LE.20). THE TOP BOUNDARY OF THE MODEL
                (FREE SURFACE)  AND THE BOTTOM BOUNDARY MAY
                ALSO BE CURVED.
                NOTE THAT THE NUMBER OF LAYERS IN THE MODEL IS
                NINT-1.
      NPNT(I)...NUMBER OF POINTS SPECIFYING THE I-TH INTERFACE.
                (NPNT(I).GE.2,NPNT(I).LE.30).
                NOTE: IN CASE OF MM=1 (SEE CARD NO.1), ALL INTERFA-
                CES REPRESENT ISOLINES OF VELOCITY. THE INTERFACES
                MAY BE OF THE FIRST ORDER WHEN THE VELOCITIES IMME-
                DIATELY ABOVE AND BELOW INTERFACE DIFFER, OR THEY
                MAY BE OF THE SECOND ORDER WHEN BOTH VELOCITIES ARE
                EQUAL.

3) NINT SYSTEMS OF CARDS, DESCRIPTIONS OF INDIVIDUAL INTERFACES.
   FOR EACH INTERFACE ONE SYSTEM OF CARDS. INTERFACES ARE TAKEN
   FROM THE TOP TO THE BOTTOM.
   EACH INTERFACE CROSSES THE WHOLE MODEL, FROM THE LEFT VERTI-
   CAL BOUNDARY TO THE RIGHT VERTICAL BOUNDARY.
   THE FIRST (TOP) INTERFACE REPRESENTS THE UPPER BOUNDARY OF THE
   MODEL (THE EARTH'S SURFACE).
   THE LAST (BOTTOM) INTERFACE REPRESENTS THE LOWER BOUNDARY OF
   THE MODEL.
   THE SYSTEM OF CARDS FOR THE I-TH INTERFACE CONSISTS OF NPNT(I)
   TRIPLETS OF NUMBERS; EACH TRIPLET CORRESPONDS TO ONE POINT OF
   THE INTERFACE (POINTS ARE TAKEN FROM THE LEFT TO THE RIGHT).
   EACH TRIPLET CONTAINS THE FOLLOWING THREE NUMBERS: A) THE X
   COORDINATE OF THE POINT, B) THE Z COORDINATE OF THE POINT,
   C) A SPECIAL INDEX III, SEE BELOW. THE TRIPLETS ARE READ IN
   INDEPENDENTLY FOR EACH INTERFACE.
   THE FIRST POINTS OF ALL INTERFACES MUST LIE ON ONE VERTICAL
   LINE, WHICH IS CALLED THE LEFT BOUNDARY OF THE MODEL.
   THE LAST POINTS OF ALL INTERFACES MUST LIE ON ANOTHER VERTI-
   CAL LINE, WHICH IS CALLED THE RIGHT BOUNDARY OF THE MODEL.
   THE PART OF THE INTERFACE BETWEEN TWO SUCCESSIVE POINTS ON
   THE INTERFACE IS CALLED THE ELEMENT OF THE INTERFACE.

   THE QUANTITY III CONTROLS THE CHARACTER OF THE INTERFACE IN
   THE NEIGHBOURHOOD OF THE CORRESPONDING POINT:
      III=0     THE INTERFACE ITSELF, AS WELL AS ITS FIRST AND
                SECOND DERIVATIVES, ARE SMOOTH AT THE CORRESPON-
                DING POINT.
      III=-1    THE CORNER POINT IN THE INTERFACE (DISCONTINUOUS
                FIRST DERIVATIVE).
      III=-2    THE CORNER POINT OF THE INTERFACE. MOREOVER, THE
                INTERFACE TO THE RIGHT SIDE OF THE POINT IS
                FICTITIOUS.
      III=K,K.GT.0   THE ELEMENT OF THE INTERFACE TO THE RIGHT
                OF THE CORRESPONDING POINT COINCIDES WITH THE
                K-TH ELEMENT OF THE ADJOINING UPPER INTERFACE,
                STARTING FROM THE K-TH POINT OF THE UPPER INTER-
                FACE. THE COORDINATES OF THE CORRESPONDING POINTS
                AT BOTH INTERFACES MUST BE THE SAME.


4) DESCRIPTION OF VELOCITY DISTRIBUTION IN INDIVIDUAL LAYERS.
   NOTE THAT THESE INPUT DATA MAY DIFFER FOR DIFFERENT KINDS
   OF APPROXIMATION OF THE VELOCITY DISTRIBUTION.

***************************************************
MM=   BICUBIC SPLINE OR PIECE-WISE BILINEAR INTERPOLATION
0,2   (MM=0 OR MM=2, SEE INPUT DATA CARD NO.1)
***************************************************

   (NINT-1) SYSTEMS OF CARDS, DESCRIPTIONS OF VELOCITY DISTRIBU-
   TION IN INDIVIDUAL LAYERS. LAYERS ARE TAKEN FROM THE TOP TO
   THE BOTTOM. STARTING FROM THE FIRST LAYER. THE DATA CAN BE
   READ IN EITHER FROM A CARD READER, TERMINAL, ETC., OR FROM
   A DATA SET STORED IN THE FILE LU (SEE DATA SUB 1). VELOCITY
   IS SPECIFIED AT GRID POINTS OF A RECTANGULAR NETWORK. THE
   NETWORK MAY BE DIFFERENT IN DIFFERENT LAYERS. THE NETWORK
   MUST ALWAYS COVER THE WHOLE CORRESPONDING LAYER. NOTE THAT
   IN THIS CASE THE INPUT DATA ARE IDENTICAL TO THOSE USED IN
   THE PROGRAM SEIS81.
   THE SYSTEM OF CARDS FOR EACH LAYER CONSISTS OF THE FOLLOWING
   CARDS:

4A) ONE CARD,NUMBER OF NETWORK LINES.
   MX,MY
      MX...     NUMBER OF VERTICAL LINES IN THE NETWORK.
      MY...     NUMBER OF HORIZONTAL LINES IN THE NETWORK.

4B) ONE CARD (OR GROUP OF CARDS), X COORDINATES OF VERTICAL NET-
    WORK LINES.
   SX(1),SX(2),...,SX(MX)
      SX(I)...  X COORDINATE OF THE I-TH VERTICAL LINE.
                SX(1) SPECIFIES THE LEFT-HAND BOUNDARY
                OF THE MODEL (THE SAME FOR ALL LAYERS).
                SX(MX) SPECIFIES THE RIGHT-HAND BOUNDARY
                OF THE MODEL (THE SAME FOR ALL LAYERS).

4C) ONE CARD (OR GROUP OF CARDS), Z COORDINATES OF HORIZONTAL
    NETWORK LINES.
   SY(1),SY(2),...,SY(MY)
      SY(I)...  Z COORDINATE OF THE I-TH HORIZONTAL LINE.
                SY(1) AND SY(MY) MUST BE CHOSEN IN SUCH A WAY
                THAT THE NETWORK COVERS THE WHOLE LAYER. THEY
                MAY BE COMPLETELY OUTSIDE THE LAYER.

4D) A GROUP OF CARDS, VELOCITY VALUES AT GRID POINTS ALONG THE
    VERTICAL NETWORK LINES. FOR EACH VERTICAL NETWORK LINE,'MY'
    VELOCITY VALUES ARE READ IN, INDEPENDENTLY FOR EACH LINE.
    THE GRID POINTS ALONG THE VERTICAL NETWORK
    LINE ARE TAKEN FROM THE TOP TO THE BOTTOM. THE VERTICAL NET-
    WORK LINES ARE TAKEN FROM LEFT TO RIGHT.

    NOTES ON MAXIMUM PERMITTED DIMENSIONS:
    THE SUM OF ALL VERTICAL LINES IN ALL LAYERS MAY NOT EXCEED
    350.
    THE SUM OF ALL HORIZONTAL LINES IN ALL LAYERS MAY NOT EXCEED
    350.
    THE TOTAL NUMBER OF ALL GRID POINTS IN ALL LAYERS (I.E., THE
    TOTAL NUMBER OF ALL VELOCITY VALUES READ IN) MAY NOT EXCEED
    1000.

*************************************************************


**********************************************************
MM=    LINEAR INTERPOLATION BETWEEN ISOVELOCITY INTERFACES
1     (MM=1, SEE INPUT DATA CARD NO.1)
**********************************************************

    THESE DATA ARE ALWAYS READ IN FROM A CARD READER OR TERMINAL.


4A) ONE CARD (OR A GROUP OF CARDS), VELOCITIES BELOW INTERFACES
   V1(1),V1(2),....,V1(NINT-1)
      V1(I)...    VELOCITY IMMEDIATELY BELOW I-TH INTERFACE,
                  I.E., ALONG THE UPPER BOUNDARY OF THE I-TH
                  LAYER.

4B) ONE CARD (OR A GROUP OF CARDS), VELOCITIES ABOVE INTERFACES
   V2(1),V2(2),...,V2(NINT-1)
      V2(I)...    VELOCITY IMMEDIATELY ABOVE (I+1)-TH INTERFACE,
                  I.E., ALONG THE LOWER BOUNDARY OF THE I-TH
                  LAYER.

5) ONE CARD, SPECIFICATION OF DENSITIES, S VELOCITIES AND QUALITY
   FACTORS.
   NRO, NVS(1), NVS(2),...,NVS(NINT-1), NABS
      NRO... NRO=0... THE DENSITY AT ANY POINT OF THE MEDIUM
                IS DETERMINED AUTOMATICALLY FROM THE RELATION
                DENSITY=1.7+0.2*VP (WHERE VP IS THE P VELOCITY).
                THE FOLLOWING CARD NO.6A IS NOT READ IN IN THIS
                CASE.
             NRO=1... THE DENSITY IN THE I-TH LAYER IS DETER-
                MINED FROM THE RELATION
                DENSITY=RHO1(I)+RHO2(I)*VP,
                WHERE RHO1(I) AND RHO2(I) ARE READ IN IN THE INPUT
                DATA CARD NO.6A.
                FOR A LIQUID LAYER, DENSITY=1 AUTOMATICALLY.
      NVS(I).. NVS(I)=0...THE VELOCITY VALUES READ IN SUB 4
                CORRESPOND TO P WAVE VELOCITY DISTRIBUTION
                IN THE I-TH LAYER. S WAVE VELOCITIES IN THE I-TH
                LAYER ARE DETERMINED FROM THE RELATION
                VS=VP/PTOS(I). THE VALUES OF PTOS(I) ARE READ IN
                IN THE INPUT CARD NO.7
               NVS(I)=1... THE VELOCITY VALUES READ IN SUB 4
                CORRESPOND TO THE S WAVE VELOCITY DISTRIBUTION IN
                THE I-TH LAYER. P WAVE VELOCITIES IN THIS LAYER
                ARE DETERMINED FROM THE RELATION
                VP=VS*PTOS(I), WHERE PTOS(I) ARE READ IN IN THE
                CARD NO.7.
                NOTE THAT THE FIRST LAYER MAY BE ALSO LIQUID
                (E.G. OCEAN). IN THIS CASE, ONLY NVS(1)=0 IS ALLO-
                WED, NVS(1)=1 IS NOT PERMITTED. PUT PTOS(1).GE.100
                IN THIS CASE, THEN THE S WAVE VELOCITY IN THE
                FIRST LAYER IS AUTOMATICALLY ZERO.
      NABS... NABS=0... NO ABSORPTION IS CONSIDERED.
              NABS=1... SLIGHT ABSORPTION IS CONSIDERED.
                FOR NABS=1, THE DATA SUB 6B ARE READ IN. THE
                QUALITY FACTORS QP AND QS ARE DETERMINED FROM
                ONE OF THE RELATIONS GIVEN IN DESCRIPTION OF
                INPUT DATA SUB 6B. NO SPECIFICATION OF THE TYPE
                OF THE ABSORPTION (FUTTERMAN'S OR NON-CAUSAL) IS
                MADE IN THIS PROGRAM. THE TYPE IS SPECIFIED LATER
                IN PROGRAMS RAYPLOT OR SYNTPL, INCLUDED IN THIS
                PACKAGE. HERE ONLY THE GLOBAL ABSORPTION FACTOR
                IS COMPUTED.

6A)ONE CARD, SPECIFICATION OF DENSITIES.
   NOTE: THIS CARD IS INCLUDED ONLY WHEN NRO=1.
   RHO1(1),RHO2(1),RHO1(2),RHO2(2),..,RHO1(NINT-1),RHO2(NINT-1)
      RHO1(I),RHO2(I)... PARAMETERS OF THE RELATION FOR THE
                DETERMINATION OF THE DENSITY FROM P WAVE VELO-
                CITIES IN THE I-TH LAYER,
                DENSITY=RHO1(I)+RHO2(I)*VP.

6B)(NINT-1) CARDS, SPECIFICATION OF QUALITY FACTORS.
   NOTE: THESE CARDS ARE INCLUDED ONLY IF NABS.NE.0.
         ONE CARD FOR EACH LAYER, FROM TOP TO THE BOTTOM.
         EACH CARD HAS A FORM:
   NQP,NQS,QP1,QP2,QP3,QS1,QS2,QS3
      NQP... NQP=0... THE QUALITY FACTOR QP IN THE CONSIDERED
                LAYER IS DETERMINED FROM THE RELATION:
                QP=QP1+QP2*VP+QP3*VP*VP,
                WHERE VP IS THE P WAVE VELOCITY;
             NQP=1... THE QUALITY FACTOR QP IN THE CONSIDERED
                LAYER IS DETERMINED FROM THE RELATION:
                QP=QP1+QP2*VP+QP3/VP.
      NQS... THE SAME AS FOR QP; IT IS ONLY NECESSARY TO RE-
                PLACE THE QUANTITIES VP,QP1,QP2,QP3 BY VS,QS1,
                QS2,QS3, WHERE VS IS THE S WAVE VELOCITY.
      QP1,QP2,QP3,QS1,QS2,QS3... COEFFICIENTS OF THE RELATIONS
                FOR THE DETERMINATION OF QUALITY FACTORS
                IN A CONSIDERED LAYER.

7) ONE CARD, SPECIFICATION OF THE RATIOS BETWEEN S AND P VELO-
   CITIES IN INDIVIDUAL LAYERS
   PTOS(1),PTOS(2),...,PTOS(NINT-1)
      PTOS(I)...RATIO OF THE P WAVE VELOCITY TO THE S WAVE VELO-
                CITY IN THE I-TH LAYER,VS=VP/PTOS(I).
                NOTE THAT PTOS(I) MUST BE ALWAYS GREATER THAN
                SQRT(2.). FOR PTOS(I).LT.SQRT(2.), E.G. WHEN
                PTOS(I) IS NOT SPECIFIED, PTOS(I)=1.732 AUTO-
                MATICALLY.
                FOR PTOS(I).GE.100., VS VELOCITY IN THE I-TH
                LAYER IS ZERO AUTOMATICALLY (LIQUID LAYER).

8) ONE CARD, CONTROLS THE PRINTER PLOT OF THE P VELOCITY DISTRI-
   BUTION IN THE MODEL, WHICH IS PRINTED ON A LINE PRINTER WHEN
   MPRINT.GE.1 (SEE CARD NO.1). THE PRINTER PLOT CONSISTS OF A
   SYSTEM OF DIGITS FROM 0 TO 9 AND ASTERISKS. THE DIGITS ARE
   DETERMINED FROM THE COMPUTED VELOCITY DISTRIBUTION USING
   THE RELATION I=IFIX(10.*(V-VMIN)/(VMAX-VMIN)). OUTSIDE THE
   RANGE VMIN, VMAX, THE ASTERISKS ARE PRINTED. THE CARD SHOULD
   BE FORMALLY INCLUDED EVEN WHEN MPRINT=0.

   VMIN,VMAX,BMIN,BMAX,BLEFT,BRIGHT
      VMIN,VMAX... MINIMUM AND MAXIMUM P VELOCITY SHOWN IN THE
                PRINTER PLOT
      BMIN,BMAX... Z COORDINATES OF THE TOP AND THE BOTTOM HORI-
                ZONTAL BOUNDARIES OF THE PRINTER PLOT.
                WHEN BMIN AND BMAX ARE NOT SPECIFIED,Z-COORDINATES
                OF THE FIRST POINTS OF THE TOP AND THE BOTTOM
                INTERFACES ARE TAKEN AS BMIN AND BMAX.
      BLEFT,BRIGHT... X COORDINATES OF THE VERTICAL BOUNDARIES
                OF THE PRINTER PLOT. THEY MAY DIFFER FROM VERTI-
                CAL BOUNDARIES OF THE MODEL. IF BLEFT AND BRIGHT
                ARE NOT SPECIFIED, THEN THEY ARE TAKEN AS VERTICAL
                BOUNDARIES OF THE MODEL.


SPECIFICATION OF PARAMETERS OF COMPUTED SEISMIC BODY WAVES

**********************************************************

9) ONE CARD, VARIOUS SWITCHES.
   ICONT,MEP,MOUT,MDIM,METHOD,IBP,IBS,IUP,IUS,IDP,IDS,IREAD,MREG,
   ITMAX,NLAY,LU1,LU2,LTAPE,KSH,ILOC
      ICONT...  CONTROLS THE CONTINUATION OF THE COMPUTATION.
             ICONT=0... TERMINATION OF COMPUTATION. THIS IS
                THE LAST CARD IN THE INPUT DATA CARD SET.
             ICONT=1... THE COMPUTATIONS CONTINUE, THE NEXT
                CARD SHOULD BE CARD NO 10.
             ICONT=-1... THE COMPUTATIONS CONTINUE FOR A NEW
                MODEL OF MEDIUM.THE NEXT CARD SHOULD BE THE CARD
                NO.1.
                IN CASE OF ICONT=0 OR ICONT=-1, THE REMAINING
                QUANTITIES IN THIS CARD MAY BE ARBITRARY.
      MEP...    ABS(MEP)... NUMBER OF RECEIVER POSITIONS. THE RE-
                CEIVERS CAN BE SITUATED EITHER ALONG THE EARTH'S
                SURFACE OR ALONG ONE OF THE INTERFACES, OR ALONG
                A VERTICAL PROFILE. IN THE FORMER TWO CASES, THE
                RECEIVERS ARE SPECIFIED BY A VARYING X COORDINATE,
                Z COORDINATES FOLLOW AUTOMATICALLY FROM THE EQUA-
                TION OF THE SURFACE OR INTERFACE. IN THE LATTER
                CASE, THE RECEIVERS ARE SPECIFIED BY A VARYING Z
                COORDINATE, X COORDINATES FOLLOW AUTOMATICALLY
                FROM THE SPECIFICATION OF THE VERTICAL PROFILE.
                THE SIGN OF MEP CONTROLS THE WAY IN WHICH THE
                SYSTEM OF RECEIVER POSITIONS IS SPECIFIED (SEE
                CARD NO.10).
             IF(MEP.GT.0),RECEIVERS ARE DISTRIBUTED REGULARLY.
                ONLY THE POSITIONS OF THE FIRST RECEIVER AND THE
                STEP IN THE VARYING COORDINATE ARE READ IN, SEE
                CARD NO.10.
             IF(MEP.LT.0), THE RECEIVERS ARE DISTRIBUTED IRRE-
                GULARLY ALONG THE PROFILE. THE VARYING COORDINATES
                OF ALL RECEIVER POSITIONS ARE READ IN, SEE CARD
                NO.10.
                ABS(MEP).LE.100.
      MOUT...   CONTROLS THE PRINTING OF THE OUTPUT ON A LINE
                PRINTER.
             MOUT=0... ONLY INPUT DATA ARE REPRODUCED ON A
                LINE PRINTER AND LISTING OF WAVE CODES IS PRINTED.
             MOUT=1... ONLY THE TABLE OF THE TERMINATION POINTS
                OF RAYS WHICH WERE SUCCESFULLY ITERATED AND TER-
                MINATE AT SPECIFIED RECEIVERS AND OF BOUNDARY RAYS
             MOUT=2... THE TABLE OF TERMINATION POINTS OF ALL
                RAYS COMPUTED IN THE PROCESS OF ITERATIONS.
      MDIM...   CONTROLS ON WHICH LEVEL AND WITH WHICH TYPE OF
                SOURCE THE PROGRAM SHOULD WORK.
             MDIM=0... ONLY THE RAYS AND TRAVEL TIMES ARE COM-
                PUTED, NO AMPLITUDES AND PHASE SHIFTS.
                FOR MDIM.GT.0, THE AMPLITUDES AND PHASE SHIFTS
                ARE COMPUTED.
                THE SPREADING MAY BE OPTIONALLY COMPUTED IN TWO
                DIFFERENT WAYS:
             MDIM=1... THE GEOMETRICAL SPREADING IS NOT CONSIDERED
                IN THE EVALUATION OF RAY AMPLITUDES.
             MDIM=2... LINE SOURCE APPROXIMATION
             MDIM=3... POINT SOURCE APPROXIMATION.
      METHOD... TO FIND THE RAYS WHICH HAVE TERMINATION POINTS
                AT A SPECIFIED SYSTEM OF RECEIVERS, THE SHOOTING
                METHOD IS USED. TO PERFORM ITERATIONS, VARIOUS
                METHODS CAN BE APPLIED:
             METHOD=0... REGULA FALSI.
             METHOD=1... COMBINATION OF REGULA FALSI METHOD
                AND METHOD OF HALVING INTERVALS.
             METHOD=2... HALVING OF INTERVALS.
             METHOD=3... METHOD BASED ON PARAXIAL RAY APPROXIMATION.
      IBP,IBS,IUP,IUS,IDP,IDS... SWITCHES WHICH CONTROL THE AUTO-
                MATIC GENERATION OF NUMERICAL CODES OF ELEMENTARY
                WAVES. ONLY DIRECT WAVES AND PRIMARY REFLECTED
                WAVES (POSSIBLY CONVERTED AT THE POINT OF REFLEC-
                TION) CAN BE GENERATED AUTOMATICALLY.
      IBP...    CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF-
                LECTED WAVES WHICH START AS A P WAVE FROM THE
                SOURCE. THE REFLECTION TAKES PLACE BELOW THE
                SOURCE.
             IBP=0... NO PRIMARY REFLECTED WAVES ARE GENERATED.
             IBP=1... PP PRIMARY REFLECTED WAVES ARE GENERATED.
             IBP=2... PS PRIMARY REFLECTED WAVES ARE ALSO GENE-
                RATED (TOGETHER WITH PP PRIMARY REFLECTED WAVES).
      IBS...    CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF-
                LECTED WAVES WHICH START AS AN S WAVE FROM THE
                SOURCE. THE REFLECTION TAKES PLACE BELOW THE SOURCE.
             IBS=0... NO PRIMARY REFLECTED WAVES ARE GENERATED.
             IBS=1... SS PRIMARY REFLECTED WAVES ARE GENERATED.
             IBS=2... SP PRIMARY REFLECTED WAVES ARE ALSO GENE-
                RATED (TOGETHER WITH SS PRIMARY REFLECTED WAVES).
      IUP...    CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF-
                LECTED WAVES WHICH START AS A P WAVE FROM THE
                SOURCE. THE REFLECTION TAKES PLACE ABOVE THE SOURCE.
             IUP=0... NO PRIMARY REFLECTED WAVES ARE GENERATED.
             IUP=1... PP PRIMARY REFLECTED WAVES ARE GENERATED.
             IUP=2... PS PRIMARY REFLECTED WAVES ARE ALSO GENE-
                RATED (TOGETHER WITH PP PRIMARY REFLECTED WAVES).
      IUS...    CONTROLS THE AUTOMATIC GENERATION OF PRIMARY REF-
                LECTED WAVES WHICH START AS AN S WAVE FROM THE
                SOURCE. THE REFLECTION TAKES PLACE ABOVE THE SOURCE.
             IUS=1... NO PRIMARY REFLECTED WAVES ARE GENERATED.
             IUS=2... SS PRIMARY REFLECTED WAVES ARE GENERATED.
             IUS=3... SP PRIMARY REFLECTED WAVES ARE ALSO GENE-
                RATED (TOGETHER WITH SS PRIMARY REFLECTED WAVES).
      IDP...    CONTROLS THE AUTOMATIC GENERATION OF THE DIRECT
                P WAVE.
             IDP=0... NO DIRECT P WAVE IS GENERATED.
             IDP=1... DIRECT P WAVE IS GENERATED.
      IDS...    CONTROLS THE AUTOMATIC GENERATION OF THE DIRECT
                S WAVE.
             IDS=0... NO DIRECT S WAVE IS GENERATED
             IDS=1... DIRECT S WAVE IS GENERATED.
      IREAD...  CONTROLS THE MANUAL GENERATION OF NUMERICAL CODES
                OF ELEMENTARY WAVES.
                IREAD=0... NO NUMERICAL CODES ARE MANUALLY GENERATED.
             IREAD=1... NUMERICAL CODES OF CERTAIN ELEMENTARY
                WAVES ARE MANUALLY GENERATED, SEE CARD NO. 13.
      MREG...   CONTROLS THE EVALUATION OF AMPLITUDES AT THE
                TERMINATION POINT OF THE RAY.
             MREG=0.. AT THE EARTH'S SURFACE, THE COEFFICIENTS
                OF CONVERSION ARE APPLIED (THEY TAKE INTO ACCOUNT
                THE EFFECTS OF THE EARTH'S SURFACE).
             MREG=1... AT THE EARTH'S SURFACE, THE COEFFICIENTS
                OF CONVERSION ARE NOT APPLIED, THE COMPONENTS OF
                THE DISPLACEMENT VECTOR ARE CALCULATED AS IF WITHIN
                AN INFINITE HALFSPACE.
                NOTE: WHEN THE TERMINATION POINT OF THE RAY IS
                SITUATED INSIDE THE MEDIUM (I.E. WHEN  IND.NE.3),
                THE HORIZONTAL AND VERTICAL COMPONENTS OF THE DIS-
                PLACEMENT ARE ALWAYS COMPUTED AS IF WITHIN AN
                INFINITE HALFSPACE (EVEN IF THE RECEIVER IS SITUATED AT
                AN INTERFACE).
                SIMILARLY, WHEN THE FIRST LAYER IS LIQUID (OCEAN),
                THE COEFFICIENTS OF CONVERSION ARE NOT USED.
      ITMAX...  NUMBER OF ITERATIONS PERMITTED IN A 2-POINT RAY
                TRACING. IF ITMAX=0, THEN ITMAX=20.
      NLAY...   CONTROLS THE GENERATION OF THE NUMERICAL CODE OF
                PRIMARY REFLECTED WAVES FROM THE LAST INTERFACE.
             NLAY=0... THE PRIMARY REFLECTED WAVES FROM THE
                 LAST INTERFACE ARE GENERATED (WHEN IBP.NE.0 OR
                 IBS.NE.0).
             NLAY=1... THEY ARE NOT GENERATED.
      LU1...    THE NUMBER OF THE FILE IN WHICH RAY DIAGRAMS,
                TRAVEL TIMES AND AMPLITUDES ARE STORED.
                FOR LU1=0, THESE QUANTITIES ARE NOT STORED.
      LU2.      THE NUMBER OF THE FILE IN WHICH DATA FOR THE
                CONSTRUCTION OF RAY SYNTHETIC SEISMOGRAMS ARE
                STORED.
                FOR LU2=0, THESE QUANTITIES ARE NOT STORED.
      LTAPE...  DUMMY PARAMETER
      KSH...    SWITCH CONTROLLING SELECTION OF THE TYPE OF COMPUTED
                WAVES.
             KSH=0... ONLY P AND SV WAVES ARE CONSIDERED. NO
                SH WAVES ARE COMPUTED.
             KSH=1... ONLY SH WAVES ARE CONSIDERED. NO P, SV
                WAVES ARE COMPUTED. ANY WAVE CONTAINING AN ELEMENT
                ALONG WHICH THE WAVE PROPAGATES AS A P WAVE IS
                AUTOMATICALLY EXCLUDED.
             KSH=2... P, SV AND SH WAVES ARE CONSIDERED TOGETHER.
                FOR ANY PURELY SHEAR WAVE SPECIFIED BY THE
                CODE, BOTH SV AND SH AMPLITUDES ARE COMPUTED.
      ILOC...   SWITCH SPECIFYING THE LOCATION OF RECEIVERS.
             ILOC=0... THE RECEIVERS ARE LOCATED AT THE EARTH'S
                SURFACE.
             ILOC=1... THE RECEIVERS ARE LOCATED ON A VERTICAL
                PROFILE.
             ILOC.GT.100... THE RECEIVERS ARE LOCATED ON THE
                (ILOC-100)TH INTERFACE.

10) ONE CARD OR GROUP OF CARDS: SPECIFICATION OF RECEIVER POSITIONS.

WHEN THE RECEIVERS ARE SITUATED ALONG THE EARTH'S SURFACE OR
AN INTERFACE (ILOC.NE.1), THEY ARE SPECIFIED BY THEIR X-COORDINATES.
WHEN THE RECEIVERS ARE SITUATED ALONG A VERTICAL PROFILE
(ILOC.EQ.1), THEY ARE SPECIFIED BY THEIR Z-COORDINATES.

WHEN MEP.LT.0: (DENOTE NDST=-MEP)
   DST(1),DST(2),...,DST(NDST),XVSP
                THE RECEIVERS ARE DISTRIBUTED IRREGULARLY ALONG
                THE EARTH'S SURFACE, AN INTERFACE, OR A VERTICAL
                PROFILE.
             DST(I) DENOTES THE VARYING COORDINATE OF THE I-TH
                RECEIVER. DST(I) SHOULD FORM A MONOTONOUSLY IN-
                CREASING (DECREASING) SERIES.
             XVSP... THE X-COORDINATE OF THE VERTICAL PROFILE
                FOR ILOC.EQ.1. OTHERWISE, DUMMY PARAMETER.

WHEN MEP.GT.0:
   RMIN,RSTEP,XVSP
                THE RECEIVERS ARE DISTRIBUTED REGULARLY ALONG
                THE EARTH'S SURFACE, AN INTERFACE, OR A VERTICAL
                PROFILE.
                THE VARYING COORDINATE OF THE I-TH RECEIVER IS
                GIVEN BY THE FORMULA DST(I)=RMIN+RSTEP*(I-1).
             XVSP... THE X-COORDINATE OF THE VERTICAL PROFILE
                FOR ILOC.EQ.1. OTHERWISE, DUMMY PARAMETER.





11) ONE CARD, POSITION OF THE SOURCE,ETC.
    XSOUR,ZSOUR,TSOUR,REPS,REPS1,BRD(1),BRD(2)
       XSOUR... THE X-COORDINATE OF THE SOURCE
       ZSOUR... THE Z-COORDINATE OF THE SOURCE
                THE SOURCE MAY BE SITUATED AT ANY POINT OF THE
                MODEL, INCLUDING ITS BOUNDARIES. THE SOURCE
                SITUATED ON AN INTERFACE (WITH EXCEPTION OF THE
                LAST INTERFACE) IS CONSIDERED AS A POINT OF THE
                LAYER LYING UNDER THE MENTIONED INTERFACE. IN CASE
                OF LAST INTERFACE, THE SOURCE IS CONSIDERED AS
                A POINT OF THE LAST LAYER. THE SOURCE CANNOT BE
                SITUATED IN THE LAYER OF ZERO THICKNESS (IN CASE
                OF COINCIDENCE OF NEIGHBOURHOOD INTERFACES).
                IF ABS(ZSOUR).LT.0.0001, THE SOURCE IS CONSIDERED
                AS LYING ON THE EARTH SURFACE AND ITS REAL
                Z COORDINATE IS DETERMINED AUTOMATICALLY (THIS IS
                VERY IMPORTANT IN CASE WHEN THE EARTH'S SURFACE IS
                CURVED AND ACTUAL Z-COORDINATE OF THE SOURCE IS NOT
                A PRIORI KNOWN).
      TSOUR...  INITIAL (HYPOCENTRAL) TIME
      REPS...   RADIUS OF THE VICINITY OF A RECEIVER. IF A RAY
                HAS ITS TERMINATION POINT WITHIN THIS VICINITY,
                THE ITERATIONS (2-POINT RAY TRACING) TERMINATE.
                IF REPS IS NOT SPECIFIED, THEN REPS=0.05.
      REPS1...  CONTROLS THE TERMINATION OF ITERATIONS TO BOUNDARY
                RAYS. IF THE DISTANCE BETWEEN TERMINATION POINTS
                OF TWO SUCCESSIVELY GENERATED RAYS IS LESS THAN
                REPS1, THE ITERATIONS TERMINATE. IF REPS1 IS NOT
                SPECIFIED, THEN REPS1=0.05*REPS.
      BRD(1),BRD(2)... X-COORDINATES OF THE LEFT AND RIGHT VER-
                TICAL BOUNDARIES OF THE PART OF THE MODEL - CALLED
                FURTHER ON SUBMODEL - IN WHICH THE COMPUTATIONS
                FOR THE GIVEN SOURCE WILL BE PERFORMED. THESE
                BOUNDARIES MAY GENERALLY DIFFER FROM VERTICAL
                BOUNDARIES OF THE WHOLE MODEL. THEY MUST BE, HOW-
                EVER, SITUATED WITHIN IT. IF ANY OF THE BOUNDARIES
                OF THE SUBMODEL IS SPECIFIED OUTSIDE THE WHOLE
                MODEL, THEN IT IS AUTOMATICALLY CORRECTED AND
                TAKEN AS THE CORRESPONDING BOUNDARY OF THE MODEL.
                IF BRD(1) AND BRD(2) ARE NOT SPECIFIED, THEN THEY
                ARE AUTOMATICALLY CONSIDERED AS VERTICAL BOUNDARIES
                OF THE WHOLE MODEL.

12) ONE CARD, QUANTITIES THAT CONTROL THE BASIC SYSTEM OF INITIAL
    ANGLES IN THE TWO-POINT RAY TRACING, ETC.
   DT,AMIN1,ASTEP1,AMAX1,AMIN2,ASTEP2,AMAX2,AC
      DT...     TIME STEP IN THE INTEGRATION OF THE RAY-TRACING
                SYSTEM. DT SHOULD BE GREATER THAN ZERO.
                IF DT.LT.0.00001, THEN DT=1.
      AMIN1,ASTEP1,AMAX1... DETERMINE THE SYSTEM OF INITIAL
                ANGLES FOR PRIMARY REFLECTED WAVES (AND POSSIBLY
                FOR OTHER MANUALLY GENERATED ELEMENTARY WAVES,THE
                FIRST ELEMENT OF WHICH STRIKES THE INTERFACE SITUATED
                BELOW THE SOURCE). SPECIFIED IN RADIANS. THE ANGLES
                ARE MEASURED IN THE PLANE (X,Z) FROM THE POSITIVE
                DIRECTION OF THE X-AXIS, POSITIVE CLOCKWISE.
      AMIN2,ASTEP2,AMAX2... DETERMINE THE SYSTEM OF INITIAL ANGLES
                FOR THE DIRECT WAVES (AND POSSIBLY FOR OTHER MA-
                NUALLY GENERATED ELEMENTARY WAVES,THE FIRST ELE-
                MENT OF WHICH STRIKES THE INTERFACE ABOVE THE
                SOURCE). SPECIFIED IN RADIANS. THE ANGLES ARE MEA-
                SURED IN THE PLANE (X,Z) FROM THE POSITIVE DIRECTION
                OF THE X-AXIS< POSITIVE CLOCKWISE.
      AC...     REQUIRED ACCURACY OF INTEGRATION OF THE RAY TRA-
                CING SYSTEM. RECOMENDED VALUES: 0.0001 - 0.001.
                IF AC NOT SPECIFIED, THEN AC=0.0001.

13) ARBITRARY NUMBER OF CARDS, ENDING BY A BLANK CARD.
    MANUAL GENERATION OF NUMERICAL CODES OF ELEMENTARY WAVES.
    THE CARDS ARE INCLUDED ONLY WHEN IREAD.NE.0, SEE CARD NO.9
    FOR IREAD. THE LAST CARD OF THE SYSTEM IS A BLANK CARD.
   KC,KREF,CODE(1),CODE(2),...,CODE(KREF)
      KC,KREF,CODE(1),CODE(2),...,CODE(KREF)   THE NUMERICAL CODE
                OF THE WAVE. THE WHOLE RAY IS DIVIDED INTO ELE-
                MENTS, EACH OF WHICH LIES BETWEEN TWO SUCCESSIVE
                POINTS AT WHICH THE RAY STRIKES AN INTERFACE.
                IF THE END POINTS OF AN ELEMENT LIE ON DIFFE-
                RENT INTERFACES, IT IS CALLED A SIMPLE ELEMENT,
                WHEN THE END POINTS LIE ON THE SAME INTERFACE,
                IT IS CALLED A COMPOUND ELEMENT. ANY COMPOUND ELE-
                MENT IS FORMALLY REGARDED AS TWO SIMPLE ELEMENTS.
                THE PART OF A RAY BETWEEN THE SOURCE AND THE FIRST
                POINT AT WHICH THE RAY STRIKES AN INTERFACE IS
                CONSIDERED AS ONE ELEMENT OF THE RAY. IF THE SOURCE
                IS SITUATED ON AN INTERFACE AND A RAY FROM THE
                SOURCE HITS THE SAME INTERFACE, THE ELEMENT OF
                THE RAY IS ALSO COMPOUND ELEMENT.
                WHEN THE RAY CROSSES A LAYER OF ZERO THICKNESS
                (COINCIDING INTERFACES), THE RAY ELEMENT WITHIN
                THE LAYER IS FORMALLY CONSIDERED IN THE SAME WAY
                AS IN THE CASE OF A LAYER OF A NON-ZERO THICKNESS.
      KC...  KC.NE.0... THE CODE CORRESPONDS TO A REFLECTED, OR
                MULTIPLY REFLECTED, POSSIBLY CONVERTED WAVE.
             KC=1... AFTER LEAVING THE SOURCE THE RAY SHOULD
                STRIKE THE INTERFACE BELOW THE SOURCE.
             KC=-1... AFTER LEAVING THE SOURCE THE RAY SHOULD
                STRIKE THE INTERFACE ABOVE THE SOURCE.
             KC=0... DIRECT WAVE (NO REFLECTIONS AND NO
                CONVERSIONS AT INTERFACES, ONLY TRANSMISSIONS).
      KREF...   THE NUMBER OF SIMPLE ELEMENTS OF THE RAY.
                IN CASE OF REFRACTED WAVE KREF=1.
      CODE(I)...ABS(CODE(I))... THE NUMBER OF THE LAYER IN WHICH
                THE I-TH SIMPLE ELEMENT OF THE RAY IS SITUATED.
             CODE(I).GT.0 FOR A P-WAVE ELEMENT,
             CODE(I).LT.0 FOR A S-WAVE ELEMENT.
                 COMPOUND ELEMENTS ARE REGARDED AS TWO SIMPLE
                 ELEMENTS AND ARE DENOTED BY TWO EQUAL NUMBERS.
                 NOTE: CODE(I) IS AN INTEGER.

******************************************************************

THE TERMINATION OF COMPUTATIONS.
********************************

THE SYSTEM OF CARDS NO. 9 - 13 CAN BE REPEATED ANY NUMBER OF
TIMES TO PERFORM VARIOUS COMPUTATIONS FOR THE SAME MODEL.
EACH SYSTEM OF CARDS NO.9-13 PRODUCES ONE FILE WITH FORTRAN
REFERENCE NUMBER LU1 (IF LU1.NE.0) AND ANOTHER WITH REFERENCE
NUMBER LU2 (IF LU2.NE.0). THIS POSSIBILITY CAN BE USED ONLY IN
OPERATING SYSTEMS ALLOWING TO SPECIFY THE NUMBER OF A FILE
UTSIDE THE SOURCE PROGRAM (E.G., IBM OPERATING SYSTEM).
TO TERMINATE THE COMPUTATIONS, PUT ICONT=0 IN THE CARD NO.9.
IF ICONT=-1, THE COMPUTATIONS CONTINUE, BUT FOR A NEW MODEL.
THE CARDS ARE THEN READ IN STARTING FROM THE CARD NO.1.

******************************************************************

MOST IMPORTANT OUTPUT INFORMATION.
**********************************

THE PARAMETER IND.
******************
THE PARAMETER IND SPECIFIES THE REASON FOR THE TERMINATION OF
THE COMPUTED RAY OR TERMINATION OF COMPUTATIONS AT ALL. PARAMETER
IND IS STORED IN COMMON/RAY/, SEE BELOW.

IND=1   TERMINATION OF THE RAY AT THE LEFT BORDER OF THE SUBMODEL.
IND=2   TERMINATION OF THE RAY AT THE RIGHT BORDER OF THE SUBMODEL.
IND=3   TERMINATION OF THE RAY AT THE UPPER BORDER OF THE MODEL
        (I.E., ON THE FREE SURFACE).
IND=4   TERMINATION OF THE RAY AT THE BOTTOM BORDER OF THE MODEL.
IND=5   IN THE RUNGE-KUTTA METHOD, THE BASIC TIME STEP DT IS
        HALVED MORE THAN TEN TIMES, WITHOUT REACHING THE REQUIRED
        ACCURACY AC, SEE CARD NO. 12. THERE IS PERHAPS SOMETHING
        WRONG WITH YOUR VELOCITY DISTRIBUTION IN THE MODEL.
IND=6   IN CARD NO. 12, DT=0  (NOT ALLOWED). THE COMPUTATION
        OF THE RAY DIAGRAM DOES NOT START.
IND=7   IN CARD NO. 12, DT.LT.0  (NOT ALLOWED). THE COMPUTATION
        OF THE RAY DIAGRAM DOES NOT START.
IND=8   THE RAY STRIKES A DIFFERENT INTERFACE THAN INDICATED IN
        THE NUMERICAL CODE OF THE WAVE.
IND=9   OVERCRITICAL INCIDENCE AT AN INTERFACE WHERE THE NUMERICAL
        CODE REQUIRES A TRANSMISSION.
IND=12  TRAVEL TIME ALONG THE RAY IS GREATER THAN TMAX. IN THE
        PROGRAM SEIS, TMAX=10000 AUTOMATICALLY.
IND=13  THE NUMBER OF POINTS ALONG THE RAY IS LARGER THAN 400.
IND=14  THE POSITION OF THE SOURCE DOES NOT CORRESPOND TO THE
        NUMERICAL CODE OF THE WAVE, I.E.  ABS(CODE(1)) IS DIFFERENT
        FROM THE NUMBER OF THE LAYER WITHIN WHICH THE SOURCE
        IS SITUATED. THE SAME INDICATION APPEARS IF THE SOURCE IS
        SITUATED IN THE LAYER OF ZERO THICKNESS. (FOR DETAILS,
        SEE CARD NO.11).
IND=15  IN THE CASE WHERE THE COMPUTATION OF THE AMPLITUDES
        IS REQUIRED (MDIM.NE.0, SEE CARD NO. 9 FOR MDIM), AND THE
        NUMERICAL CODE OF THE WAVE REQUIRES THE COMPUTATION OF A
        WAVE REFLECTED FROM THE BOTTOM BOUNDARY OF THE MODEL, THE
        COMPUTATION OF THE RAY IS TERMINATED, AS THE VELOCITIES
        AND THE DENSITY BELOW THIS BOUNDARY ARE NOT KNOWN.
IND=16  THE NUMERICAL CODE REQUIRES A COMPUTATION OF A WAVE
        REFLECTED FROM A FICTITIOUS INTERFACE (III=-2). THE
        AMPLITUDE OF SUCH A WAVE EQUALS ZERO, THEREFORE THE
        COMPUTATION OF THE RAY IS TERMINATED. THE SAME INDICATION
        APPEARS IF THE TERMINATION POINT OF THE RAY LIES ON
        A FICTITIOUS PART OF AN INTERFACE.
IND=17  THE NUMERICAL CODE REQUIRES A REFLECTION FROM THE SECOND
        INTERFACE OF TWO COINCIDING INTERFACES. THE COMPUTATION
        OF THE RAY IS TERMINATED. THE SAME APPLIES TO OTHER IN-
        TERFACES, WHEN THE NUMBER OF COINCIDING INTERFACES IS
        LARGER THAN TWO.
IND=18  THE NUMERICAL CODE PRESCRIBES CONVERSION OF WAVE TRANS-
        MITTED ON AN INTERFACE WHICH COINCIDES WITH ANOTHER INTER-
        FACE. THE COMPUTATION OF THE RAY IS TERMINATED SINCE
        KINEMATICALLY AND DYNAMICALLY IDENTICAL RAY CAN BE OBTAI-
        NED WITH NUMERICAL CODE WHICH DOES NOT INCLUDE THE CORRES-
        PONDING CONVERSION.
IND=19  THE ROUTINE MODEL USED FOR THE COMPUTATIONS DOES NOT COR-
        RESPOND TO THAT SPECIFIED IN THE INPUT DATA, SEE PARAMETER
        MM IN INPUT DATA CARD NO.1. THE COMPUTATION DOES NOT START.
IND=20  TWO NEIGHBOURHOOD INTERFACES IN THE MODEL INTERSECT EACH
        OTHER. THE COMPUTATION TERMINATES.
IND=21  A RECTANGULAR NETWORK OF VELOCITY VALUES FOR ONE LAYER
        DOES NOT COVER THE WHOLE LAYER. THE COMPUTATION TERMINATES.
IND=22  TERMINATION OF THE RAY AT THE VERTICAL PROFILE.
IND=50  THE SOURCE IS SITUATED OUTSIDE THE MODEL, THE COMPUTATION
        OF THE RAY DIAGRAM DOES NOT START.
IND.GT.100  IF THE NUMERICAL CODE OF THE WAVE REQUIRES THE TER-
        MINATION OF THE COMPUTATION AT SOME INNER INTERFACE AND
        THE COMPUTATION OF THE RAY IS SUCCESSFULLY FINISHED, THEN
        (IND-100) GIVES THE NUMBER OF THE INTERFACE AT WHICH THE
        COMPUTATION OF THE RAY IS TERMINATED.

THE HISTORY OF THE RAY.
***********************

ALL THE IMPORTANT INFORMATION ABOUT ANY COMPUTED RAY IS STORED IN
COMMON/RAY/AY(12,400),DS(11,50) KINT(50),MREG,N,IREF,IND,IND1.
THIS INFORMATION IS LOST AS SOON AS A NEW RAY BEGINS TO BE COM-
PUTED. A MAXIMUM NUMBER OF 400 POINTS ALONG THE RAY IS ALLOWED.

THE MEANING OF THE ABOVE QUANTITIES IS THE FOLLOWING:

   N...      NUMBER OF POINTS ALONG THE RAY.
   IREF...   NUMBER OF REFLECTIONS/TRANSMISSIONS, INCLUDING
             REFLECTIONS/TRANSMISSIONS AT FICTITIOUS INTERFACES.
   IND...    EXPLAINED ABOVE.
   IND1...   IABS(IND1)... THE NUMBER OF THE LAST INTERFACE
             HIT BY THE RAY BEFORE THE COMPUTATION OF THE
             RAY WAS TERMINATED (OR THE INTERFACE AT WHICH
             THE COMPUTATION WAS TERMINATED).
         IND1=0... THE COMPUTATION OF THE RAY TERMINATED
             ON ONE OF THE VERTICAL BOUNDARIES OF THE MODEL.
             THE RAY DID NOT CROSS ANY INTERFACE ON ITS WAY
             FROM THE SOURCE TO THE TERMINATION POINT.
         IND1.LT.0...  THERE IS AT LEAST ONE COMPOUND ELEMENT
             ALONG THE RAY.
   MREG...   CONTROLS THE COMPUTATION OF AMPLITUDES AT THE TER-
             MINATION POINT OF THE RAY,SEE INPUT CARD NO. 9.

THE QUANTITIES AY(J,I) GIVE SOME INFORMATION ABOUT THE I-TH COMPUTED
POINT OF THE RAY. I=1 CORRESPONDS TO THE SOURCE, I=N TO
THE TERMINATION POINT.
   AY(1,I)...THE TIME CORRESPONDING TO THE POINT.
   AY(2,I)...X-COORDINATE OF THE POINT.
   AY(3,I)...Z-COORDINATE OF THE POINT.
   AY(4,I)...X COMPONENT OF THE SLOWNESS VECTOR AT THE POINT.
   AY(5,I)...Z COMPONENT OF THE SLOWNESS VECTOR AT THE POINT.
   AY(6,I)...VELOCITY V AT THE POINT.
   AY(7,I)...DV/DX
   AY(8,I)...DV/DZ.
   AY(9,1)...SECOND DERIVATIVE OF V WITH RESPECT TO X.
   AY(10,I)..SECOND DERIVATIVE OF V WITH RESPECT TO X AND Z.
   AY(11,I)..SECOND DERIVATIVE OF V WITH RESPECT TO Z.
   AY(12,I)..QUALITY FACTOR AT THE POINT.

THE QUANTITIES DS(J,K) AND KINT(K) GIVE SOME INFORMATION ABOUT THE
K-TH POINT OF REFLECTION/TRANSMISSION. THE CASE K=1 CORRESPONDS
TO THE INTERFACE WHICH WAS HIT BY THE RAY FIRST.

   DS(1,K)...THE RADIUS OF THE CURVATURE OF THE INTERFACE.
   DS(2,K)...X COMPONENT OF THE UNIT VECTOR TANGENT TO THE RAY
             AT THE POINT OF INCIDENCE, IN LOCAL COORDINATE
             SYSTEM CONNECTED WITH THE POINT OF INCIDENCE.
             Z-AXIS OF THIS SYSTEM IS PERPENDICULAR TO THE
             INTERFACE AND POINTS INTO THE MEDIUM IN WHICH
             INCIDENT WAVE PROPAGATES. Y-AXIS IS PERPENDICULAR
             TO THE PLANE OF THE MODEL, ITS POSITIVE DIRECTION
             POINTS BEHIND THE (X,Z) PLANE. X-AXIS IS TANGENT
             TO THE INTERFACE, ITS POSITIVE DIRECTION IS
             CHOSEN SO THAT THE LOCAL COORDINATE SYSTEM IS
             RIGHT-HANDED.
             IN CASE OF TERMINATION POINT OF THE RAY ON ONE
             OF THE VERTICAL BOUNDARIES OF THE MODEL,
             DS(2,K)=AY(4,N).
   DS(3,K)...Z COMPONENT OF THE UNIT VECTOR TANGENT TO THE RAY
             AT THE POINT OF INCIDENCE, IN LOCAL COORDINATE
             SYSTEM CONNECTED WITH THE POINT OF INCIDENCE
             (FOR DETAILS SEE ABOVE).
             IN CASE OF TERMINATION POINT OF THE RAY ON ONE
             OF THE VERTICAL BOUNDARIES OF THE MODEL,
             DS(3,K)=AY(5,N).
   DS(4,K)...VP AT THE POINT OF INCIDENCE.
   DS(5,K)...VS AT THE POINT OF INCIDENCE.
   DS(6,K)...VP AT OTHER SIDE OF INTERFACE.
   DS(7,K)...VS AT OTHER SIDE OF INTERFACE.
   DS(8,K)...DENSITY AT THE POINT OF INCIDENCE.
   DS(9,K)...DENSITY AT THE OTHER SIDE OF INTERFACE.
   DS(10,K)..THE SAME AS FOR DS(2,K) BUT FOR GENERATED WAVE.
   DS(11,K)..THE SAME AS FOR DS(3,K) BUT FOR GENERATED WAVE.

   KINT(K)...KINT(K).GT.0... GIVES THE SEQUENTIAL NUMBER OF
             THE POINT OF THE RAY WHICH CORRESPONDS TO THE
             K-TH POINT OF REFLECTION/TRANSMISSION.
         KINT(K)=0...    THE K-TH POINT OF REFLECTION/
             TRANSMISSION HAS A FORMAL MEANING, AS IT CORRESPONDS
             TO A COMPOUND  ELEMENT OF THE RAY, WITH BOTH END
             POINTS OF THE ELEMENT AT THE SAME INTERFACE.
         KINT(K)=-1...   THE K-TH POINT OF REFLECTION/
             TRANSMISSION CORRESPONDS TO A POINT ON THE INTERFACE
             WHICH COINCIDES WITH THE INTERFACE SITUATED ABOVE IT.

DESCRIPTION OF INTERFACES.
**************************

THE INFORMATION ABOUT THE INTERFACES IS STORED IN
COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20),
BRD(2),III(30,20),NPNT(20),NINT.

   NINT...   SPECIFIED IN INPUT DATA CARD NO.2.
   NPNT(I)...SPECIFIED IN INPUT DATA CARD NO.2.
   BRD(1)... X-COORDINATE OF THE LEFT VERTICAL BOUNDARY
             OF THE SUBMODEL. SPECIFIED IN INPUT CARD NO.11.
   BRD(2)... X-COORDINATE OF THE RIGHT VERTICAL BOUNDARY
             OF THE SUBMODEL. SPECIFIED IN INPUT CARD NO.11.

BETWEEN TWO SUCCESSIVE POINTS OF THE INTERFACE, WITH THE X-COORDINATE
OF THE LEFT-HAND POINT DENOTED BY X1, THE CUBIC PARABOLA
APPROXIMATING THE INTERFACE HAS A FORM:
Z = A+B*(X-X1)+C*(X-X1)**2+D*(X-X1)**3

   X1(K,I)...THE X-COORDINATE OF THE K-TH POINT OF THE I-TH
             INTERFACE.
   A1(K,I)...THE Z-COORDINATE OF THE K-TH POINT OF THE I-TH
             INTERFACE (CORRESPONDS TO A IN THE ABOVE FORMULA).
   B1(K,I),C1(K,I),D1(K,I)...  THE COEFFICIENTS OF THE CUBIC
             PARABOLA IN THE K-TH ELEMENT OF THE I-TH INTERFACE.
             THEY CORRESPOND TO B,C,D IN THE ABOVE FORMULA. THE
             K-TH ELEMENT IS BETWEEN THE K-TH AND (K+1)TH POINT.
   III(K,I)..HAS THE SAME MEANING AS III READ FROM THE INPUT
             CARD NO.3.

DESCRIPTION OF VELOCITY, DENSITY AND QUALITY FACTOR DISTRIBUTIONS.
******************************************************************

INFORMATION ABOUT THE VELOCITY, DENSITY AND QUALITY FACTOR
DISTRIBUTION IS STORED IN THE FOLLOWING COMMONS:
COMMON/MEDIM/V(1000),SX(350),SY(350),NX(20),NY(20),NVS(19),
   PTOS(19)  FOR MM=0,2, SEE INPUT CARD NO.1,
COMMON/MEDIM/V1(19),V2(19),NVS(19),PTOS(19) FOR MM=1,
COMMON/VCOEF/A02(1000),A20(1000),A22(1000)  FOR MM=0,2,
COMMON/DEN/RHO1(19),RHO2(19),NRHO.
COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19),
   NQP(19),NQS(19),NABS,
COMMON/SOUR/ROS,VPS,VSS.

   V...      VELOCITY VALUES SPECIFIED IN INPUT DATA CARD NO.4D
             FOR MM=0,2.
   SX,SY...  COORDINATES OF VERTICAL AND HORIZONTAL LINES IN THE
             VELOCITY NETWORK, SPECIFIED IN INPUT DATA CARDS
             NO.4B AND 4C FOR MM=0,2.
   NX,NY...  NUMBER OF VERTICAL AND HORIZONTAL LINES IN THE
             VELOCITY NETWORK, SPECIFIED IN INPUT DATA CARDS
             NO.4A (SEE MX,MY) FOR MM=0,2.
   NVS...    SEE A DETAILED DESCRIPTION IN DATA CARD NO.5.
   PTOS...   RATIOS OF P AND S VELOCITIES IN INDIVIDUAL LAYERS,
             SEE DETAILS IN CARD NO. 7.
   V1,V2...  VELOCITY VALUES SPECIFIED IN INPUT DATA CARDS NO.
             4A AND 4B FOR MM=1.
   A02,A20,A22.. BASIC COEFFICIENTS OF THE BICUBIC SPLINE OR BI-
             LINEAR INTERPOLATION AT GRID POINTS SPECIFIED IN THE
             SAME WAY AS FOR V (SEE DESCRIPTION UNDER THE INPUT
             DATA CARD NO.4D FOR MM=0,2). A02 CORRESPONDS TO THE
             FIRST (SECOND) DERIVATIVE OF VELOCITY WITH RESPECT
             TO Z, A20 WITH RESPECT TO X AND A22 TO THE SECOND
             (FOURTH) DERIVATIVE WITH RESPECT TO X AND Z IN CASE
             OF BILINEAR (BICUBIC SPLINE) INTERPOLATION.
   RHO1,RHO2... QUANTITIES THAT SPECIFY THE RELATION BETWEEN
             THE DENSITIES AND P WAVE VELOCITIES IN INDIVIDUAL
             LAYERS, DENSITY=RHO1(I)+RHO2(I)*VP IN THE I-TH LAYER,
             SEE CARD NO.6A. FOR NRO=0 (SEE CARD NO.5 FOR NRO)
             RHO1(I)=1.7, RHO2(I)=0.2 IN ALL LAYERS.
             FOR PTOS(1).GE.100 (SEE CARD NO. 7 FOR PTOS(1)),
             RHO1(1)=1,RHO2(1)=0 (A LIQUID LAYER).
   NRHO... NRHO=0, THE UPPERMOST LAYER IS SOLID, FOR
           NRHO=1, THE UPPERMOST LAYER IS LIQUID. IN THE
             LATTER CASE, THE S VELOCITY=0 AND THE DENSITY=1
             IN THE LAYER. ALSO MREG=1 IN THIS CASE.
   QP1,QP2,QP3,QS1,QS2,QS3... QUANTITIES WHICH SPECIFY THE RE-
             LATION BETWEEN THE QUALITY FACTOR AND THE CORRES-
             PONDING VELOCITY, SEE INPUT DATA CARD NO.6B.
   NQP,NQS...SWITCHES WHICH DETERMINE THE TYPE OF THE RELATION
             BETWEEN THE QUALITY FACTOR AND THE CORRESPONDING
             VELOCITY, SEE INPUT DATA CARD NO.6B.
   ROS,VPS,VSS... VALUES OF DENSITY, P-WAVE AND S-WAVE VELOCITY
             AT THE SOURCE.
NOTE: FOR MM=0,2 THE INFORMATION ABOUT VELOCITY DISTRIBUTION
IN ONE LAYER IS STORED IN THE ABOVE VARIABLES IN THE SAME
WAY AS DESCRIBED IN INPUT DATA CARDS NO.4. THE DATA FOR INDI-
VIDUAL LAYERS ARE STORED SUCCESSIVELY FROM THE TOP TO THE
BOTTOM, ONE AFTER ANOTHER.

DESCRIPTION OF CERTAIN AUXILIARY QUANTITIES
*******************************************

AUXILIARY QUANTITIES WHICH ARE USED IN VARIOUS PARTS OF THE PROGRAM
ARE STORED IN COMMON /AUXI/ AND COMMON /AUXX/.

COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT,NTR

   INTR...   NUMBER OF THE INTERFACE JUST STRUCK BY THE RAY.
   INT1...   NUMBER OF THE INTERFACE WHICH WAS HIT BEFORE INTR.
   IOUT...   INDEX USED WHEN THE POINT OF INCIDENCE OF A RAY
             AT AN INTERFACE IS DETERMINED IN AN ITERATIVE WAY
             (NOT USED IN THIS VERSION).
   KRE...    INDEX SPECIFYING WHAT TYPE OF WAVE IS CONSIDERED.
         KRE=0... REFRACTED WAVE IS CONSIDERED.
         KRE=KREF... (SEE CARD NO.14 FOR KREF), MULTIPLY
             REFLECTED WAVE IS CONSIDERED.
   IREFR...  INDEX SPECIFYING IF THERE IS AT LEAST ONE COMPOUND
             ELEMENT OF THE RAY (IREFR=1) OR NOT (IREFR=0).
   LAY...    NUMBER OF LAYER JUST PASSED THROUGH THE RAY.
   ITYPE...  INDEX SPECIFYING THE TYPE OF THE WAVE ALONG THE
             INVESTIGATED ELEMENT OF THE RAY. ITYPE=1 FOR P
             WAVE, ITYPE=-1 FOR S WAVE.
   NDER...   INDEX SPECIFYING WHETHER THE SECOND DERIVATIVES OF
             VELOCITY SHOULD BE EVALUATED (NDER=1  FOR MDIM.GT.0)
             OR NOT (NDER=0  FOR MDIM=0).
   IPRINT... BLANK PARAMETER IN THIS VERSION.
   MPRINT... CONTROLS THE PRINTING OF THE DESCRIPTION OF THE
             MODEL ON THE LINE PRINTER, SEE INPUT DATA CARD NO.1.
   NTR...    SERVES FOR THE PURPOSE OF A MORE DETAILED DETER-
             MINATION OF THE REASON OF THE TERMINATION OF THE RAY.
             IT SHOWS THE PLACE IN THE ROUTINE OUTP WHERE THE
             COMPUTATION OF THE RAY WAS TERMINATED (SEE THE LISTING
             OF THE ROUTINE OUTP).

COMMON/AUXX/MMX(20),MMY(20),MMXY(20),MAUX

USED ONLY IN THE ROUTINE MODEL AND RELATED ROUTINES IN CASE
THAT MM=0 OR 2.

   MMX(I)... SPECIFIES THE NUMBER OF THE FIRST COMPONENT OF THE
             VECTOR SX (SEE /MEDIM/) CORRESPONDING TO THE I-TH
             LAYER.
   MMY(I)... SPECIFIES THE NUMBER OF THE FIRST COMPONENT OF THE
             VECTOR SY (SEE /MEDIM/) CORRESPONDING TO THE I-TH
             LAYER.
   MMXY(I)...SPECIFIES THE NUMBER OF THE FIRST COMPONENTS OF
             THE VECTORS V (SEE /MEDIM/) AND AO2,A20,A22 (SEE
             /VCOEF/) CORRESPONDING TO I-TH LAYER.
   MAUX...   INDEX SPECIFYING IF IT IS NECESSARY TO EVALUATE
             THE COEFFICIENTS OF BICUBIC SPLINE OR BILINEAR
             INTERPOLATION,OR WHETHER IT IS POSSIBLE TO USE THE
             COEFFICIENT DETERMINED IN THE PREVIOUS STEP.

******************************************************************

OUTPUT TABLES.
**************

ALL THE INPUT DATA ARE REPRODUCED ON THE LINE PRINTER. THE
OTHER OUTPUT IS CONTROLLED BY THE PARAMETERS MPRINT AND
MOUT, SEE INPUT CARDS NO.1 AND 9.

THE INDEX MPRINT:
*****************
THE INDEX MPRINT CONTROLS THE PRINTING OF THE DESCRIPTION OF THE
MODEL ON THE LINE PRINTER.
FOR MPRINT=0:
*************
INPUT DATA FOR THE MODEL ARE REPRODUCED.
FOR MPRINT=1:
*************
PRINTER PLOT OF THE VELOCITY DISTRIBUTION WITHIN THE PART OF THE
MODEL BOUNDED BY HORIZONTAL BOUNDARIES BMIN AND BMAX AND BY VERTICAL
BOUNDARIES BLEFT AND BRIGHT, SEE DETAILS IN THE DESCRIPTION OF
THE INPUT CARD NO.8. THE PRINTERPLOT SHOWS APPROXIMATELY THE
POSITION OF INTERFACES AND/OR THE POSITION OF VELOCITY ISOLINES
IN THE MODEL.
FOR MPRINT=2:
*************
THE SAME PRINTOUT AS IN THE CASE OF MPRINT=1, BUT WITH A DETAILED
DESCRIPTION OF INDIVIDUAL INTERFACES. FOR EACH INTERFACE, THE
FOLLOWING QUANTITIES ARE PRINTED:
 A) COEFFICIENTS OF CUBIC PARABOLAS CORRESPONDING TO INDIVIDUAL
ELEMENTS OF THE INTERFACE (FROM THE LEFT TO THE RIGHT VERTICAL
BOUNDARY OF THE WHOLE MODEL).
 B) COORDINATES OF 26 POINTS ALONG EACH INTERFACE, WITH REGULAR
STEP IN THE X-COORDINATE, FROM THE LEFT BOUNDARY 'BLEFT' TO THE
RIGHT BOUNDARY 'BRIGHT', SEE INPUT DATA CARD NO.8.

THE INDEX MOUT:
***************

THE INDEX MOUT CONTROLS THE PRINTING OF A DESCRIPTION OF THE WHOLE
PROCESS OF ITERATIONS (TWO-POINT RAY TRACING) AND OF ITS RESULTS
ON A LINE PRINTER.
FOR MOUT=0:
***********
NO PRINT OF RESULTS. ONLY THE LISTING OF ELEMENTARY WAVES, FOR
WHICH THE COMPUTATIONS ARE PERFORMED, IS GIVEN. THE "EXTERNAL
CODE" (FULL CODE) CORRESPONDS TO THE NUMERICAL CODE OF THE WAVE
DESCRIBED IN INPUT DATA CARD NO.14. THE "INTERNAL CODE" IS THE
SUCCESSIVE NUMBER OF THE ELEMENTARY WAVES GENERATED BY THE ROUTINE
GENER. FOR MORE DETAILS, SEE DESCRIPTION OF THE ROUTINE GENER.
FOR MOUT=1 OR MOUT=2:
*********************
IN ADDITION TO THE NUMERICAL CODES OF INDIVIDUAL ELEMENTARY WAVES,
CERTAIN RESULTS OF ITERATIONS ARE PRINTED AT THE LINE PRINTER.
THE PRINTING IS DIFFERENT FOR MOUT=1 (ONLY BASIC PRINT) AND
MOUT=2 (HEAVY PRINTING).
THE FOLLOWING QUANTITIES ARE PRINTED IN VARIOUS SITUATIONS:
   IND, IND1... SEE EXPLANATION IN THE TABLE "THE PARAMETER
                IND" AND IN "THE HISTORY OF RAY".
   ITER...      THE SUCCESSIVE NUMBER OF THE ITERATION. IN THE
                PRINT OF THE FINAL RESULT OF ITERATIONS, IT
                GIVES THE FULL NUMBER OF ITERATIONS.
   II...        THE SUCCESSIVE NUMBER OF THE RECEIVER POSITION
                CORRESPONDING TO THE TERMINATION POINT OF THE
                RAY. THE RECEIVERS ARE NUMBERED FROM THE LEFT
                TO THE RIGHT.
   INDEX...     THE NUMBER OF SUCCESSFULLY ITERATED RAYS FOR A
                GIVEN ELEMENTARY WAVE.
   IRES...   IRES=1... BOUNDARY RAY TERTMINATING ON THE RECEIVER
                PROFILE WAS FOUND.
             IRES=0... BOUNDARY RAY TERMINATING ON THE RECEIVER
                PROFILE WAS NOT FOUND.
   KMAH...      KMAH INDEX, SEE DETAILS IN THE DESCRIPTION OF
                THE ROUTINE JPAR.
   X,Z...       COORDINATES OF THE TERMINATION POINT OF THE RAY.
   T...         TRAVEL TIME AT THE TERMINATION POINT OF THE RAY.
   TAST...      GLOBAL ABSORPTION FACTOR AT THE TERMINATION POINT
                OF THE RAY.
   AA...        PARAMETER OF THE RAY (INITIAL ANGLE) IN THE BASIC
                SYSTEM OF INITIAL ANGLES IN THE RAY SHOOTING, SEE
                DATA CARD NO.12.
   PNEW...      THE PARAMETER OF THE RAY (INITIAL ANGLE) IN ITERATIONS.
   DD...        COORDINATE OF THE RECEIVER (X-COORDINATE FOR RECEIVERS
                SITUATED ALONG THE EARTH'S SURFACE OR AN INTERFACE,
                Z-COORDINATE FOR RECEIVERS SITUATED ALONG A VERTICAL
                PROFILE).
   AX,AY,AZ...  AMPLITUDES OF RADIAL (HORIZONTAL ALONG THE PROFILE,
                POSITIVE TO THE RIGHT), TRANSVERSE (HORIZONTAL,
                PERPENDICULAR TO THE PROFILE) AND VERTICAL (POSITIVE
                UPWARDS) COMPONENTS OF DISPLACEMENT VECTOR AT THE
                TERMINATION POINT OF RAY (AX,AY,AZ SPECIFIED IN
                A RIGHT-HANDED COORDINATE SYSTEM).
   PHX,PHY,PHZ... PHASE SHIFTS CORRESPONDING TO THE ABOVE COMPONENTS
                AT THE TERMINATION POINT OF THE RAY.
   SPR...       THE GEOMETRICAL SPREADING AT THE TERMINATION
                POINT OF THE RAY.

ALL TABLES ARE PRINTED WITHOUT HEADINGS, BUT CERTAIN LINES ARE
SUPPLEMENTED BY THE ADDITIONAL TEXT: "ITERATIONS", "SUCCESSFUL
ITERATION", "BOUNDARY RAY", "CRITICAL RAY".
   GENERALLY, FOR EACH RAY ONE LINE IS PRINTED. THE EXCEPTION IS
THE SUCCESSFULLY ITERATED RAY, WHICH MAY BE SUPPLEMENTED BY SOME
ADDITIONAL INFORMATION ON AMPLITUDES (WHEN MDIM.GT.0,SEE LATER).
IN THE FOLLOWING, THE SYMBOL WW MEANS EITHER X (WHEN RECEIVERS
ARE SITUATED ALONG THE EARTH'S SURFACE OR ALONG AN INTERFACE),
OR Z (WHEN RECEIVERS ARE SITUATED ALONG A VERTICAL PROFILE).

IN THE BASIC SYSTEM OF SHOOTING, THE FOLLOWING QUANTITIES
ARE PRINTED WHEN MOUT=2:
             IND,IND1,W,T,AA.
      IN ITERATIONS, FOR MOUT=2:
             "ITERATIONS", IND,IND1,ITER,DD,W,T,PNEW.
      IN SUCCESSFUL ITERATIONS, WHEN MOUT=2:
             "SUCCESSFUL ITERATION",IND,IND1,ITER,II,INDEX,DD,W,T,PNEW,
             AX,PHX,AY,PHY,AZ,PHZ,SPR,TAST,KMAH.
      IN SUCCESSFUL ITERATIONS, WHEN MOUT=1 AND MDIM=0:
              DD,X,Z,T,PNEW,IND,IND1,ITER,II.
      IN SUCCESSFUL ITERATIONS, WHEN MOUT=1 AND MDIM.NE.0:
              DD,X,Z,T,TAST,PNEW,SPR,AX,PHX,AY,PHY,AZ,PHZ,IND,IND1,ITER,
              II,KMAH.
      IN CASE OF BOUNDARY RAYS, ALL ITERATIONS ARE PRINTED WHEN MOUT=2:
              "BOUNDARY RAY",IND,IND1,ITER,W,T,PNEW.
      IN CASE OF BOUNDARY RAYS, ONLY THE FINAL ITERATION IS PRINTED WHEN
      MOUT=1:
              X,Z,T,PNEW,IND,IND1,IRES.
      IN CASE OF CRITICAL RAYS, ALL ITERATIONS ARE PRINTED WHEN MOUT=2:
              "CRITICAL RAY", IND,IND1,ITER,W,T,PNEW.
      IN CASE OF CRITICAL RAYS, ONLY THE FINAL ITERATION IS PRINTED
      WHEN MOUT=1:
              X,Z,T,PNEW,IND,IND1,IRES.


OUTPUT TO THE FILE LU1.
***********************

IN THE CASE THAT LU1.NE.0, CERTAIN INFORMATION ABOUT THE MODEL,
THE MOST IMPORTANT RESULTS OF COMPUTATION ARE STORED IN THE
FILE LU1 FOR POSSIBLE FURTHER PROCESSING IN OTHER PROGRAMS.
THE PROGRAM RAYPLOT, INCLUDED IN THIS PACKAGE, MAY BE USED FOR
PLOTTING OF DATA STORED IN THE FILE LU1.

WHEN LU1=0, NO DATA ARE STORED.

THE DATA CONSIST OF THE FOLLOWING PARTS:

A) INFORMATION ABOUT POSSIBLE TERMINATION OF COMPUTATION, ABOUT
   THE SOURCE AND ABOUT INTERFACES:
1) ICONT.
   ICONT=0... INDICATION OF THE END OF DATA STORED IN THE FILE
   LU1. THIS IS THE LAST QUANTITY STORED IN THE FILE LU1.
   ICONT=1... INDICATION THAT A NEW DATA SET (2-8) FOR A NEW RAY
   DIAGRAM MAY FOLLOW.
2) NINT, (NPNT(I),I=1,NINT).
   FOR THE MEANING OF THESE SYMBOLS SEE INPUT DATA CARD NO.2.
3) FOR I=1,2,...,NINT:
   (A1(J,I),B1(J,I),C1(J,I),D1(J,I),X1(J,I),III(J,I),J=1,NC),
   WHERE NC=NPNT(I)-1.
   FOR THE MEANING OF THESE SYMBOLS SEE INPUT DATA CARD NO.3
   AND THE DESCRIPTION OF INTERFACES.

B) INFORMATION ABOUT THE SOURCE POSITION AND PARAMETERS OF THE
   MEDIUM AT THE SOURCE:
4) X0,Z0,ROS,VPS,VSS.
   COORDINATES OF THE SOURCE, THE DENSITY AND P- AND S-WAVE VELOCITY
   AT THE SOURCE.

C) INFORMATION ABOUT RAYS:
   FOR EACH SUCCESSFULLY ITERATED RAY WITH THE TERMINATION
   POINT AT SOME RECEIVER AT THE EARTH'S SURFACE,THE QUAN-
   TITIES SPECIFIED SUB 5 AND 6 ARE STORED. THE DATA SPECIFIED
   IN 5 AND 6 ARE STORED SUCCESSIVELY FOR ALL RAYS WITHIN
   ONE RAY DIAGRAM, AN ARBITRARY NUMBER OF TIMES. AFTER THE
   LAST RAY, THE PARAMETER N IS SET TO ZERO IN THE MAIN
   PROGRAM. THEN THE QUANTITIES DESCRIBED SUB D ARE STORED.
5) N,IND
   IF(N.NE.0)... NUMBER OF POINTS ALONG THE RAY.
   IF N=0... INDICATION OF THE END OF THE RAY DIAGRAM.
   IND... THE INDEX SPECIFYING THE REASON OF THE TERMINATION
   OF THE RAY. SEE DETAILED EXPLANATION ABOVE. WHEN THE RECEIVERS
   ARE SITUATED ALONG THE EARTH'S SURFACE, THE RAYS WITH IND=3
   ARE STORED IN LU1. WHEN THE RECEIVERS ARE SITUATED ALONG THE
   (IND-100)TH INTERFACE, THE RAYS WITH IND.GT.100 ARE STORED IN
   LU1. WHEN THE RECEIVERS ARE SITUATED ALONG A VERTICAL PROFILE,
   THE RAYS WITH IND=22 ARE STORED IN LU1.
6) (AY(2,I),AY(3,I),I=1,N).
   FOR THE MEANING OF THESE SYMBOLS SEE THE DESCRIPTION OF
   COMMON/RAY/IN THE HISTORY OF THE RAY.

D) INFORMATION ABOUT TRAVEL-TIMES AND AMPLITUDES
7) NS.
   NUMBER OF RAYS WHICH TERMINATE AT THE SPECIFIED SYSTEM OF
   RECEIVERS.
8) (INDI(I),XCOOR(I),TIME(I),TAS(I),ANG(I),AMPX(I),AMPY(I),
   AMPZ(I),PHAX(I),PHAY(I),PHAZ(I),I=1,NS).
   THE MEANING OF THE ABOVE SYMBOLS IS AS FOLLOWS:
   A SPECIAL PARAMETER IND WHICH SPECIFIES THE REASON OF THE TERMINATION
   OF THE COMPUTATION OF THE RAY, THE X-COORDINATE OF THE TERMINATION
   POINT OF THE RAY, THE CORRESPONDING TRAVEL TIME, GLOBAL ABSORPTION
   FACTOR T*, PARAMETER OF THE RAY (INITIAL ANGLE), AMPLITUDES OF RADIAL,
   TRANSVERSE AND VERTICAL COMPONENTS AND CORRESPONDING PHASE SHIFTS.
   ALL THESE DATA FOLLOW SUCCESSIVELY FOR ALL THE RAYS WHICH REACHED
   THE RECEIVER PROFILE (I=1,2,....,NS).

OUTPUT TO THE FILE LU2.
***********************

IN THE FILE LU2, THE DATA FOR COMPUTATION OF RAY SYNTHETIC
SEISMOGRAMS FOR INDIVIDUAL RECEIVER POSITIONS ARE STORED TOGETHER
WITH SOME OTHER RELEVANT INFORMATION. ALL THE DATA CORRESPOND
TO ONE SYSTEM OF CARDS NO.9-13.
TO PERFORM THE COMPUTATION OF SYNTHETIC SEISMOGRAMS AND THEIR
PLOTTING, THE PROGRAMS SYNTPL AND SEISPLOT INCLUDED IN THIS PACKAGE
MAY BE USED.

THE DATA ARE STORED IN THE FOLLOWING SUCCESSION:

1) MTEXT                                 FORMAT(17A4)
   ARBITRARY ALPHANUMERIC TEXT DESCRIBING THE MODEL.
2) NDST                                  FORMAT(26I3)
   NDST IS A NUMBER OF RECEIVER POSITIONS.
3) XSOUR,ZSOUR,TSOUR,RSTEP,ROS,VPS,VSS   FORMAT(8F10.5)
   XSOUR,ZSOUR... X-COORDINATE AND Z-COORDINATE OF THE SOURCE
   TSOUR...  INITIAL (HYPOCENTRAL) TIME.
   RSTEP...  (AVERAGE) DIFFERENCE BETWEEN THE X-COORDINATES
             OF TWO NEIGHBOURHOOD RECEIVER POSITIONS.
   ROS...    THE DENSITY AT THE SOURCE.
   VPS,VSS...THE P- AND S-WAVE VELOCITY AT THE SOURCE.
4) DST(I),I=1,NDST                        FORMAT(8F10.5)
   COORDINATES OF RECEIVER POSITIONS (X-COORDINATES FOR RECEIVERS
SITUATED ALONG THE EARTH'S SURFACE OR AN INTERFACE, Z-COORDINATES
FOR RECEIVERS SITUATED ALONG A VERTICAL PROFILE).
5) NCD,II,T,AX,AY,AZ,PHX,PHY,PHZ,PNEW,TAST
                                FORMAT(2I3,F10.5,2E12.6,4F10.6)
   NCD...    NCODE=IABS(NCD), NCODE IS THE INTERNAL CODE OF THE
             WAVE (NCODE=1,2,...), SEE DETAILED DESCRIPTION OF
             THE INTERNAL CODE IN THE ROUTINE GENER.
         NCD.GT.0... P WAVE.
         NCD.LT.0... S WAVE.
   II...     SUCCESSIVE NUMBER OF THE RECEIVER POSITION.
   T...      ARRIVAL TIME.
   AX...     AMPLITUDE OF THE RADIAL COMPONENT.
   AY...     AMPLITUDE OF THE TRANSVERSE COMPONENT.
   AZ...     AMPLITUDE OF THE VERTICAL COMPONENT.
   PHX...    PHASE SHIFT OF THE RADIAL COMPONENT.
   PHY...    PHASE SHIFT OF THE TRANSVERSE COMPONENT.
   PHZ...    PHASE SHIFT OF THE VERTICAL COMPONENT.
   PNEW...   THE PARAMETER OF THE RAY (INITIAL ANGLE).
   TAST...   THE GLOBAL ABSORPTION FACTOR T*.

THE DATA 5 ARE REPEATED AT LU2 MANY TIMES, DEPENDING ON THE
NUMBER OF ELEMENTARY WAVES AND ON THE NUMBER OF RECEIVER POSITIONS.

*****************************************************************

A SHORT DESCRIPTION OF THE UNIVERSAL 2-D RAY TRACING ROUTINE
RAY2(X0,Z0,T0,FI0,X,Z,T,FI,TMAX,DT,AC)
**************************************

THE ROUTINE RAY2 COMPUTES THE RAY SPECIFIED BY INITIAL CONDITIONS
X0,Z0,T0,FI0. IT RETURNS THE QUANTITIES X,Z,T,FI AT THE POINT OF
TERMINATION OF THE RAY. THE COMPUTATION IS CONTROLLED BY
THE QUANTITIES TMAX,DT,AC.
  THE MEANING OF THESE PARAMETERS IS AS FOLLOWS:

   X0,Z0...  COORDINATES OF THE INITIAL POINT OF THE RAY.
   T0...     INITIAL TIME.
   FI0...    INITIAL ANGLE.
   X,Z...    COORDINATES OF THE POINT AT WHICH THE COMPUTATION
             OF THE RAY IS TERMINATED.
   T...      CORRESPONDING TIME.
   FI...     CORRESPONDING ANGLE OF THE RAY.
   DT,AC...  THE SAME MEANING AS DESCRIBED IN THE INPUT
             DATA CARD NO. 12.
   TMAX...   THE MAXIMUM PERMITTED TRAVEL TIME ALONG THE RAY.
             IN THIS VERSION OF THE PROGRAM, TMAX=10000.

THE HISTORY OF THE RAY IS STORED IN COMMON/RAY/, SEE ABOVE. THE
QUANTITY IND IN THIS COMMON SPECIFIES THE REASON FOR THE TERMINATION
OF THE RAY.
  ROUTINES USED: RKGS,FCT,OUT,VELOC,ROOT.
  THE INFORMATION ABOUT THE MODEL IS TRANSMITTED TO THE RAY TRACING
ROUTINE BY COMMON/INTRF/ AND COMMON/DEN/. ROUTINE VELOC IS USED
TO DETERMINE THE VELOCITY AND ITS DERIVATIVES AS WELL AS QUALITY
FACTOR AT ANY ARBITRARY POINT IN THE MEDIUM.
  THE NUMERICAL CODE OF THE WAVE IS TRANSMITTED TO THE RAY TRACING
ROUTINE BY COMMON/COD/.

****************************************************************

A SHORT DESCRIPTION OF THE UNIVERSAL TWO-POINT RAY TRACING
ROUTINE
TRAMP(XSOUR,ZSOUR,TSOUR,DT,AC,TMAX,ASTOP,ITMAX,MOUT,MDIM,MPSOUR,
MSSOUR,NCODE,LU1,LU2,METHOD,ISOUR,LTAPE,ITPR)
*********************************************

THE ROUTINE TRAMP IS DESIGNED FOR A TWO-POINT RAY TRACING. IT
COMPUTES THE RAYS WHICH ARRIVE AT A SYSTEM OF RECEIVERS DISTRIBUTED
REGULARLY OR IRREGULARLY ALONG A PROFILE OF RECEIVERS. THE RECEIVERS
CAN BE SITUATED ON THE EARTH'S SURFACE, ALONG INTERFACES, OR ALONG
VERTICAL PROFILES. THE MODIFIED SHOOTING METHOD IS USED. THE RAYS
TERMINATING AT SPECIFIED RECEIVERS ARE SOUGHT BY ITERATIVE USE OF
INITIAL-VALUE RAY TRACING. EACH RAY IS SPECIFIED BY ITS INITIAL
ANGLE AT THE SOURCE. THE BASIC SYSTEM OF THE INITIAL ANGLES OF
RAYS AT THE SOURCE IS SPECIFIED IN THE INPUT DATA CARD NO.12. TO
BE SURE THAT NO RECEIVERS ARE OUTSIDE THE RANGE OF INITIAL ANGLES,
THE WHOLE RANGE OF ANGLES (E.G. FROM 3.1416 TO -3.1416) SHOULD BE
USED. TO SAVE THE COMPUTING TIME, THE RANGE MAY BE SELECTED NARROWER.
FROM THE COMPUTED TERMINATION POINTS OF THESE RAYS, THE RAYS ARRIVING
AT SPECIFIED RECEIVERS ARE ITERATIVELY DETERMINED. THE TRAVEL TIME
CURVE MAY BE ARBITRARY, WITH LOOPS, CAUSTIC, SHADOWS, ETC. A SPECIAL
CARE IS DEVOTED TO CERTAIN IRREGULAR SITUATIONS IN THE RAY FIELD.
THE ROUTINE IS AUTOMATICALLY LOOKING FOR ALL BOUNDARY RAYS
(I.E. THE RAYS CORRESPONDING TO THE BOUNDARIES OF SHADOW ZONES AT
THE EARTHS SURFACE, TO THE LEFT-HAND SIDE AND RIGHT-HAND SIDE UPPER
CORNERS OF THE MODEL, ETC.). A SPECIAL CARE IS ALSO DEVOTED TO
THE PROBLEM OF SLIGHTLY REFLECTED WAVES (WHICH HAVE RAYS SIMILAR TO
HEAD WAVES). AS THE ROUTINE USES ONLY A SINGLE PRECISION, ALL THE
RAYS OF THIS WAVE, ESPECIALLY THOSE RAYS WHICH ARE VERY CLOSE TO
THE CRITICAL RAY, CANNOT BE SOMETIMES DETERMINED.

  THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS:

   XSOUR,ZSOUR,TSOUR,DT,AC,ASTOP... SEE INPUT CARD NO. 11 AND 12
             (WITH ASTOP=ASP).
   TMAX...   MAXIMUM TRAVEL TIME PERMITTED. IN SEIS, TMAX=10000.
   ITMAX,MOUT,MPSOUR,MSSOUR,METHOD,LU1,LU2,LTAPE... SEE INPUT
             DATA CARD NO.9.
   NCODE...  INTERNAL WAVE CODE. SEE DETAILS IN THE DESCRIPTION
             OF ROUTINE GENER.
   ISOUR...  NUMBER OF THE LAYER, IN WHICH THE SOURCE IS SITUATED.
   ITPR...   FOR BURIED RECEIVERS, THE NUMBER OF THE INTERFACE
             ON WHICH THE RECEIVERS ARE SITUATED.

THE INFORMATION ABOUT THE RECEIVER POSITIONS, ABOUT THE BASIC
SYSTEM OF INITIAL ANGLES FOR RAY SHOOTING, AND ABOUT THE REQUIRED
ACCURACY OF COMPUTATIONS (SEE REPS, INPUT CARD NO.11) IS TRANSFERRED
TO TRAMP BY COMMON /DIST/.
  ROUTINES USED: RAY2, AMPL,JPAR.
THE ROUTINE GENERALLY DOES NOT NEED ANY INFORMATION ON THE VELOCITY
DISTRIBUTION IN THE MODEL AND ON THE HISTORY OF RAYS. IT USES ONLY
THE INFORMATION ON THE X-COORDINATE OF THE TERMINATION POINT FOR
RECEIVERS DISTRIBUTED ALONG THE EARTH'S SURFACE OR ALONG AN INTERFACE,
OR ON THE Z-COORDINATE OF THE TERMINATION POINT FOR RECEIVERS
DISTRIBUTED ALONG A VERTICAL PROFILE, AND THE INFORMATION ON THE
INITIAL ANGLES OF RAYS AND ON THE PARAMETER IND. IN OTHER WORDS,
THE ROUTINE CAN WORK WITH ANY OTHER INITIAL VALUE RAY TRACING ROUTINE
WHICH PROVIDES THESE QUANTITIES. IN FACT, IT HAS BEEN SUCCESSFULLY
APPLIED IN VARIOUS OTHER PROGRAMS, EVEN WHERE THE RAY PARAMETER HAD
A DIFERENT MEANING THAN THE INITIAL ANGLE.
  IN THE ROUTINE TRAMP PARTS OF TWO FILES ARE GENERATED. THE FIRST
IS THE FILE LU1 (RAY DIAGRAMS TO SPECIFIED RECEIVERS, TRAVEL TIMES
AND AMPLITUDES), AND THE SECOND IS THE FILE LU2 (THE DATA FOR COMPUTATION
OF RAY SYNTHETIC SEISMOGRAMS). THE CONTENT OF FILES IS DESCRIBED ABOVE.
LIKEWISE, THE OUTPUT ON THE LINE PRINTER IS DESCRIBED ABOVE.

******************************************************************

A SHORT DESCRIPTION OF THE UNIVERSAL ROUTINE FOR 2-D DYNAMIC
RAY TRACING ALONG A KNOWN RAY
JPAR(Q0,P0,Q,P,IPRINT,KMAH)
***************************

WHEN THE RAY IS KNOWN AND ITS HISTORY IS STORED IN COMMON/RAY/,
2-D DYNAMIC RAY TRACING CAN BE PERFORMED ALONG THE RAY BY
THE PROCEDURE JPAR.
  THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS:

   Q0,P0...  INITIAL VALUES FOR THE DYNAMIC RAY TRACING
             SYSTEM. IN THE PRESENT VERSION Q0=0 AND P0=1./V0,
             WHERE V0 IS THE CORRESPONDING WAVE VELOCITY AT THE
             SOURCE.
   Q,P...    THE SOLUTIONS OF THE DYNAMIC RAY TRACING SYSTEM
             AT THE TERMINATION POINT OF THE RAY. THE QUANTITY Q
             REPRESENTS THE SPREADING FUNCTION FOR A LINE SOURCE.
             THE SPREADING FUNCTION FOR THE POINT SOURCE CAN BE
             EASILY OBTAINED BY WELL-KNOWN METHODS.
   IPRINT..  CONTROLS THE PRINTING OF DYNAMIC RAY TRACING ALONG
             THE RAY ON A LINE PRINTER. IN SEIS83, NO PRINT OF
             DYNAMIC RAY TRACING RESULTS IS PERFORMED, IPRINT=0.
   KMAH...   GIVES A NUMBER OF CHANGES OF SIGN OF THE QUANTITY
             Q ALONG THE RAY (IN OTHER WORDS: IT SPECIFIES HOW
             MANY TIMES THE RAY TOUCHED A CAUSTIC). IT IS USED
             TO DETERMINE THE PHASE SHIFT CAUSED BY CAUSTICS.

THE NUMERICAL CODE OF THE WAVE IS TRANSMITTED TO THE ROUTINE
JPAR BY COMMON/COD/.
ROUTINES USED: RK,FCTA.
  NOTE THAT ALL THE INFORMATION ABOUT THE MODEL USED IN JPAR,RK
AND FCTA IS TRANSMITTED TO THESE ROUTINES ONLY BY COMMON/RAY/.
IN OTHER WORDS, THE DYNAMIC RAY TRACING PACKAGE CAN ALSO BE
USED WITH OTHER RAY TRACING ROUTINES, WHEN THEY GENERATE
COMMON/RAY/ AS DESCRIBED ABOVE.

****************************************************************

A SHORT DESCRIPTION OF THE UNIVERSAL ROUTINE FOR THE EVALUATION
OF RAY AMPLITUDES OF SEISMIC BODY WAVES ALONG A KNOWN RAY
AMPL(FJ,AX,AY,AZ,PHX,PHY,PHZ)
*****************************

THE ROUTINE AMPL COMPUTES THE RAY AMPLITUDES AND PHASE SHIFTS
OF ARBITRARY MULTIPLY-REFLECTED (POSSIBLY CONVERTED) SEISMIC
BODY WAVES IN A 2-D LATERALLY INHOMOGENEOUS MEDIA WITH CURVED
INTERFACES. THE COMPUTATION IS PERFORMED ALONG A KNOWN RAY, AND
ALL THE INFORMATION ABOUT THE RAY IS TRANSFERRED TO AMPL BY
COMMON/RAY/. IT IS ALSO ASSUMED THAT THE GEOMETRICAL SPREADING
IS KNOWN (E.G., IT COULD BE DETERMINED INDEPENDENTLY BY THE
PROCEDURE JPAR, SEE ABOVE). THE GEOMETRICAL SPREADING MAY
CORRESPOND TO A POINT SOURCE OR A LINE SOURCE.
  THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS:

   FJ...     GEOMETRICAL SPREADING.
   AX,AY,AZ.. AMPLITUDES OF RADIAL, TRANSVERSE AND VERTICAL
             COMPONENTS.
   PHX,PHY,PHZ...PHASE SHIFTS OF RADIAL, TRANSVERSE AND VERTICAL
             COMPONENTS.

THE NUMERICAL CODE OF THE WAVE IS TRANSFERRED TO THE ROUTINE
AMPL BY COMMON/COD/.
  ROUTINES USED: COEF8.
NOTE THAT THE MODEL INFORMATION USED IN AMPL IS TRANSFERRED TO
THIS ROUTINE BY COMMON/RAY/ AND /SOUR/. IN OTHER WORDS, ROUTINE
AMPL CAN BE USED WITH ANY OTHER RAY TRACING ROUTINE, WHEN THE
ROUTINE GENERATES THESE COMMONS, AND WHEN THE GEOMETRICAL SPREADING
IS ALSO DETERMINED IN ADVANCE.

*****************************************************************

A SHORT DESCRIPTION OF THE UNIVERSAL ROUTINE FOR DETERMINATION
OF REFLECTION/TRANSMISSION COEFFICIENTS
COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH)
**************************************************

THE ROUTINE COEF8 IS A STANDARD ROUTINE FOR THE COMPUTATION OF
P, SV AND SH REFLECTION/TRANSMISSION COEFFICIENTS OF PLANE WAVES
AT A PLANE INTERFACE BETWEEN TWO HOMOGENEOUS MEDIA OR AT A FREE
SURFACE OF A HOMOGENEOUS HALFSPACE. IN THE SECOND CASE, THE
CONVERSION COEFFICIENTS AT RECEIVERS CAN ALSO BE OPTIONALLY COMPUTED.
CONVERSION COEFFICIENTS (COMPONENTS OF A CONVERSION VECTOR) ARE
DETERMINED WITH RESPECT TO A LOCAL CARTESIAN COORDINATE SYSTEM,
DEFINED AS FOLLOWS: POSITIVE X-AXIS POINTS ALONG THE PROFILE TO
THE RIGHT, POSITIVE Z-AXIS POINTS UP INTO THE FREE SPACE, THE
Y-AXIS IS CHOSEN SUCH AS TO MAKE THE SYSTEM RIGHT-HANDED (IT POINTS
BEHIND THE PROFILE).

THE CODES OF INDIVIDUAL COEFFICIENTS ARE SPECIFIED BY THE
FOLLOWING NUMBERS
A/ INTERFACE BETWEEN TWO SOLID HALFSPACES
   P1P1....1   P1SV1....2   P1P2....3   P1SV2....4   SH1SH1...9
   SV1P1...5   SV1SV1...6   SV1P2...7   SV1SV2...8   SH1SH2..10
B/ FREE SURFACE (FOR RO2.LT.0.00001)
   PP.....1       PX.....5       PX,PZ...X- AND Z- COMPONENTS OF THE
   PSV....2       PZ.....6       COEF.OF CONVERSION,INCIDENT P WAVE.
   SVP....3       SVX....7       SVX,SVZ...X- AND Z- COMPONENTS OF
   SVSV...4       SVZ....8       COEF.OF CONVERSION,INCIDENT SV WAVE.
                                 SHY...Y-COMPONENT OF SH COEF.OF CONV.

I N P U T   P A R A M E T E R S
   P...RAY PARAMETER
   VP1,VS1,RO1..PARAMETERS OF THE FIRST HALFSPACE
   VP2,VS2,RO2..PARAMETERS OF SECOND HALFSPACE. FOR THE FREE
                SURFACE TAKE RO2.LT.0.00001,EG.RO2=0., AND
                ARBITRARY VALUES OF VP2 AND VS2
   NCODE...     CODE OF THE COMPUTED COEFFICIENT
             ND=0  WHEN THE POSITIVE DIRECTION OF THE RAY
                AND THE X-AXIS MAKE AN ACUTE ANGLE
             ND=1  WHEN THE WAVE IMPINGES ON THE INTERFACE
                AGAINST THE POSITIVE DIRECTION OF THE X-AXIS

O U T P U T   P A R A M E T E R S
   RMOD,RPH...MODUL AND ARGUMENT OF THE COEFFICIENT

N O T E S
1/ POSITIVE P...IN THE DIRECTION OF PROPAGATION.
2/ POSITIVE SV..TO THE LEFT FROM P.
3/ TIME FACTOR OF INCIDENT WAVE ... EXP(-I*OMEGA*T)
4/ FORMULAE ARE TAKEN FROM CERVENY ,MOLOTKOV, PSENCIK, RAY METHOD
   IN SEISMOLOGY, PAGES 30-35. DUE TO THE ABOVE DEFINITION OF
   THE LOCAL CARTESIAN COORDINATE SYSTEM, IN WHICH THE CONVERSION
   COEFFICIENTS ARE EVALUATED, THE SIGNS OF CERTAIN COEFFICIENTS
   ARE OPPOSITE.
***************************************************************

A SHORT DESCRIPTION OF THE ROUTINE FOR THE GENERATION OF NUMERICAL
CODES OF ELEMENTARY WAVES,
GENER(N,IDP,IDS,IBP,IBS,IUP,IUS,IREAD,IS,IR,NS,IFIN,NCODE)
**********************************************************

THE ROUTINE GENER GENERATES AUTOMATICALLY THE NUMERICAL CODES
OF DIRECT AND PRIMARY REFLECTED P, S AND CONVERTED WAVES. THE
CONVERSIONS CAN TAKE PLACE ONLY AT THE REFLECTION POINT. IN
ADDITION, ANY OTHER MULTIPLY REFLECTED WAVES (POSSIBLY CONVERTED
ARBITRARY NUMBER OF TIMES) CAN BE ADDED MANUALLY, USING IREAD.NE.0.
  THE SOURCE MAY BE SITUATED IN AN ARBITRARY LAYER (NO. IS), AND
SIMILARLY THE RECEIVER MAY BE ALSO SITUATED  IN AN ARBITRARY
LAYER (NO.IR).
  THE GENERATION STARTS WITH NAWX=0. THEN THE NUMERICAL CODE OF
THE FIRST WAVE IS GENERATED. AFTER THIS, NAWX MUST BE PUT =1
IN THE CALLING PROGRAM, AND ROUTINE GENERATES SUCCESSIVELY NEW
CODES AT EACH CALLING.
  THE SYSTEM OF GENERATED WAVES IS CONTROLLED BY THE PARAMETERS
IDP,IDS,IBP,IBS,IUP,IUS. FOR IDP,IDS,IBP,IBS, SEE EXPLANATIONS
IN THE DESCRIPTION OF THE INPUT DATA CARD NO.9. THE INDICES
IUP AND IUS CONTROL THE GENERATION OF PRIMARY REFLECTIONS WITH
THE REFLECTION POINTS ABOVE THE SOURCE AND THE RECEIVER.
THE QUANTITY NCODE IS THE INTERNAL WAVE CODE, WHICH IS INTRODUCED
IN ADDITION TO THE FULL NUMERICAL CODE (CALLED HERE ALSO EXTERNAL
WAVE CODE).THE INTERNAL WAVE CODE IS JUST ONE INTEGER, WHICH
GIVES THE SUCCESSIVE NUMBER OF THE ELEMENTARY WAVE, AS IT WAS
GENERATED BY THE ROUTINE GENER. THE FIRST GENERATED WAVE HAS
INTERNAL CODE 1, ETC. THIS INTERNAL WAVE CODE HAS AN IMPORTANCE
IN THE COMPUTATION OF SYNTHETIC SEISMOGRAMS IN THE PROGRAM SYNTPL,
AS ELEMENTARY WAVES CAN BE SELECTED FOR THE CONSTRUCTION OF THESE
SEISMOGRAMS.
  THE INDEX NLAYER GIVES A FULL NUMBER OF LAYERS OF THE MODEL FOR
WHICH THE NUMERICAL CODES OF ELEMENTARY WAVES ARE GENERATED.
IT MIGHT BE, OF COURSE, LESS THAN THE REAL NUMBER OF LAYERS IN
THE MODEL USED FOR THE COMPUTATION. FOR EXAMPLE, THE INDEX NLAY
(SEE INPUT CARD NO.9) CAN BE USED TO DECREASE ARTIFICIALLY THE
NUMBER OF LAYERS IN THE MODEL USED FOR GENERATION.
FOR IREAD=0, ONLY THE AUTOMATICALLY GENERATED DIRECT AND
PRIMARY REFLECTED WAVES ARE USED. FOR IREAD=1, ARBITRARY NUMBER
OF OTHER ELEMENTARY WAVES (MULTIPLE REFLECTIONS) CAN BE OPTIONALLY
ADDED, SEE THE DESCRIPTION OF MANUAL GENERATION IN INPUT CARD
NO.9 AND 13.
  THE PARAMETER IFIN=0 WHEN THE ROUTINE RETURNS A NEW NUMERICAL CODE.
THE END OF GENERATION IS SHOWN BY IFIN=1.
  THE REFRACTED WAVES ARE AUTOMATICALLY INCLUDED IN THE NUMERICAL
CODES CORRESPONDING TO REFLECTED WAVES, AS WAVES WITH COMPOUND
ELEMENTS, SEE DETAILS IN INPUT DATA CARD NO.13.
  NOTE THAT IN THE LAYER CONTAINING THE SOURCE, THE WAVES, WHICH
STRIKE FIRST THE INTERFACE ABOVE THE SOURCE ARE COMPUTED SEPARATELY
FROM THOSE, WHICH STRIKE FIRST THE INTERFACE BELOW THE SOURCE.
  THE ROUTINE IS QUITE UNIVERSAL, IT CAN BE USED WITH ANY OTHER
PROGRAM WHICH WORKS WITH THE SAME NUMERICAL CODES OF ELEMENTARY
WAVES AND IN WHICH THE INTERFACES CROSS THE WHOLE MODEL, FROM
THE LEFT TO THE RIGHT.

****************************************************************

A SHORT DESCRIPTION OF THE ROUTINE FOR THE APPROXIMATION OF
INTERFACES AND VELOCITY DISTRIBUTION WITHIN INDIVIDUAL LAYERS
MODEL(MTEXT,LU,MM)
******************

***********************************************************
BICUBIC SPLINE INTERPOLATION OF VELOCITIES SPECIFIED IN THE
GRID POINTS OF A RECTANGULAR NETWORK:  MM=0
***********************************************************

THE ROUTINE MODEL CONSISTS OF TWO PRINCIPAL PARTS:
1) IN THE FIRST PART, THE INPUT DATA CONCERNING INTERFACES
   AND VELOCITY DISTRIBUTION INSIDE INDIVIDUAL LAYERS ARE
   READ IN. INTERFACES ARE INTERPOLATED BY CUBIC SPLINES
   (ROUTINE SPLIN), THE VELOCITY DISTRIBUTION BY BICUBIC
   SPLINES (ROUTINES BIAP AND SPLIN). THE COEFFICIENTS SO
   OBTAINED ARE STORED IN COMMONS /INTRF/, /MEDIM/ AND
   /VCOEF/. THE INPUT DATA CONCERNING THE DENSITY AND
   QUALITY FACTOR DISTRIBUTION ARE STORED IN COMMON/DEN/
   AND COMMON/ABSR/.

2) IN THE SECOND PART, CERTAIN TABLES DESCRIBING THE MODEL ARE
   OPTIONALLY PRINTED ON THE LINE PRINTER, SEE THE DETAILED
   DESCRIPTION IN THE SECTION ON OUTPUT TABLES.

THE MEANING OF THE INPUT PARAMETERS MTEXT, LU AND MM IS EXPLAINED IN
INPUT DATA CARD NO.1. NOTE THAT MM=0 FOR THIS ROUTINE.
  ROUTINES USED: BIAP, SPLIN, VELOC. THE ROUTINE VELOC (SEE ITS
DETAILED DESCRIPTION FOLLOWING) IS USED ONLY IN THE SECOND
PART OF THE PROGRAM, FOR THE CONSTRUCTION OF THE PRINTER PLOT.
THE QUANTITIES STORED IN COMMON/AUXX/AND THE QUANTITIES LAY,
ITYPE,NDER FROM COMMON/AUXI/ ARE USED FOR COMMUNICATION BETWEEN
THE ROUTINES MODEL, BIAP, SPLIN, VELOC, SEE THE DESCRIPTION OF
COMMONS /AUXI/ AND /AUXX/.

****************************************************************

A SHORT DESCRIPTION OF THE ROUTINE FOR THE DETERMINATION OF
VELOCITY AND ITS DERIVATIVES AT AN ARBITRARY POINT OF THE MEDIUM
VELOC(Y,VEL)
************

***********************************************************
BICUBIC SPLINE INTERPOLATION OF VELOCITIES SPECIFIED IN THE
GRID POINTS OF A RECTANGULAR NETWORK:  MM=0
***********************************************************

IN THIS ROUTINE, THE SPECIFIED POINT AT WHICH THE VELOCITY
IS TO BE DETERMINED IS FIRST LOCALIZED IN THE GIVEN VELOCITY
NETWORK. THEN THE VALUES OF VELOCITY AND OF ITS FIRST
PARTIAL DERIVATIVES ARE DETERMINED BY EVALUATING BICUBIC
SPLINE POLYNOMIALS. WHEN MDIM.GT.1, THE SECOND DERIVATIVES
ARE ALSO DETERMINED.
  THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS:

   Y(1),Y(2)...X AND Z COORDINATES OF THE POINT AT WHICH THE
             VELOCITY IS TO BE DETERMINED.
   Y(3),Y(4)...VALUES OF THE COMPONENTS OF THE SLOWNESS VECTOR.
             THEY HAVE ONLY A FORMAL MEANING IN THE ROUTINE.
   VEL(1)... VELOCITY AT THE POINT Y(1),Y(2).
   VEL(2),VEL(3)... FIRST PARTIAL DERIVATIVES VX AND VZ
   VEL(4),VEL(5),VEL(6)... SECOND PARTIAL DERIVATIVES VXX, VXZ
             AND VZZ.

  NOTE THAT FOR THE PROPER OPERATION OF ROUTINE VELOC, ALL THE
PARAMETERS OF THE COMMON/AUXX/ MUST BE SPECIFIED A PRIORI
(SEE THE DESCRIPTION OF COMMON/AUXX/). ALSO THE PARAMETERS
LAY,ITYPE AND NDER,(SEE THE DESCRIPTION OF COMMON/AUXI/), MUST
BE SPECIFIED BEFORE THE USE OF ROUTINE VELOC.

**************************************************************

A SHORT DESCRIPTION OF THE ROUTINE FOR THE APPROXIMATION OF
INTERFACES AND VELOCITY DISTRIBUTION WITHIN INDIVIDUAL LAYERS
MODEL(MTEXT,LU,MM)
******************

***************************************************
LINEAR VERTICAL INTERPOLATION OF VELOCITIES BETWEEN
ISOVELOCITY INTERFACES:  MM=1
***************************************************

THE ROUTINE MODEL CONSISTS OF TWO PRINCIPAL PARTS:
1) IN THE FIRST PART, THE INPUT DATA CONCERNING INTERFACES
   AND VELOCITY DISTRIBUTION INSIDE INDIVIDUAL LAYERS ARE
   READ IN. INTERFACES ARE INTERPOLATED BY CUBIC SPLINES
   (ROUTINE SPLIN). THE COEFFICIENTS THUS OBTAINED ARE STORED
   IN COMMON /INTRF/. FOR EACH LAYER TWO VELOCITIES ARE SPECIFIED,
   ONE CONSTANT ALONG THE UPPER INTERFACE, THE OTHER CONSTANT
   ALONG THE LOWER INTERFACE. THEY ARE STORED IN COMMON/MEDIM/.
   THE INPUT DATA CONCERNING THE DENSITY AND QUALITY FACTOR
   DISTRIBUTION ARE STORED IN COMMON/DEN/ AND COMMON/ABSR/.
2) IN THE SECOND PART, CERTAIN TABLES DESCRIBING THE MODEL ARE
   OPTIONALLY PRINTED ON THE LINE PRINTER, SEE THE DETAILED
   DESCRIPTION IN THE SECTION ON OUTPUT TABLES.

  THE MEANING OF THE INPUT PARAMETERS MTEXT, LU AND MM IS EXPLAINED
IN INPUT DATA CARD NO.1. LU, HOWEVER, HAS ONLY FORMAL MEANING HERE.
INTPUT DATA IN THIS ROUTINE ARE ALWAYS READ IN FROM CARDS OR
TERMINAL. NOTE THAT MM=1 FOR THIS ROUTINE.
  ROUTINES USED: SPLIN, VELOC. THE ROUTINE VELOC (SEE ITS
DETAILED DESCRIPTION FOLLOWING) IS USED ONLY IN THE SECOND
PART OF THE PROGRAM, FOR THE CONSTRUCTION OF THE PRINTER PLOT.
THE QUANTITIES LAY, ITYPE,NDER FROM COMMON/AUXI/ ARE USED
FOR COMMUNICATION BETWEEN THE ROUTINES MODEL, SPLIN, AND VELOC,
SEE THE DESCRIPTION OF COMMON /AUXI/.

**************************************************************

A SHORT DESCRIPTION OF THE ROUTINE FOR THE DETERMINATION OF
VELOCITY AND ITS DERIVATIVES AT AN ARBITRARY POINT OF THE
MEDIUM  VELOC(Y,VEL)
********************

***************************************************
LINEAR VERTICAL INTERPOLATION OF VELOCITIES BETWEEN
ISOVELOCITY INTERFACES:  MM=1
***************************************************

IN THIS ROUTINE, THE VALUES OF VELOCITY AND OF ITS FIRST
PARTIAL DERIVATIVES ARE DETERMINED BY INTERPOLATING LINEARLY
BETWEEN VELOCITY VALUES SPECIFIED ALONG THE UPPER AND LOWER
INTERFACE OF THE GIVEN LAYER. THE INTERPOLATION IS PERFORMED
ALONG A VERTICAL LINE PASSING THROUGH THE POINT IN WHICH
VELOCITY IS TO BE DETERMINED. WHEN MDIM.GT.1, THE SECOND
DERIVATIVES ARE ALSO DETERMINED.
  THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS
FOLLOWS:

   Y(1),Y(2)...X AND Z COORDINATES OF THE POINT AT WHICH THE
             VELOCITY IS TO BE DETERMINED.
   Y(3),Y(4)...VALUES OF THE COMPONENTS OF THE SLOWNESS VECTOR.
             THEY HAVE ONLY A FORMAL MEANING IN THE ROUTINE.
   VEL(1)... VELOCITY AT THE POINT Y(1),Y(2).
   VEL(2),VEL(3)... FIRST PARTIAL DERIVATIVES VX AND VZ
   VEL(4),VEL(5),VEL(6)... SECOND PARTIAL DERIVATIVES VXX, VXZ
             AND VZZ.

NOTE THAT FOR THE PROPER OPERATION OF ROUTINE VELOC, THE
PARAMETERS LAY, ITYPE AND NDER,(SEE THE DESCRIPTION OF
COMMON/AUXI/), MUST BE SPECIFIED BEFORE THE USE OF ROUTINE VELOC.

**************************************************************

A SHORT DESCRIPTION OF THE ROUTINE FOR THE APPROXIMATION OF
INTERFACES AND VELOCITY DISTRIBUTION WITHIN INDIVIDUAL LAYERS
MODEL(MTEXT,LU,MM)
******************

**********************************************************
BILINEAR INTERPOLATION OF VELOCITIES SPECIFIED IN THE GRID
POINTS OF A RECTANGULAR NETWORK:  MM=2
**********************************************************

THE ROUTINE MODEL CONSISTS OF TWO PRINCIPAL PARTS:
1) IN THE FIRST PART, THE INPUT DATA CONCERNING INTERFACES
   AND VELOCITY DISTRIBUTION INSIDE INDIVIDUAL LAYERS ARE
   READ IN. INTERFACES ARE INTERPOLATED BY CUBIC SPLINES
   (ROUTINE SPLIN), THE VELOCITY DISTRIBUTION BY BILINEAR
   INTERPOLATION (ROUTINE BILIN). THE COEFFICIENTS SO OBTAINED
   ARE STORED IN COMMONS /INTRF/, /MEDIM/ AND /VCOEF/.
   THE INPUT DATA CONCERNING THE DENSITY AND QUALITY FACTOR
   DISTRIBUTION ARE STORED IN COMMON/DEN/ AND COMMON/ABSR/.

2) IN THE SECOND PART, CERTAIN TABLES DESCRIBING THE MODEL ARE
   OPTIONALLY PRINTED ON THE LINE PRINTER,SEE THE DETAILED
   DESCRIPTION IN THE SECTION ON OUTPUT TABLES.

THE MEANING OF THE INPUT PARAMETERS MTEXT, LU AND MM IS EXPLAINED
IN INPUT DATA CARD NO.1. NOTE THAT MM=2 FOR THIS ROUTINE.
  ROUTINES USED: BIAP, SPLIN, VELOC. THE ROUTINE VELOC (SEE ITS
DETAILED DESCRIPTION FOLLOWING) IS USED ONLY IN THE SECOND
PART OF THE PROGRAM, FOR THE CONSTRUCTION OF THE PRINTER PLOT.
THE QUANTITIES STORED IN COMMON/AUXX/AND THE QUANTITIES LAY,
ITYPE,NDER FROM COMMON/AUXI/ ARE USED FOR COMMUNICATION BETWEEN
THE ROUTINES MODEL, BIAP, SPLIN, VELOC, SEE THE DESCRIPTION OF
COMMONS /AUXI/ AND /AUXX/.

****************************************************************

A SHORT DESCRIPTION OF THE ROUTINE FOR THE DETERMINATION OF
VELOCITY AND ITS DERIVATIVES AT AN ARBITRARY POINT OF THE MEDIUM
VELOC(Y,VEL)
************

**********************************************************
BILINEAR INTERPOLATION OF VELOCITIES SPECIFIED IN THE GRID
POINTS OF A RECTANGULAR NETWORK:  MM=2
**********************************************************

IN THIS ROUTINE, THE SPECIFIED POINT AT WHICH THE VELOCITY
IS TO BE DETERMINED IS FIRST LOCALIZED IN THE GIVEN VELOCITY
NETWORK. THEN THE VALUES OF VELOCITY AND OF ITS FIRST
PARTIAL DERIVATIVES ARE DETERMINED BY EVALUATING BILINEAR
POLYNOMIALS. WHEN MDIM.GT.1, THE SECOND DERIVATIVES ARE
ALSO DETERMINED.
  THE MEANING OF THE INPUT AND OUTPUT PARAMETERS IS AS FOLLOWS:

   Y(1),Y(2)...X AND Z COORDINATES OF THE POINT AT WHICH THE
             VELOCITY IS TO BE DETERMINED.
   Y(3),Y(4)...VALUES OF THE COMPONENTS OF THE SLOWNESS VECTOR.
             THEY HAVE ONLY A FORMAL MEANING IN THE ROUTINE.
   VEL(1)... VELOCITY AT THE POINT Y(1),Y(2).
   VEL(2),VEL(3)... FIRST PARTIAL DERIVATIVES VX AND VZ
   VEL(4),VEL(5),VEL(6)... SECOND PARTIAL DERIVATIVES VXX, VXZ
             AND VZZ.

NOTE THAT FOR THE PROPER OPERATION OF ROUTINE VELOC, ALL THE
PARAMETERS OF THE COMMON/AUXX/ MUST BE SPECIFIED A PRIORI
(SEE THE DESCRIPTION OF COMMON/AUXX/). ALSO THE PARAMETERS
LAY,ITYPE AND NDER,(SEE THE DESCRIPTION OF COMMON/AUXI/), MUST
BE SPECIFIED BEFORE THE USE OF ROUTINE VELOC.

**************************************************************