C FILE 'RAY.FOR' FOR THE COMPLETE TRACING OF A RAY FROM A GIVEN POINT
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C THIS FILE CONSISTS OF:
C     RAYB... BLOCK DATA SUBROUTINE DEFINING COMMON BLOCKS /RAYT/ AND
C             /DCRT/ TO STORE THE INPUT DATA FOR COMPLETE RAY TRACING.
C     DESCRIPTION OF THE QUANTITIES COMPUTED ALONG A RAY (SEE
C             C.R.T.5.2).
C     RAY1... SUBROUTINE DESIGNED TO READ THE INPUT DATA FOR COMPLETE
C             RAY TRACING AND TO STORE THEM IN THE COMMON BLOCKS /RAYT/
C             AND /DCRT/ (SEE C.R.T.5.6).
C     RAY2... SUBROUTINE DESIGNED TO CONTINUE IN THE COMPLETE RAY
C             TRACING OF A RAY FROM THE GIVEN POINT (SEE C.R.T.5.7).
C
C INPUT DATA FOR COMPLETE RAY TRACING:
C     THE DATA ARE READ IN BY THE LIST DIRECTED INPUT (FREE FORMAT).  IN
C     THE LIST OF INPUT DATA BELOW, EACH NUMBERED PARAGRAPH INDICATES
C     THE BEGINNING OF A NEW INPUT OPERATION (NEW READ STATEMENT).  IF
C     THE FIRST LETTER OF THE SYMBOLIC NAME OF THE INPUT VARIABLE IS
C     I-N, THE CORRESPONDING VALUE IN INPUT DATA MUST BE OF THE TYPE
C     INTEGER.  OTHERWISE (EXCEPT TEXTR), THE INPUT PARAMETER IS OF THE
C     TYPE REAL.
C (1) TEXTR
C     STRING DESCRIBING THE DATA.  ONLY THE FIRST 80 CHARACTERS OF THE
C     STRING ARE SIGNIFICANT.
C (2) KSTORE,NEXPS,NHLF,MODCRT
C     QUANTITIES CONTROLING NUMERICAL INTEGRATION.
C     KSTORE..SPECIFIES WHETHER THE CONVERSION COEFFICIENTS ARE TO BE
C             CONSIDERED (SEE C.R.T.5.5.4 AND 5.6E):
C             KSTORE.LE.1... NO AMPLITUDE CONVERSION COEFFICIENTS ARE
C               APPLIED AT THE POINT OF INTERSECTION WITH AN INTERFACE.
C             KSTORE.GE.2... AMPLITUDE CONVERSION COEFFICIENTS ARE
C               APPLIED AT THE POINT OF INTERSECTION WITH AN INTERFACE.
C             WHENEVER KSTORE IS CHANGED, DATA (8) AND (9) SHOULD BE
C             ADJUSTED PROPERLY.
C     NEXP... SPECIFIES INDEPENDENT VARIABLE ALONG RAYS, SEE
C             C.R.T.(5.1).
C     NHLF... MAXIMUM ALLOWED NUMBER OF HALVINGS (BISECTIONS) OF INITIAL
C             INCREMENT OF INDEPENDENT VARIABLE DURING NUMERICAL
C             INTEGRATION.
C     MODCRT..SPECIFIES THE MODE OF COMPLETE RAY TRACING, AND THE RAY
C             ELEMENTS ALONG WHICH THE COMPUTED QUANTITIES ARE STORED.
C             0...INITIAL MODE,
C               STORING THE COMPUTED QUANTITIES ALONG THE WHOLE RAYS,
C               NO STORING AT THE ENDPOINTS OF THE RAY ELEMENTS.
C             1...INITIAL MODE,
C               STORING THE QUANTITIES JUST ALONG NEW RAY ELEMENTS,
C               NO STORING AT THE ENDPOINTS OF THE RAY ELEMENTS.
C             2...INITIAL MODE,
C               STORING THE QUANTITIES JUST ALONG NEW RAY ELEMENTS,
C               STORING AT THE ENDPOINTS OF THE NEW RAY ELEMENTS.
C             3...INTERPOLATION MODE,
C               STORING THE QUANTITIES JUST ALONG NEW RAY ELEMENTS,
C               STORING AT THE ENDPOINTS OF THE NEW RAY ELEMENTS.
C               *** MODCRT=3 IS NOT ENABLED BY THE CURRENT VERSION OF
C               'INIT.FOR' ***.
C (3) STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT
C     QUANTITIES CONTROLING NUMERICAL INTEGRATION.
C     STORE...STEP OF INDEPENDENT VARIABLE FOR STORING THE COMPUTED
C             QUANTITIES ALONG A RAY (C.R.T.5.5.1).  FOR STORE=0 THE
C             QUANTITIES ARE NOT STORED ALONG RAYS.
C     STEP... INITIAL INCREMENT OF INDEPENDENT VARIABLE FOR NUMERICAL
C             INTEGRATION.
C     UEB...  UPPER ERROR BOUND OF TRAVEL TIME PER ONE STEP OF NUMERICAL
C             INTEGRATION.  ERRORS IN THE COORDINATES OF POINTS ALONG
C             THE RAY ARE APPROXIMATELY TRANSFORMED TO UNITS OF TRAVEL
C             TIME AND ARE ALSO BOUNDED BY UEB.  THE ERROR PER STEP OF
C             NUMERICAL INTEGRATION IS AUTOMATICALLY KEPT WITHIN THE
C             LIMIT UEB IF THIS DOES NOT REQUIRE MORE THAN NHLF
C             BISECTIONS OF THE INITIAL INCREMENT STEP.  IN THE OPPOSITE
C             CASE THE UPPER ERROR BOUND IS 2,4,8... TIMES GREATER FOR
C             THE COMPUTATION OF THE REST OF THE RAY.  THUS THE
C             COMPUTATION OF EACH RAY IS COMPLETED.
C     UEBPP,UEBPH,UEBHH... MAXIMUM ALLOWED ACCUMULATED DEVIATIONS OF THE
C             TWO POLARIZATION VECTORS FROM THE CONDITIONS OF
C             ORTHONORMALITY (C.R.T.5.8.2D AND 5.8.3I).  THE ACCUMULATED
C             DEVIATIONS ARE EXPRESSED IN TIME UNITS.  IF ANY OF THE
C             ACCUMULATED DEVIATIONS OVERRFLOWS ITS UPPER BOUND, A
C             WARNING IS GENERATED.
C     UEBDRT..MAXIMUM ALLOWED DEVIATION OF L.H.S. AND R.H.S. IN
C             C.R.T.,EQ.(5.11).  IF THE DEVIATION OF ANY COMPONENT OF
C             THE 4*4 MATRIX EXCEEDS UEBDRT, A WARNING IS GENERATED.
C (4) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,TMAX
C     THE BOUNDARIES OF THE COMPUTATIONAL VOLUME.
C     X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX... COORDINATES SPECIFYING THE
C             COORDINATE PLANES BOUNDING THE COMPUTATIONAL VOLUME.  THE
C             COORDINATE PLANES ARE INDEXED 101, 102, 103, 104, 105 AND
C             106.  DEFAULT VALUES ARE THE BOUNDARIES OF THE MODEL, SEE
C             SUBROUTINE MODEL1 AND THE INPUT DATA FOR THE MODEL.
C     TMAX... MAXIMUM TRAVEL TIME FOR THE COMPUTATION OF RAYS. THE
C             CORRESPONDING ISOCHRONE IS INDEXED BY 107.  DEFAULT VALUE
C             IS 999999.
C (5) NSRFCA
C     NUMBER OF AUXILIARY SURFACES, SEE C.R.T.5.3.  THE SURFACES ARE
C     INDEXED SEQUENTIALLY BY POSITIVE INTEGERS FROM NSRFC+1 TO
C     NSRFC+NSRFCA.  DEFAULT VALUE IS 0.
C (6) THE INDICES OF END SURFACES.
C     EACH INDEX MUST HAVE THE SIGN CORRESPONDING TO THE SIDE OF THE
C     SURFACE AT WHICH THE COMPUTATIONAL VOLUME IS SITUATED.
C     THE INDICES MAY BE SPECIFIED IN AN ARBITRARY ORDER AND MUST BE
C     TERMINATED BY A SLASH.  SEE C.R.T.5.4F.
C (7) THE INDICES OF SURFACES FOR STORING COMPUTED QUANTITIES.
C     THE INDICES OF SURFACES 1 TO NSRFC COVERING STRUCTURAL INTERFACES
C     MUST INCLUDE THE PROPER SIGN.  THE INDICES MAY BE SPECIFIED IN AN
C     ARBITRARY ORDER AND MUST BE TERMINATED BY A SLASH.  NOTE THAT THE
C     SAME SURFACE MAY BE SPECIFIED SEVERAL TIMES, ESPECIALLY IF THE
C     QUANTITIES AT DIFFERENT POINTS OF INTERSECTIONS WITH THE SURFACE
C     ARE STORED IN DIFFERENT FILES, SEE (8) AND (9) BELOW.
C     FOR DETAILS SEE C.R.T.5.5.2.
C (8) FOR EACH ABOVE SURFACE, THE SEQUENTIAL NUMBER OF THE FIRST POINT
C     OF INTERSECTION OF A RAY WITH THE SURFACE IN WHICH THE COMPUTED
C     QUANTITIES ARE STORED INTO THE CORRESPONDING FILE.  A RAY MAY HAVE
C     SEVERAL POINTS OF INTERSECTION WITH A SPECIFIED SURFACE.  THE
C     QUANTITIES NEED NOT BE STORED AT ALL POINTS OF INTERSECTION. THEY
C     MAY BE STORED JUST FROM THE POINT OF INTERSECTION SPECIFIED IN
C     THIS INPUT (8) TO THE ONE SPECIFIED IN THE INPUT (9) BELOW.
C     USUALLY, USER WILL WISH TO STORE THE QUANTITIES AT ALL POINTS OF
C     INTERSECTION, THAT CORRESPONDS TO THE VALUES OF 1 IN THIS INPUT
C     AND IS THE DEFAULT (INDICATED BY A SLASH IN THE INPUT DATA).
C     FOR A STRUCTURAL INTERFACE (NOT AUXILIARY SURFACE), THE POSITIVE
C     OR NEGATIVE SIDE OF THE SURFACE ISRFR IS SPECIFIED, THEN:
C       NO POINT OF INTERSECTION IS ACCOUNTED IF THE REFLECTED RAY IS
C         SITUATED AT THE OTHER SIDE OF THE SURFACE ISRFR THAN THE
C         SPECIFIED ONE, AND
C       ONE POINT OF INTERSECTION IS ACCOUNTED IF THE REFLECTED RAY IS
C         SITUATED AT THE SPECIFIED SIDE OF THE SURFACE ISRFR AND
C         KSTORE.GE.2 IN THE DATA (2).
C       TWO POINTS OF INTERSECTION ARE ACCOUNTED IF THE REFLECTED RAY IS
C         SITUATED AT THE SPECIFIED SIDE OF THE SURFACE ISRFR AND
C         KSTORE.LE.1 IN THE DATA (2).
C (9) FOR EACH ABOVE SURFACE, THE SEQUENTIAL NUMBER OF THE LAST POINT
C     OF INTERSECTION OF A RAY WITH THE SURFACE IN WHICH THE COMPUTED
C     QUANTITIES ARE STORED INTO THE CORRESPONDING FILE.  DEFAULT VALUES
C     OF THIS INPUT ARE 999999.
C (10) THE DATA SPECIFYING NSRFCA FUNCTIONS DESCRIBING THE AUXIIARY
C     SURFACES. THE DATA ARE READ BY SUBROUTINE SRFC1.  FOR THEIR
C     DESCRIPTION REFER TO SUBROUTINE SRFC1.
C     THIS INPUT IS PERFORMED ONLY IF NSRFCA.GT.0, SEE (5).
C
C STORAGE IN THE MEMORY:
C     THE INPUT DATA (1) TO (9) ARE STORED IN COMMON BLOCKS /RAYT/ AND
C     /DCRT/ DEFINED IN THE FOLLOWING SUBROUTINE:
C     ------------------------------------------------------------------
      BLOCK DATA RAYB
      CHARACTER*80 TEXTR
      COMMON/RAYT/TEXTR
      SAVE /RAYT/
      INTEGER MEND,MSTOR
      PARAMETER (MEND=128)
      PARAMETER (MSTOR=128)
      INTEGER KSTORE,NEXPS,NHLF,MODCRT
      REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7)
      INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR)
      COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH,
     *             UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR
      SAVE /DCRT/
      END
C     ------------------------------------------------------------------
C     TEXTR...THE NAME OF THE DATA. STRING OF 80 CHARACTERS.
C     KSTORE,NEXPS,NHLF,MODCRT... INPUT DATA (2).
C     STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT... INPUT DATA (3).
C     BOUNDR..BOUNDARIES X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,TMAX OF THE
C             COMPUTATIONAL VOLUME, SEE INPUT DATA (4).
C     NSRFCA..NUMBER OF AUXILIARY SURFACES, SEE INPUT DATA (5).
C     NEND... NUMBER OF END SURFACES, SEE INPUT DATA (6).
C     KEND... CONTAINS THE INDICES OF END SURFACES, SEE INPUT DATA (6).
C     NSTOR...NUMBER OF SURFACES FOR STORING COMPUTED QUANTITIES, SEE
C             INPUT DATA (7).
C     KSTOR...CONTAINS THE INDICES OF SURFACES FOR STORING COMPUTED
C             QUANTITIES, SEE INPUT DATA (7), THE INPUT DATA (8) AND
C             FINALY, THE INPUT DATA (9).  THE TOTAL OF 3*NSTOR
C             INTEGERS.
C     COMMON BLOCK /DCRT/ IS ALSO INCLUDED IN SUBROUTINE FILES
C     'RAYCB.FOR', 'INIT.FOR', 'WRIT.FOR' AND 'SCROPC.FOR'.
C     ALL THE INPUT DATA ARE STORED SEQUENTIALLY IN THE SAME ORDER AS
C     THEY WERE READ.  THE ONLY EXCEPTION ARE LOCATIONS NEND AND NSTOR
C     WHICH ARE INSERTED WHEN READING THE INPUT DATA.  THE INDEX OF THE
C     LAST ALLOCATED NUMERIC UNIT OF ARRAY KEND IS NAMED MEND.  THE
C     INDEX OF THE LAST ALLOCATED NUMERIC UNIT OF ARRAY KSTOR IS NAMED
C     MSTOR.  IF MEND OR MSTOR IS CHANGED, IT MUST BE ADJUSTED IN ALL
C     SUBROUTINES WHICH INCLUDE COMMON BLOCK /DCRT/.
C
C DATE: 1994, JANUARY 15
C CODED BY LUDEK KLIMES
C
C=======================================================================
C
C DESCRIPTION OF THE QUANTITIES COMPUTED ALONG A RAY (SEE C.R.T.5.2)
C
C LOCAL QUANTITIES (C.R.T.5.5.4):
C     YL(1)=VP... VELOCITY OF P-WAVES AT THE POINT.
C     YL(2)=VS... VELOCITY OF S-WAVES AT THE POINT.
C     YL(3)=RO... DENSITY AT THE POINT.
C     YL(4)=V1,YL(5)=V2,YL(6)=V3... VELOCITY DERIVATIVES IN GENERAL
C             COORDINATES.
C
C QUANTITIES COMPUTED ALONG A RAY (C.R.T.5.2):
C BASIC QUANTITIES COMPUTED ALONG A RAY (C.R.T.5.2.1):
C     Y(1)... TRAVEL TIME.
C     Y(2)... IMAGINARY PART OF THE COMPLEX-VALUED TRAVEL TIME.
C     Y(3),Y(4),Y(5)... COORDINATES OF POINTS ALONG THE RAY.
C     Y(6),Y(7),Y(8)... COVARIANT COMPONENTS OF THE SLOWNESS VECTOR.
C     Y(9),Y(10),Y(11)... COVARIANT COMPONENTS OF THE POLARIZATION
C             VECTOR E1 PERPENDICALAR TO THE RAY.
C     ( Y(12),Y(16),Y(20),Y(24) )    RAY PROPAGATOR MATRIX (I.E. THE
C     ( Y(13),Y(17),Y(21),Y(25) ) ...MATRIX OF FUNDAMENTAL SOLUTIONS
C     ( Y(14),Y(18),Y(22),Y(26) )    OF THE DYNAMIC RAY TRACING SYSTEM.
C     ( Y(15),Y(19),Y(23),Y(27) )
C     Y(28) TO Y(NY), WHERE NY=27+NAMPL... NAMPL REAL QUANTITIES
C             REPRESENTING COMPLEX-VALUED VECTORIAL REDUCED AMPLITUDES.
C             THE VECTORIAL REDUCED AMPLITUDES ARE SPECIFIED IN THE
C             RAY-CENTRED COORDINATE SYSTEM.
C             P WAVE AT THE INITIAL POINT OF THE RAY,
C             P WAVE AT THE POINT UNDER CONSIDERATION:
C               NAMPL=2,
C               Y(28)=REAL(A33), Y(29)=AIMAG(A33).
C             P WAVE AT THE INITIAL POINT OF THE RAY,
C             S WAVE AT THE POINT UNDER CONSIDERATION:
C               NAMPL=4,
C               Y(28)=REAL(A13), Y(29)=AIMAG(A13),
C               Y(30)=REAL(A23), Y(31)=AIMAG(A23).
C             S WAVE AT THE INITIAL POINT OF THE RAY,
C             P WAVE AT THE POINT UNDER CONSIDERATION:
C               NAMPL=4,
C               Y(28)=REAL(A31), Y(29)=AIMAG(A31),
C               Y(30)=REAL(A32), Y(31)=AIMAG(A32).
C             S WAVE AT THE INITIAL POINT OF THE RAY,
C             S WAVE AT THE POINT UNDER CONSIDERATION:
C               NAMPL=8,
C               Y(28)=REAL(A11), Y(29)=AIMAG(A11),
C               Y(30)=REAL(A21), Y(31)=AIMAG(A21),
C               Y(32)=REAL(A12), Y(33)=AIMAG(A12),
C               Y(34)=REAL(A22), Y(35)=AIMAG(A22).
C     Y(NY+1) TO Y(35), WHERE NY=27+NAMPL... UNDEFINED.
C AUXILIARY QUANTITIES COMPUTED ALONG A RAY (C.R.T.5.2.2):
C     YY(1)...INDEPENDENT VARIABLE ALONG A RAY.
C     YY(2)=UEBRAY... UPPER ERROR BOUND FOR RAY TRACING, WHICH IS EQUAL
C             TO THE INPUT VALUE UEB AT THE INITIAL POINT OF THE RAY
C             (C.R.T.5.6J).  IT IS ALWAYS DOUBLED WHEN THE NUMERICAL
C             INTEGRATION REQUIRES MORE THAN NHLF BISECTIONS OF THE
C             INITIAL INCREMENT STEP (C.R.T.5.6G).  UEBRAY.GT.UEB AT
C             THE ENDPOINT OF THE RAY INDICATES A DECREASED ACCURACY
C             OF COMPUTATION.
C     YY(3)=ERRPP,YY(4)=ERRPH,YY(5)=ERRHH... DEVIATIONS (IN ABSOLUTE
C             VALUES OF THE TWO COMPUTED BASIS VECTORS OF THE
C             RAY-CENTRED COORDINATE SYSTEM FROM THE CONDITIONS OF
C             ORTHONORMALITY (C.R.T.5.8.2D AND 5.8.3I), ACCUMULATED
C             ALONG THE RAY.  ANY OF THEM MAY BE COMPARED WITH THE
C             CORRESPONDING SPECIFIED LIMIT UEBPP, UEBPH, UEBHH AT THE
C             ENDPOINT OF THE RAY.
C     IY(1)=NY=27+NAMPL... NUMBER  OF  THE  BASIC  QUANTITIES
C             DESCRIBING THE POINT OF A RAY, SEE ABOVE.
C     IY(2)=KODIND... POSITION IN THE CODE (INDEX IN ARRAY KODE)
C             CORRESPONDING TO THE CONSIDERED ELEMENT OF A RAY. ITS
C             VALUE IS DETERMINED BY SUBROUTINE CODE (C.R.T.4.3).
C     IY(3)=ICB0... INDEX OF THE COMPLEX BLOCK, FROM WHICH THE RAY
C             ENTERED THE COMPLEX BLOCK IN WHICH THE COMPUTED
C             ELEMENT OF THE RAY IS SITUATED.
C             IY(3)=0 BEFORE LEAVING THE COMPLEX BLOCK, IN WHICH THE
C             INITIAL POINT OF THE RAY IS SITUATED.
C     IY(4)=ISB1... INDEX OF THE SIMPLE BLOCK CONTAINING THE
C             COMPUTED ELEMENT OF THE RAY.
C     IY(5)=ICB1... INDEX OF THE COMPLEX BLOCK CONTAINING THE
C             COMPUTED ELEMENT OF THE RAY, SUPPLEMENTED BY A SIGN '+'
C             FOR P WAVE AND SIGN '-' FOR S WAVE.
C     IY(6)=ISRF... INDEX OF THE SURFACE AT WHICH THE ENDPOINT OF THE
C             COMPUTED ELEMENT OF THE RAY IS SITUATED, SUPPLEMENTED BY
C             A SIGN '+' OR '-' FOR THE ENDPOINT SITUATED AT THE
C             POSITIVE OR NEGATIVE SIDE OF THE SURFACE, RESPECTIVELY.
C             ZERO INSIDE THE COMPLEX BLOCK.  NOTE: THE SIGN OF THIS
C             QUANTITY IS THE EXEPTION TO THE DEFINITION IN THE ORIGINAL
C             PAPER ON C.R.T.
C     IY(7)=ISB2... INDEX OF THE SIMPLE BLOCK TOUCHING THE COMPLEX
C             BLOCK ICB1 FROM THE OTHER SIDE OF THE SURFACE ISRF AT
C             THE ENDPOINT OF THE COMPUTED ELEMENT OF THE RAY.
C             ISB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF ISRF.
C             UNDEFINED INSIDE THE COMPLEX BLOCK, DEFINED ONLY AT
C             THE ENDPOINT OF THE ELEMENT OF THE RAY.
C     IY(8)=ICB2... INDEX OF THE COMPLEX BLOCK TOUCHING THE COMPLEX
C             BLOCK ICB1 FROM THE OTHER SIDE OF THE SURFACE ISRF AT
C             THE ENDPOINT OF THE COMPUTED ELEMENT OF THE RAY.
C             ICB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF ISRF.
C             UNDEFINED INSIDE THE COMPLEX BLOCK, DEFINED ONLY AT
C             THE ENDPOINT OF THE ELEMENT OF THE RAY.
C     IY(9)=IFCT... NUMBER OF INVOCATIONS OF SUBROUTINE FCT EVALUATING
C             THE RIGHT-HAND SIDES OF THE ORDINARY DIFFERENTIAL
C             EQUATIONS, ALONG THE COMPUTED PART OF THE RAY.
C     IY(10)=IOUTP... NUMBER OF SUCCESSFUL STEPS OF THE NUMERICAL
C             INTEGRATION, ALONG THE RAY (I.E. THE NUMBER OF INVOCATIONS
C             OF SUBROUTINE OUTP DECREASED BY THE NUMBER OF INVOCATIONS
C             OF SUBROUTINE HPCG).
C     IY(11)=ITRANS... NUMBER OF TRANSFORMATIONS AT AN INTERFACE (I.E.
C             THE NUMBER OF INVOCATIONS OF SUBROUTINE TRANS).
C     IY(12)=KMAH... NUMBER OF CAUSTIC POINTS ALONG THE RAY (THE INDEX
C             OF THE RAY TRAJECTORY).
C
C DATE: 1990, JANUARY 23
C WRITTEN BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C=======================================================================
C
      SUBROUTINE RAY1(LUN)
      INTEGER LUN
C
C SUBROUTINE RAY1 READS THE INPUT DATA FOR COMPLETE RAY TRACING (SEE
C C.R.T.5.6) AND STORES THEM IN COMMON BLOCKS /RAYT/ AND /DCRT/.  IT IS
C CALLED ONCE, BEFORE THE COMPLETE TRACING OF INDIVIDUAL RAYS BEGINS.
C SUBROUTINE RAY1 CALLS SUBROUTINE SRFC1 IN ORDER TO READ THE INPUT DATA
C FOR THE AUXILIARY SURFACES, AND TO STORE THEM IN THE MEMORY.
C
C INPUT:
C     LUN...  LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE
C             CONTAINING THE INPUT DATA.
C THE INPUT PARAMETER IS NOT ALTERED.
C
C NO OUTPUT.
C
C COMMON BLOCK /MODELC/ (SUBROUTINE FILE 'MODEL.FOR'):
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RAYT/:
      CHARACTER*80 TEXTR
      COMMON/RAYT/TEXTR
C COMMON BLOCK /DCRT/:
      INTEGER MEND,MSTOR
      PARAMETER (MEND=128)
      PARAMETER (MSTOR=128)
      INTEGER KSTORE,NEXPS,NHLF,MODCRT
      REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7)
      INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR)
      COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH,
     *             UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR
C ALL THE STORAGE LOCATIONS OF THE COMMON BLOCKS ARE DEFINED IN THIS
C SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      INTEGER NSRFC
      EXTERNAL RAYB,NSRFC,SRFC1
C     RAYB... BLOCK DATA SUBROUTINE OF THIS FILE.
C     NSRFC...FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'.
C     SRFC1 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C
C ERROR MESSAGES:
C     561...  INSUFFICIENT MEMORY IN /DCRT/:
C             INSUFFICIENT MEMORY FOR THE INPUT DATA IN COMMON BLOCK
C             /DCRT/.  THE DIMENSION MEND OF ARRAY KEND MUST BE
C             ENLARGED.  SEE THE BLOCK DATA SUBROUTINE RAYB.
C     562...  INSUFFICIENT MEMORY IN /DCRT/:
C             INSUFFICIENT MEMORY FOR THE INPUT DATA IN COMMON BLOCK
C             /DCRT/.  THE DIMENSION MSTOR OF ARRAY KSTOR MUST BE
C             ENLARGED.  SEE THE BLOCK DATA SUBROUTINE RAYB.
C     563...  INADMISSIBLE END SURFACE:
C             THE INDEX OF AN END SURFACE IN INPUT DATA IS GREATER THAN
C             THE NUMBER OF SPECIFIED SURFACES.
C     564...  INADMISSIBLE STORE SURFACE:
C             THE INDEX OF A SURFACE FOR STORING THE COMPUTED QUANTITIES
C             IS GREATER THAN THE NUMBER OF SPECIFIED SURFACES AND WHILE
C             NOT SPECIFYING THE COMPUTATIONAL VOLUME BOUNDARY, OR THE
C             INDEX OF AN AUXILIARY SURFACE FOR STORING IS NEGATIVE.
C
C DATE: 1993, DECEMBER 18
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I
C
C     (1) STRING DESCRIBING THE DATA
      READ(LUN,*) TEXTR
C     (2,3) QUANTITIES CONTROLING NUMERICAL INTEGRATION
      MODCRT=0
      READ(LUN,*) KSTORE,NEXPS,NHLF,MODCRT
      READ(LUN,*) STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT
C     (4) COORDINATE PLANES AND ISOCHRONE BOUNDING THE COMPUTATIONAL
C     VOLUME
      DO 10 I=1,6
        BOUNDR(I)=BOUNDM(I)
   10 CONTINUE
      BOUNDR(7)=999999.
      READ(LUN,*) BOUNDR
C     (5) NUMBER OF AUXILIARY SURFACES
      NSRFCA=0
      READ(LUN,*) NSRFCA
C
C     (6) END SURFACES:
C     INITIALIZING MEMORY FOR INDICES OF END SURFACES
      DO 11 I=1,MEND
        KEND(I)=0
   11 CONTINUE
C     READING INDICES OF END SURFACES:
      READ(LUN,*) (KEND(I),I=1,MEND)
      DO 12 I=1,MEND
        IF(IABS(KEND(I)).GT.NSRFC()+NSRFCA) THEN
          PAUSE 'ERROR 563 IN RAY1: INADMISSIBLE END SURFACE'
          STOP
        ELSE IF(KEND(I).EQ.0) THEN
          NEND=I-1
          GO TO 13
        END IF
   12 CONTINUE
      PAUSE 'ERROR 561 IN RAY1: INSUFFICIENT MEMORY IN /DCRT/'
      STOP
   13 CONTINUE
C
C     (7) SURFACES FOR STORING COMPUTED QUANTITIES:
C     INITIALIZING MEMORY FOR INDICES OF THE SURFACES
      DO 21 I=1,MSTOR
        KSTOR(I)=0
   21 CONTINUE
C     READING INDICES OF THE SURFACES:
      READ(LUN,*) (KSTOR(I),I=1,MSTOR)
      DO 22 I=1,MSTOR/3+1
        IF(KSTOR(I).LT.-NSRFC().OR.KSTOR(I).GT.NSRFC()+NSRFCA) THEN
          IF(KSTOR(I).LT.101.OR.KSTOR(I).GT.107) THEN
            PAUSE 'ERROR 564 IN RAY1: INADMISSIBLE STORE SURFACE'
            STOP
          END IF
        ELSE IF(KSTOR(I).EQ.0) THEN
          NSTOR=I-1
          GO TO 23
        END IF
   22 CONTINUE
      PAUSE 'ERROR 562 IN RAY1: INSUFFICIENT MEMORY IN /DCRT/'
      STOP
   23 CONTINUE
C     (8) SEQUENTIAL NUMBERS OF THE FIRST POINTS OF INTERSECTION
      DO 28 I=NSTOR+1,2*NSTOR
        KSTOR(I)=1
   28 CONTINUE
      READ(LUN,*) (KSTOR(I),I=NSTOR+1,2*NSTOR)
C     (9) SEQUENTIAL NUMBERS OF THE LAST POINTS OF INTERSECTION
      DO 29 I=2*NSTOR+1,3*NSTOR
        KSTOR(I)=999999
   29 CONTINUE
      READ(LUN,*) (KSTOR(I),I=2*NSTOR+1,3*NSTOR)
C
C     (10) AUXILIARY SURFACES:
      IF(NSRFCA.GT.0) THEN
        CALL SRFC1(LUN,NSRFCA)
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE RAY2(YL,Y,YY,IY,IEND)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12),IEND
C
C THIS SUBROUTINE CONTINUES IN THE COMPLETE TRACING OF A RAY FROM THE
C GIVEN POINT (SEE C.R.T.5.7).  FOR MODCRT.LE.2 (SEE THE INPUT DATA
C (2)), THE WHOLE RAY IS COMPUTED.  FOR MODCRT.GE.3, JUST ONE RAY
C ELEMENT IS COMPUTED.
C
C INPUT:
C     YL...   ARRAY CONTAINING LOCAL QUANTITIES AT THE GIVEN POINT.
C     Y...    ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY
C             AT THE GIVEN POINT.
C     YY...   ARRAY CONTAINING REAL AUXILIARY QUANTITIES COMPUTED ALONG
C             THE RAY AT THE GIVEN POINT.
C     IY...   ARRAY CONTAINING INTEGER AUXILIARY QUANTITIES COMPUTED
C             ALONG THE RAY AT THE GIVEN POINT.
C
C OUTPUT:
C     YL...   ARRAY CONTAINING LOCAL QUANTITIES AT THE ENDPOINT OF THE
C             RAY.
C     Y...    ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY
C             AT THE ENDPOINT OF THE RAY.
C     YY...   ARRAY CONTAINING REAL AUXILIARY QUANTITIES COMPUTED ALONG
C             THE RAY AT THE ENDPOINT OF THE RAY.
C     IY...   ARRAY CONTAINING INTEGER AUXILIARY QUANTITIES COMPUTED
C             ALONG THE RAY AT THE ENDPOINT OF THE RAY.
C     IEND... REASON OF THE TERMINATION OF THE COMPUTATION OF A RAY (SEE
C               C.R.T.5.4).  IEND=10,21,22,23,30 ARE INDICATED IN
C               SUBROUTINE CODE (FILE 'CODE.FOR').  IEND=24,25,26,32,33
C               ARE INDICATED IN SUBROUTINE TRANS (FILE 'TRANS.FOR').
C               THE REASONS IEND=40,50,60 ARE CHECKED AND INDICATED BY
C               IY(6) IN SUBROUTINE RAYCB (FILE 'RAYCB.FOR').  IEND=61
C               IS INDICATED IN THIS SUBROUTINE RAY2 USING IY(6).
C             0...  THE COMPUTATION OF THE RAY MAY CONTINUE.
C             10... RAY SATISFIES THE WHOLE CODE.
C             21... THE POINT OF INCIDENCE IS SITUATED AT A DIFFERENT
C               SURFACE THAN THAT SPECIFIED BY THE CODE.
C             22... THE NEXT ELEMENT OF THE RAY IS REQUIRED BY THE CODE
C               TO BE SITUATED IN A COMPLEX BLOCK THAT DOES NOT TOUCH
C               THE POINT OF INCIDENCE.
C             23,24... TRANSMISSION IS REQUIRED BY THE CODE AT A FREE
C               SURFACE.
C             25... RAY OF THE REQUIRED REFLECTED OR TRANSMITTED WAVE
C               IS NOT REAL-VALUED (OVERCRITICAL REFLECTION OR
C               TRANSMISSION).
C             26... S WAVE IN A LIQUID BLOCK IS REQUIRED BY THE CODE.
C             30,32... REFLECTION OR TYPE CONVERSION AT THE FICTITIOUS
C               PART OF THE INTERFACE.
C             33... ZERO REFLECTION/TRANSMISSION COEFFICIENT.
C             40... TRAVEL TIME GREATER THAN THE SPECIFIED LIMIT.
C             50... RAY HAS INTERSECTED A COORDINATE PLANE LIMITING THE
C               COMPUTATIONAL VOLUME FOR COMPLETE RAY TRACING.
C             60,61... RAY HAS INTERSECTED A SMOOTH SURFACE LIMITING THE
C               COMPUTATIONAL VOLUME FOR COMPLETE RAY TRACING.
C             OTHER VALUES OF IEND MAY INDICATE THE PROGRAM FAILURE.
C
C COMMON BLOCK /DCRT/:
      INTEGER MEND,MSTOR
      PARAMETER (MEND=128)
      PARAMETER (MSTOR=128)
      INTEGER KSTORE,NEXPS,NHLF,MODCRT
      REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7)
      INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR)
      COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH,
     *             UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL NSRFC,KOOR,CODE,RAYCB,TRANS,CONV,RPAR31,RPAR32,RPAR33
      EXTERNAL WRIT31,WRIT32,WRIT33
      INTEGER NSRFC,KOOR
C     NSRFC,BLOCK,VELOC... FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'.
C     KOOR,METRIC... FILE 'METRIC.FOR' OF THE PACKAGE 'MODEL'.
C     SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     CDE,CROSS,HIVD2,SMVPRD... FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C     HPCG... IBM SCIENTIFIC SUBROUTINE PACKAGE (FILE 'HPCG.FOR' OF THE
C             PACKAGE 'MODEL').
C     CODE... FILE 'CODE.FOR'.
C     RAYCB AND SUBSEQUENT ROUTINES... FILE 'RAYCB.FOR'.
C     TRANS,CONV AND SUBSEQUENT ROUTINES... FILE 'TRANS.FOR'.
C     RPAR31,RPAR32,RPAR33... FILE 'RPAR.FOR'.
C     WRIT31,WRIT32,WRIT33 AND SUBSEQUENT ROUTINES... FILE 'WRIT.FOR'.
C     FCT,OUTP,CDEF,PSHIFT... THIS FILE.
C
C ERROR MESSAGES:
C     570...  WRONG FUNCTION OF SUBROUTINE CODE:
C             REASON 37 OR 38 OF THE TERMINATION OF THE RAY COMPUTATION
C             IS REPORTED BY THE SUBROUTINE TRANS, WHILE IT SHOULD BE
C             REPORTED PREVIOUSLY BY THE SUBROUTINE CODE AS THE REASON
C             22 OR 23, RESPECTIVELY.
C
C DATE: 1990, NOVEMBER 18
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER KODIND,ICBNEW,I,NAMPL
      REAL YC(12)
C
C     KODIND..POSITION IN THE CODE (INDEX IN THE ARRAY KODE) CORRES-
C             PONDING TO THE NEXT ELEMENT OF THE RAY.
C     ICBNEW..THE INDEX OF THE COMPLEX BLOCK IN WHICH THE NEXT
C             ELEMENT OF THE RAY IS TO BE SITUATED, SUPPLEMENTED
C             BY THE SIGN "+" FOR P WAVE OR "-" FOR S WAVE.
C     I...    AUXILIARY LOOP VARIABLE.
C     NAMPL...NUMBER OF REAL QUANTITIES REPRESENTING COMPLEX-VALUED
C             VECTORIAL REDUCED AMPLITUDES.  SEE THE SUBROUTINE CONV OF
C             THE SUBROUTINE FILE 'TRANS.FOR'.
C     YC...   ARRAY CONTAINING REAL QUANTITIES REPRESENTING COMPLEX-
C             -VALUED VECTORIAL REDUCED AMPLITUDES.  SEE C.R.T.5.5.4 AND
C             THE SUBROUTINE CONV OF THE SUBROUTINE FILE 'TRANS.FOR'.
C
C.......................................................................
C
   10 CONTINUE
      IF(IY(6).EQ.0) THEN
C
C       (A1) ELEMENT OF THE RAY INSIDE A COMPLEX BLOCK
        CALL RAYCB(YL,Y,YY,IY)
C
C       (B1) STORAGE OF THE COMPUTED QUANTITIES AT THE INTERFACE
        CALL RPAR31(YL,Y,YY,IY)
        IF(STORE.NE.0) THEN
          CALL WRIT31(YL,Y,YY,IY)
        END IF
        DO 21 I=1,NSTOR
          IF(KSTOR(I).EQ.IY(6)) THEN
            CALL RPAR32(I,YL,Y,YY,IY)
            CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC)
            CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC)
          END IF
          IF(KSTOR(I).EQ.-IY(6).AND.KSTORE.GE.2) THEN
            IY(6)=-IY(6)
            CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC)
            CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC)
            IY(6)=-IY(6)
          END IF
   21   CONTINUE
        IF(IABS(IY(6)).LE.NSRFC()) THEN
          CALL RPAR33(YL,Y,YY,IY)
          CALL WRIT33(YL,Y,YY,IY)
        END IF
C
C       (A2) CHECK FOR THE BOUNDARY OF THE COMPUTATIONAL VOLUME
        IF(IABS(IY(6)).GE.107) THEN
          IEND=40
          GO TO 99
        ELSE IF(IABS(IY(6)).GE.101) THEN
          IEND=50
          GO TO 99
        ELSE IF(IABS(IY(6)).GT.NSRFC()) THEN
          IEND=60
          GO TO 99
        ELSE IF(MODCRT.GE.3) THEN
          IEND=0
          GO TO 99
        END IF
C
      ELSE
C
C       (B2) INTERPRETATION OF THE ELEMENTARY WAVES CODE
        CALL CODE(IY,KODIND,ICBNEW,IEND)
        IF(IEND.NE.0) THEN
          GO TO 99
        END IF
C
C       (B3) CHECK FOR THE END SURFACES
C       NOTE: REFLECTION FROM AN END SURFACE IS ALLOWED.
        IF(IABS(IY(5)).NE.IABS(ICBNEW)) THEN
          DO 23 I=1,NEND
            IF(IABS(KEND(I)).EQ.IABS(IY(6))) THEN
              IEND=61
              GO TO 99
            END IF
   23     CONTINUE
        END IF
C
C       (B4) TRANSFORMATION ACROSS THE INTERFACE
        CALL TRANS(YL,Y,YY,IY,KODIND,ICBNEW,IEND)
        IF(IEND.GE.37) THEN
          PAUSE 'ERROR 570 IN RAY2: WRONG FUNCTION OF SUBROUTINE CODE'
          STOP
        ELSE IF(IEND.NE.0) THEN
          GO TO 99
        END IF
C
C       (B5) STORAGE OF THE COMPUTED QUANTITIES AT THE INTERFACE
        CALL RPAR31(YL,Y,YY,IY)
        IF(STORE.NE.0) THEN
          CALL WRIT31(YL,Y,YY,IY)
        END IF
        DO 25 I=1,NSTOR
          IF(KSTOR(I).EQ.IY(6)) THEN
            CALL RPAR32(I,YL,Y,YY,IY)
            IF(KSTORE.LE.1) THEN
              CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC)
              CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC)
            END IF
          END IF
   25   CONTINUE
C
C       (B6) THE RAY MAY CONTINUE IN THE NEXT COMPLEX BLOCK
        IY(6)=0
C
      END IF
      GO TO 10
C
   99 CONTINUE
      RETURN
      END
C
C=======================================================================
C
