C SUBROUTINE FILE 'AP.FOR': APPLICATIONS AND PROCESSING OF THE RESULTS
C OF COMPLETE RAY TRACING.
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C THIS FILE CONTAINS SUBROUTINES RELATED TO THE CHAPTER 7 OF THE PAPER
C ON C.R.T., DESIGNED TO READ FROM THE FILES THE QUANTITIES DESCRIBING
C THE POINTS OF RAYS, AND TO EVALUATE MANY OTHER QUANTITIES USED IN
C SEISMOLOGY AND DISCUSSED IN THE CHAPTER C.R.T.7.  THESE SUBROUTINES
C MAY BE INLUDED IN USER'S APPLICATION PROGRAMS FOLLOWING THE COMPLETE
C RAY TRACING PROGRAM.  INDIVIDUAL SUBROUTINES CORRESPOND TO THE
C INDIVIDUAL SECTIONS OF THE CHAPTER C.R.T.7, AND MAY CALL OTHER
C SUBROUTINES OF THE FILES OF THE PACKAGES 'MODEL' AND 'CRT' COMPOSING
C THE COMPLETE RAY TRACING PROGRAM.
C
C ATTENTION:  THE LINES DEALING WITH CURVILINEAR COORDINATES ARE DENOTED
C BY '*' IN THE FIRST COLUMN.  THAT IS WHY THIS CODE WORKS PROPERLY ONLY
C IN CARTESIAN COORDINATES.  TO CONSIDER CURVILINEAR COORDINATES, EACH
C '*' IN THE FIRST COLUMN HAS TO BE REPLACED BY A SPACE.  IN SUCH A
C CASE, SUBROUTINE METR1 HAS TO BE CALLED FIRST TO SPECIFY THE
C COORDINATE SYSTEM.
C
C THIS IS A PRELIMINARY VERSION, CONTAINING ONLY SOME OF THE ROUTINES
C CORRESPONDING TO SECTIONS OF THE CHAPTER C.R.T.7.  THE ROUTINES
C DENOTED BY * IN THE SECOND COLUMN WILL LIKELY BE CODED IN THE NEAR
C FUTURE, OTHERS LATER ON.
C
C THIS FILE CONSISTS OF THE FOLLOWING EXTERNAL PROCEDURES:
C     POINTB..BLOCK DATA SUBROUTINE DEFINING COMMON BLOCK /POINTC/ TO
C             STORE THE QUANTITIES DESCRIBING A POINT ON A RAY.
C     AP00... SUBROUTINE DESIGNED TO READ FROM THE FILES THE QUANTITIES
C             DESCRIBING A POINT ON A RAY, AND TO STORE THEM INTO THE
C             COMMON BLOCK /POINTC/.  THE INPUT FILES ARE ASSUMED TO
C             HAVE THE STRUCTURE OF THE OUTPUT FILES OF THE COMPLETE RAY
C             TRACING PROGRAM 'CRT', WRITTEN BY THE SUBROUTINES OF THE
C             FILE 'WRIT.FOR'.
C     AP01... SUBROUTINE DESIGNED TO EVALUATE THE TRAVEL TIME FROM THE
C             INITIAL POINT OF A RAY TO A POINT SITUATED ON A RAY, SEE
C             C.R.T.7.1.
C     AP02... SUBROUTINE DESIGNED TO EVALUATE THE COMPONENTS OF THE
C             SLOWNESS VECTOR EITHER AT THE INITIAL POINT OF A RAY OR AT
C             A POINT SITUATED ON A RAY, SEE C.R.T.7.2.
C     AP03... SUBROUTINE DESIGNED TO EVALUATE THE COVARIANT COMPONENTS
C             OF THE BASIS VECTORS OF THE RAY-CENTRED COORDINATE SYSTEM
C             AT THE INITIAL POINT OF A RAY OR AT A POINT SITUATED ON A
C             RAY, SEE C.R.T.7.3.
C     AP03A...AUXILIARY SUBROUTINE TO AP03, EVALUATING THE BASIS OF THE
C             INTRINSIC RAY-CENTRED COORDINATE SYSTEM AT THE GIVEN
C             POINT.
C*    AP04
C     AP05... SUBROUTINE DESIGNED TO EVALUATE THE COMPONENTS OF THE
C             MATRIX OF GEOMETRICAL SPREADING AT A POINT SITUATED ON A
C             RAY, SEE C.R.T.7.5.
C     AP06... SUBROUTINE DESIGNED TO EVALUATE THE COMPONENTS OF THE
C             TRANSFORMATION MATRIX P AT A POINT SITUATED ON A RAY, SEE
C             C.R.T.7.6.
C     AP07... SUBROUTINE DESIGNED TO EVALUATE THE GEOMETRICAL SPREADING
C             AT A POINT SITUATED ON A RAY, SEE C.R.T.7.7.
C     AP08... SUBROUTINE DESIGNED TO EVALUATE THE COMPONENTS OF THE
C             SYMMETRIC 3*3 MATRICES M AND N OF SECOND DERIVATIVES OF
C             THE TRAVEL-TIME FIELD AT A POINT SITUATED ON A RAY, SEE
C             C.R.T.7.8.  SUBROUTINE AP03 SHOULD BE CALLED BEFORE THE
C             INVOCATION OF AP08 TO DEFINE THE BASIS OF R.C.C.S.
C*    AP09
C*    AP10(XI,X)
C*    AP11(XI,PI)
C*    AP12(XI,X)
C*    AP13(XI,XF,X)
C*    AP14
C     AP15... SUBROUTINE DESIGNED TO EVALUATE THE RAY AMPLITUDES AT A
C             POINT SITUATED ON A RAY, SEE C.R.T.7.15.  SUBROUTINE AP03
C             SHOULD BE CALLED BEFORE THE INVOCATION OF AP15 TO DEFINE
C             THE BASIS OF R.C.C.S.
C*    AP16
C
C NOTE: THERE ARE NO APPLICATION ROUTINES CORRESPONDING TO THE SECTIONS
C 7.18, 7.23 AND 7.27 OF THE PAPER ON C.R.T.  AP28 AND SUBSEQUENT
C APPLICATIONS DO NOT CORRESPOND TO ANY SECTION OF THE PAPER ON C.R.T.
C
C STORAGE IN THE MEMORY:
C     WHEN PROCESSING OF THE RESULTS OF COMPLETE RAY TRACING, THE
C     QUANTITIES DESCRIBING SOME POINTS ON A RAY ARE REQUIRED TO BE
C     KNOWN.  IN THE CHAPTER 7 OF THE PAPER ON C.R.T., THREE DIFFERENT
C     POINTS SITUATED ON A RAY ARE INTRODUCED:
C     O/O (O SUBSCRIPT O)... INITIAL POINT OF THE RAY.
C     O/S (O SUBSCRIPT S)... ANOTHER POINT SITUATED ON THE RAY.  IN SOME
C             APPLICATIONS IT MAY BE TREATED AS THE ENDPOINT OF THE RAY.
C     O/F (O SUBSCRIPT F)... ANOTHER POINT SITUATED ON THE RAY, USUALLY
C             SITUATED BETWEEN THE POINTS O/O AND O/S, SEE C.R.T.7.13:
C             FRESNEL VOLUMES.  THIS POINT IS REQUIRED JUST BY FEW
C             APPLICATIONS AND USUALLY NEED NOT BE DEFINED.
C     THE QUANTITIES DESCRIBING THE POINTS O/O, O/S (AND, POSSIBLY, O/F)
C     OF A RAY ARE STORED BY THE SUBROUTINE AP00 INTO THE COMMON BLOCK
C     /POINTC/ DEFINED IN THE FOLLOWING SUBROUTINE:
C     ------------------------------------------------------------------
      BLOCK DATA POINTB
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
      SAVE /POINTC/
      DATA IWAVE/0/,IRAY/0/
      END
C     ------------------------------------------------------------------
C     IWAVE...INDEX OF THE ELEMENTARY WAVE (OUTPUT OF THE SUBROUTINE
C             CODE1).  ELEMENTARY WAVES ARE INDEXED BY 1,2,3,... .
C             IWAVE=0 IF THE END OF THE INPUT FILE IS ENCOUNTERED.
C     IRAY... INDEX OF THE RAY (OUTPUT OF THE SUBROUTINE RPAR2).  RAYS
C             WITHIN EACH ELEMENTARY WAVE ARE INDEXED BY 1,2,3,... .
C             NOTE THAT SOME OF THE INDEXED RAYS NEED NOT EXIST.
C             IRAY=0 IF THE END OF THE INPUT FILE IS ENCOUNTERED.
C     IPT...  INDEX OF THE POINT ON A RAY.  THE POINTS O/S ARE INDEXED
C             SEQUENTIALLY ALONG EACH RAY, TAKING THE VALUES 1,2,3,...
C             IF THE POINT IS THE FIRST CONSIDERED POINT OF A NEW
C             ELEMENTARY WAVE, IPT=0 INSTEAD OF 1.
C     NYF,ICB1F,ISRFF,YLF,YF... QUANTITIES AT THE POINT O/F.  THEY NEED
C             NOT BE DEFINED (OPTION IN THE SUBROUTINE AP00).  THEIR
C             MANING IS THE SAME AS OF NY,ICB1,ISRF,YL,Y BELOW.
C             NYF=0... INDICATION THAT THESE QUANTITIES ARE NOT DEFINED,
C               I.E. THAT THE POINT O/F HAS NOT BEEN DEFINED.
C     NY,ICB1,ISRF,YL,Y... QUANTITIES AT THE POINT O/S.
C             THE DETAILED DESCRIPTION OF THE QUANTITIES YL AND Y MAY BE
C             FOUND IN THE FILE 'RAY.FOR'.  A BRIEF DESCRIPTION FOLLOWS:
C     NY=IY(1)=27+NAMPL... NUMBER OF THE BASIC QUANTITIES DESCRIBING THE
C             POINT OF THE RAY, SEE THE FILE 'RAY.FOR'.
C     ICB1=IY(5)... INDEX OF THE COMPLEX BLOCK IN WHICH THE STORED POINT
C             OF THE RAY IS SITUATED, SUPPLEMENTED BY A SIGN '+' FOR P
C             WAVE AND SIGN '-' FOR S WAVE.
C     ISRF=IY(6)... 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     YL...   ARRAY CONTAINING LOCAL QUANTITIES AT THE POINT OF THE RAY
C             (C.R.T.5.5.4).
C     Y...    ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY
C             (C.R.T.5.2.1).
C     ICB1I,YLI,YI... QUANTITIES AT THE INITIAL POINT O/I OF THE RAY:
C     ICB1I...INDEX OF THE COMPLEX BLOCK IN WHICH THE INITIAL POINT OF
C             THE RAY IS SITUATED, SUPPLEMENTED BY A SIGN '+' FOR P WAVE
C             AND SIGN '-' FOR S WAVE (SEE C.R.T.6.1).
C     IEND... REASON OF THE TERMINATION OF THE COMPUTATION OF A RAY (SEE
C             C.R.T.5.4).  FOR A DETAILED DESCRIPTION SEE SUBROUTINES
C             RAY2 (FILE 'RAY.FOR') AND INIT2 (FILE 'INIT.FOR').
C     ISHEET..RAY-HISTORY INDEX.  THE DIFFERENT RAY HISTORIES ARE
C             CONSECUTIVELY INDEXED BY POSITIVE INTEGERS 1,2,3,...
C             ACCORDING TO THEIR APPEARANCE DURING RAY TRACING.
C             THE RAY HISTORIES ARE INDEXED INDEPENDENTLY WITHIN EACH
C             ELEMENTARY WAVE.  THE INDICES ARE THE OUTPUT OF SUBROUTINE
C             RPAR4 OF THE FILE 'RPAR.FOR'.
C             THE RAY-HISTORY INDICES ARE COMPLEMENTED WITH SIGN:
C             POSITIVE - SUCCESSFUL RAY (CROSSING REFERENCE SURFACE),
C             NEGATIVE - UNSUCCESSFUL RAY (TERMINATING BEFORE CROSSING
C             REFERENCE SURFACE).
C     YLI...  ARRAY CONTAINING THE VALUES OF THE QUANTITIES YL(1)-YL(6),
C             SEE C.R.T.5.5.4, DESCRIBING THE LOCAL PROPERTIES OF THE
C             MODEL AT THE INITIAL POINT OF THE RAY.  SEE THE LIST OF
C             YL(1) TO YL(6) IN THE FILE 'RAY.FOR'.
C     YI...   ARRAY CONTAINING THE QUANTITIES DESCRIBING THE PROPERTIES
C             OF THE RAYS AND OF THE TRAVEL-TIME FIELD AT THE INITIAL
C             POINT OF THE RAY, SEE C.R.T.6.1.  THESE QUANTITIES ARE
C             LISTED IN THE FILE 'INIT.FOR'.
C     THE STORAGE LOCATIONS OF THE COMMON BLOCK /POINTC/ ARE REDEFINED
C     IN THE EXTERNAL PROCEDURE AP00.  COMMON BLOCK /POINTC/ IS INCLUDED
C     IN THE APPLICATION SUBROUTINES OF THIS FILE AND MAY ALSO BE
C     INCLUDED IN ANY USER'S SUBROUTINE PROCESSING THE RESULTS OF THE
C     COMPLETE RAY TRACING PROGRAM.  THE QUANTITIES DESCRIBING ANOTHER
C     POINT ALONG A RAY CAN BE STORED INTO THE COMMON BLOCK /POINTC/ BY
C     ANOTHER INVOCATION OF THE SUBROUTINE AP00.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C=======================================================================
C
      SUBROUTINE AP00(LU1,LU2,LU3)
      INTEGER LU1,LU2,LU3
C
C THIS SUBROUTINE READS FROM THE GIVEN FILES THE QUANTITIES DESCRIBING
C SOME POINTS ON A RAY, AND STORES THEM INTO THE COMMON BLOCK /POINTC/.
C THE INPUT FILES ARE ASSUMED TO HAVE THE STRUCTURE OF THE OUTPUT FILES
C OF THE COMPLETE RAY TRACING PROGRAM 'CRT', WRITTEN BY THE SUBROUTINES
C OF THE FILE 'WRIT.FOR'.  WHEN ALL POINTS OF THE GIVEN FILES ARE READ
C OVER, IWAVE AND IRAY OF THE COMMON BLOCK /POINTC/ ARE SET TO ZEROS.
C THE LOCATIONS IWAVE AND IRAY OF THE COMMON BLOCK /POINTC/ SHOULD NOT
C BE CHANGED BUT WITH THE FOLLOWING EXEPTION: THEY MUST BE SET TO ZEROS
C BEFORE THIS SUBROUTINE IS CALLED WITH ALTERED (REOPENED) FILES
C CORRESPONDING TO THE INPUT PARAMETERS LU1, LU2 AND LU3.
C IN THE CHAPTER 7 OF THE PAPER ON C.R.T., THE INITIAL POINT OF THE RAY
C IS DENOTED O/O (O SUBSCRIPT O) AND THE POINTS ON THE RAY ARE DENOTED
C O/F AND O/S.  AFTER THE INVOCATION OF THIS SUBROUTINE, A USER MAY WANT
C TO CALL HIS OWN SUBROUTINE DECIDING WHETHER THE POINT ON A RAY IS TO
C BE IGNORED (E.G. WHEN IT IS TOO FAR FROM THE RECEIVERS) OR PROCESSED
C ON.
C
C INPUT:
C     LU1...  ZERO IF THE POINT O/F OF THE RAY NEED NOT BE DEFINED.
C             OTHERWISE THE LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT
C             DEVICE CONTAINING A FILE WITH THE QUANTITIES ALONG RAYS
C             (SEE C.R.T.5.5.1) OR A FILE WITH THE QUANTITIES AT A
C             SPECIFIED SURFACE (SEE C.R.T.5.5.2).  THE POINTS O/F WILL
C             BE READ FROM THIS FILE.  IF THERE WILL BE A POINT O/S OF
C             THE SAME RAY IN THE FILE LU2, THE QUANTITIES CORRESPONDING
C             TO THESE POINTS WILL BE STORED INTO THE COMMON BLOCK
C             /POINTC/.
C     LU2...  LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE
C             CONTAINING A FILE WITH THE QUANTITIES ALONG RAYS (SEE
C             C.R.T.5.5.1, JUST FOR LU1=0) OR A FILE WITH THE QUANTITIES
C             AT A SPECIFIED SURFACE (SEE C.R.T.5.5.2).  FOR LU1=0, ALL
C             POINTS O/S FROM THIS FILE WILL BE SUCCESSIVELY STORED INTO
C             THE COMMON BLOCK /POINTC/.
C     LU3...  LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE
C             CONTAINING THE FILE WITH THE QUANTITIES AT THE INITIAL
C             POINTS OF RAYS, CORRESPONDING TO THE ABOVE FILE LU1 OR LU2
C             (SEE C.R.T.6.1).
C THE INPUT PARAMETERS ARE NOT ALTERED.
C
C NO OUTPUT.
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C ALL THE STORAGE LOCATIONS OF THE COMMON BLOCKS ARE DEFINED IN THIS
C SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL POINTB
C     POINTB..BLOCK DATA SUBROUTINE OF THIS FILE.
C
C ERROR MESSAGES:
C     701...  WRONG FILE WITH COMPUTED POINTS:
C             INDEX OF THE ELEMENTARY WAVE IS NOT POSITIVE.  MAYBE, THE
C             FILE LU1 HAS THE STRUCTURE OF A FILE LU3.
C     702...  WRONG FILE WITH COMPUTED POINTS:
C             INDEX OF THE ELEMENTARY WAVE IS NOT POSITIVE.  MAYBE, THE
C             FILE LU2 HAS THE STRUCTURE OF A FILE LU3.
C     703...  WRONG FILE WITH INITIAL POINTS:
C             INDEX OF THE ELEMENTARY WAVE IS NOT SUPPLIED WITH A MINUS
C             SIGN. MAYBE, THE FILE LU3 HAS THE STRUCTURE OF A FILE LU2.
C     704...  THE WAVE NOT FOUND:
C             THE INITIAL POINT OF A RAY FROM FILE LU2 IS NOT FOUND IN
C             THE FILE LU3.
C     705...  INITIAL POINT OF THE RAY NOT FOUND:
C             THE INITIAL POINT OF A RAY FROM FILE LU2 IS NOT FOUND IN
C             THE FILE LU3.
C     706...  END OF THE FILE WITH INITIAL POINTS:
C             THE INITIAL POINT OF A RAY FROM FILE LU2 IS NOT FOUND IN
C             THE FILE LU3.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER IWAVEF,IRAYF,IWAVEI,IRAYI,I
C     IWAVEF,IRAYF... VALUES OF IWAVE,IRAY CORRESPONDING TO THE LAST
C             POINT O/F OF A RAY READ IN DURING THE CURRENT INVOCATION.
C     IWAVEI,IRAYI... VALUES OF IWAVE,IRAY CORRESPONDING TO THE LAST
C             READ IN INITIAL POINT OF A RAY.
C     I...    AUXILIARY LOOP VARIABLE.
C
C     READING THE RESULTS OF THE COMPLETE RAY TRACING:
      IWAVEF=-1
      IRAYF=0
      NYF=0
      IWAVEI=IWAVE
      IRAYI=IRAY
C
C     POINTS O/F AND O/S ON A RAY:
C     FOR LU1=0 JUST THE POINT O/S IS BEING READ,
C     OTHERWISE, POINTS O/F AND O/S MUST BE SITUATED ON THE SAME RAY.
   10 CONTINUE
      IF(LU1.NE.0.AND.(IWAVEF.LT.IWAVE.OR.
     *                (IWAVEF.EQ.IWAVE.AND.IRAYF.LT.IRAY))) THEN
C       NEW POINT O/F ON A RAY:
        READ(LU1,END=80)
     *               IWAVEF,IRAYF,NYF,ICB1F,ISRFF,XF,YLF,(YF(I),I=1,NYF)
        IF (IWAVEF.LE.0) THEN
          PAUSE 'ERROR 701 IN AP00: WRONG FILE WITH COMPUTED POINTS'
          STOP
        END IF
        GO TO 10
      END IF
      IF(LU1.EQ.0.OR.IWAVEF.GT.IWAVE.OR.
     *              (IWAVEF.EQ.IWAVE.AND.IRAYF.GT.IRAY)) THEN
C       NEW POINT O/S ON A RAY:
        READ(LU2,END=80) IWAVE,IRAY,NY,ICB1,ISRF,X,YL,(Y(I),I=1,NY)
        IF (IWAVE.LE.0) THEN
          PAUSE 'ERROR 702 IN AP00: WRONG FILE WITH COMPUTED POINTS'
          STOP
        END IF
        IF(LU1.NE.0) THEN
          GO TO 10
        END IF
      END IF
C
C     DEFINING IPT:
      IF(IWAVE.NE.IWAVEI) THEN
        IPT=0
      ELSE IF(IRAY.NE.IRAYI) THEN
        IPT=1
      ELSE
        IPT=MAX0(1,IPT)+1
      END IF
C
C     INITIAL POINT O/O OF A RAY:
   20 CONTINUE
      IF(IWAVE.NE.IWAVEI.OR.IRAY.NE.IRAYI) THEN
        IF(IWAVE.LT.IWAVEI) THEN
          WRITE(*,'('' WAVE:'',I4,'',  RAY:'',I6)') IWAVE,IRAY
          PAUSE 'ERROR 704 IN AP00: THE WAVE NOT FOUND'
          STOP
        ELSE IF(IWAVE.EQ.IWAVEI.AND.IRAY.LT.IRAYI) THEN
          WRITE(*,'('' WAVE:'',I4,'',  RAY:'',I6)') IWAVE,IRAY
          PAUSE 'ERROR 705 IN AP00: INITIAL POINT OF THE RAY NOT FOUND'
          STOP
        ELSE
C         INITIAL POINT O/O OF A RAY
          READ(LU3,END=90) IWAVEI,IRAYI,ICB1I,IEND,ISHEET,YLI,YI
          IF(IWAVEI.LT.0) THEN
            IWAVEI=-IWAVEI
          ELSE
            PAUSE 'ERROR 703 IN AP00: WRONG FILE WITH INITIAL POINTS'
            STOP
          END IF
          GO TO 20
        END IF
      END IF
      RETURN
C
C     END OF THE FILE WITH THE COMPUTED POINTS:
   80 CONTINUE
      IWAVE=0
      IRAY=0
      IPT=0
      RETURN
C
C     END OF THE FILE WITH THE INITIAL POINTS OF RAYS:
   90 CONTINUE
      PAUSE 'ERROR 706 IN AP00: END OF THE FILE WITH INITIAL POINTS'
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE AP01(TT,TTIM)
      REAL TT,TTIM
C
C THIS SUBROUTINE EVALUATES THE TRAVEL TIME FROM THE INITIAL POINT OF A
C RAY TO A POINT SITUATED ON A RAY, SEE C.R.T.7.1.
C
C NO INPUT.
C
C OUTPUT:
C     TT...   TRAVEL TIME FROM THE INITIAL POINT O/O OF A RAY TO THE
C             POINT O/S READ INTO THE COMMON BLOCK /POINTC/ BY THE LAST
C             INVOCATION OF THE SUBROUTINE AP00.
C     TTIM... IMAGINARY TRAVEL TIME FROM THE INITIAL POINT O/O OF A RAY
C             TO THE POINT O/S (SEE C.R.T.7.1).
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
      TT  =Y(1)-YI(1)
      TTIM=Y(2)-YI(2)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP02(PI,P)
      REAL PI(6),P(6)
C
C THIS SUBROUTINE EVALUATES THE COMPONENTS OF THE SLOWNESS VECTOR EITHER
C AT THE INITIAL POINT OF A RAY OR AT A POINT SITUATED ON A RAY, SEE
C C.R.T.7.2.
C
C NO INPUT.
C
C OUTPUT:
C     PI...   THREE COVARIANT (PI(1:3)) AND THREE CONTRAVARIANT
C             (PI(4:6)) COMPONENTS OF THE GRADIENT OF THE TRAVEL TIME
C             FIELD AT THE INITIAL POINT O/O OF A RAY.
C     P...    THREE COVARIANT (P(1:3)) AND THREE CONTRAVARIANT (P(4:6))
C             COMPONENTS OF THE GRADIENT OF THE TRAVEL TIME FIELD AT THE
C             POINT O/S READ INTO THE COMMON BLOCK /POINTC/ BY THE LAST
C             INVOCATION OF THE SUBROUTINE AP00.
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
*     INTEGER KOOR
*     EXTERNAL KOOR,METRIC,SMVPRD
C     KOOR,METRIC... FILE 'METRIC.FOR' OF THE PACKAGE 'MODEL'.
C     SMVPRD... FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR THE METRIC TENSOR ETC.:
C
*     REAL G(12),GAMMA(18),GSQRD
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /APCAUX/ BUT TO SHARE THE MEMORY.
*     COMMON/APCAUX/ G,GAMMA,GSQRD
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C
C.......................................................................
C
*     REAL PI4,PI5,PI6
*     SAVE PI4,PI5,PI6
*     INTEGER JWAVE,JRAY
*     SAVE    JWAVE,JRAY
*     DATA    JWAVE,JRAY/0,0/
C
C.......................................................................
C
C     COVARIANT COMPONENTS:
      PI(1)=YI(6)
      PI(2)=YI(7)
      PI(3)=YI(8)
      P (1)=Y (6)
      P (2)=Y (7)
      P (3)=Y (8)
C
C     CONTRAVARIANT COMPONENTS:
*     IF(KOOR().NE.0) THEN
C       CURVILINEAR COORDINATES:
*       IF(IWAVE.NE.JWAVE.OR.IRAY.NE.JRAY) THEN
*         CALL METRIC(YI(3),GSQRD,G,GAMMA)
*         CALL SMVPRD(G(7),YI(6),YI(7),YI(8),PI4,PI5,PI6)
*         JWAVE=IWAVE
*         JRAY=IRAY
*       END IF
*       PI(4)=PI4
*       PI(5)=PI5
*       PI(6)=PI6
*       CALL METRIC(Y (3),GSQRD,G,GAMMA)
*       CALL SMVPRD(G(7),Y(6),Y(7),Y(8),P(4),P(5),P(6))
*     ELSE
C       CARTESIAN COORDINATES:
        PI(4)=YI(6)
        PI(5)=YI(7)
        PI(6)=YI(8)
        P (4)=Y (6)
        P (5)=Y (7)
        P (6)=Y (8)
*     END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP03(IUSER,HI,H,HUI)
      INTEGER IUSER
      REAL HI(18),H(18),HUI(9)
C
C THIS SUBROUTINE EVALUATES THE COVARIANT COMPONENTS OF THE BASIS
C VECTORS OF THE RAY-CENTRED COORDINATE SYSTEM AT THE INITIAL POINT OF A
C RAY OR AT A POINT SITUATED ON A RAY, SEE C.R.T.7.3.
C
C INPUT:
C     IUSER...IUSER=0... INTRINSIC CHOICE OF POLARIZATION VECTORS.
C               ANY OTHER INPUT NEED NOT BE SPECIFIED.
C             IUSER=1... USER'S CHOICE OF POLARIZATION VECTORS AT THE
C               INITIAL POINT O/O OF THE RAY.
C               HI(1:9) MUST BE SPECIFIED.
C               IF HUI(1:9) HAS ALREADY BEEN EVALUATED FOR THE RAY, IT
C               MUST BE SPECIFIED ON THE INPUT TOO.
C             IUSER=2... USER'S CHOICE OF POLARIZATION VECTORS AT THE
C               INITIAL POINT O/O OF THE RAY.
C               HI(10:18) MUST BE SPECIFIED.
C               IF HUI(1:9) HAS ALREADY BEEN EVALUATED FOR THE RAY, IT
C               MUST BE SPECIFIED ON THE INPUT TOO.
C             IUSER=3... TRANSFORMATION MATRIX FROM THE INTRINSIC
C               RAY-CENTRED TO THE USER'S RAY-CENTRED COORDINATE SYSTEM
C               IS GIVEN.
C               HI(1:9) NEED NOT BE SPECIFIED.
C               HUI(1:9) HAS TO BE SPECIFIED.
C     HI(1:9)... COVARIANT COMPONENTS OF THE BASIS VECTORS OF THE USER'S
C             RAY-CENTRED COORDINATE SYSTEM AT THE INITIAL POINT O/O OF
C             A RAY.
C     HUI...  COMPONENTS OF THE 3*3 TRANSFORMATION MATRIX FROM THE
C             INTRINSIC TO THE USER'S RAY-CENTRED COORDINATE SYSTEM.
C
C OUTPUT:
C     HI...   COVARIANT (HI(1:9)) AND CONTRAVARIANT (HI(10:18))
C             COMPONENTS OF THE BASIS VECTORS OF THE USER'S RAY-CENTRED
C             COORDINATE SYSTEM AT THE INITIAL POINT O/O OF A RAY.
C             IF HI(1:18) HAS ALREADY BEEN EVALUATED FOR THE RAY, JUST
C             THE COPY OF INPUT VALUES.
C     H...    COVARIANT (H(1:9)) AND CONTRAVARIANT (H(10:18)) COMPONENTS
C             OF THE BASIS VECTORS OF THE USER'S RAY-CENTRED COORDINATE
C             SYSTEM AT THE POINT O/S READ INTO THE COMMON BLOCK
C             /POINTC/ BY THE LAST INVOCATION OF THE SUBROUTINE AP00.
C     HUI...  COMPONENTS OF THE 3*3 TRANSFORMATION MATRIX FROM THE
C             INTRINSIC TO THE USER'S RAY-CENTRED COORDINATE SYSTEM.
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL AP03A
C     KOOR,METRIC... FILE 'METRIC.FOR' OF THE PACKAGE 'MODEL'.
C     SMVPRD..FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C     AP03A...THIS FILE.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
      REAL HA(18),HA01,HA02,HA03,HA04,HA05,HA06,HA07,HA08,HA09
      REAL        HA10,HA11,HA12,HA13,HA14,HA15,HA16,HA17,HA18
      EQUIVALENCE (HA(01),HA01),(HA(02),HA02),(HA(03),HA03)
      EQUIVALENCE (HA(04),HA04),(HA(05),HA05),(HA(06),HA06)
      EQUIVALENCE (HA(07),HA07),(HA(08),HA08),(HA(09),HA09)
      EQUIVALENCE (HA(10),HA10),(HA(11),HA11),(HA(12),HA12)
      EQUIVALENCE (HA(13),HA13),(HA(14),HA14),(HA(15),HA15)
      EQUIVALENCE (HA(16),HA16),(HA(17),HA17),(HA(18),HA18)
C
C     HA...   AUXILIARY STORAGE LOCATION FOR COVARIANT (HA(1:9)) AND
C             CONTRAVARIANT (HA(10:18)) COMPONENTS OF THE BASIS VECTORS
C             OF THE INTRINSIC RAY-CENTRED COORDINATE SYSTEM.
C
      INTEGER JWAVE,JRAY
      SAVE    JWAVE,JRAY
      DATA    JWAVE,JRAY/0,0/
C
C.......................................................................
C
      IF(IWAVE.NE.JWAVE.OR.IRAY.NE.JRAY) THEN
        IF(IUSER.EQ.0) THEN
          CALL AP03A(YI,HI)
          HUI(1)=1.
          HUI(2)=0.
          HUI(3)=0.
          HUI(4)=0.
          HUI(5)=1.
          HUI(6)=0.
          HUI(7)=0.
          HUI(8)=0.
          HUI(9)=1.
        ELSE
          CALL AP03A(YI,HA)
          IF(IUSER.NE.2) THEN
            IF(IUSER.EQ.1) THEN
              HUI(1)=HI(1)*HA10+HI(2)*HA11+HI(3)*HA12
              HUI(2)=HI(4)*HA10+HI(5)*HA11+HI(6)*HA12
              HUI(3)=HI(7)*HA10+HI(8)*HA11+HI(9)*HA12
              HUI(4)=HI(1)*HA13+HI(2)*HA14+HI(3)*HA15
              HUI(5)=HI(4)*HA13+HI(5)*HA14+HI(6)*HA15
              HUI(6)=HI(7)*HA13+HI(8)*HA14+HI(9)*HA15
              HUI(7)=HI(1)*HA16+HI(2)*HA17+HI(3)*HA18
              HUI(8)=HI(4)*HA16+HI(5)*HA17+HI(6)*HA18
              HUI(9)=HI(7)*HA16+HI(8)*HA17+HI(9)*HA18
            END IF
            HI(10)=HUI(1)*HA10+HUI(4)*HA11+HUI(7)*HA12
            HI(11)=HUI(2)*HA10+HUI(5)*HA11+HUI(8)*HA12
            HI(12)=HUI(3)*HA10+HUI(6)*HA11+HUI(9)*HA12
            HI(13)=HUI(1)*HA13+HUI(4)*HA14+HUI(7)*HA15
            HI(14)=HUI(2)*HA13+HUI(5)*HA14+HUI(8)*HA15
            HI(15)=HUI(3)*HA13+HUI(6)*HA14+HUI(9)*HA15
            HI(16)=HUI(1)*HA16+HUI(4)*HA17+HUI(7)*HA18
            HI(17)=HUI(2)*HA16+HUI(5)*HA17+HUI(8)*HA18
            HI(18)=HUI(3)*HA16+HUI(6)*HA17+HUI(9)*HA18
          END IF
          IF(IUSER.NE.1) THEN
            IF(IUSER.EQ.2) THEN
              HUI(1)=HI(10)*HA01+HI(11)*HA02+HI(12)*HA03
              HUI(2)=HI(13)*HA01+HI(14)*HA02+HI(15)*HA03
              HUI(3)=HI(16)*HA01+HI(17)*HA02+HI(18)*HA03
              HUI(4)=HI(10)*HA04+HI(11)*HA05+HI(12)*HA06
              HUI(5)=HI(13)*HA04+HI(14)*HA05+HI(15)*HA06
              HUI(6)=HI(16)*HA04+HI(17)*HA05+HI(18)*HA06
              HUI(7)=HI(10)*HA07+HI(11)*HA08+HI(12)*HA09
              HUI(8)=HI(13)*HA07+HI(14)*HA08+HI(15)*HA09
              HUI(9)=HI(16)*HA07+HI(17)*HA08+HI(18)*HA09
            END IF
            HI( 1)=HUI(1)*HA01+HUI(4)*HA02+HUI(7)*HA03
            HI( 2)=HUI(2)*HA01+HUI(5)*HA02+HUI(8)*HA03
            HI( 3)=HUI(3)*HA01+HUI(6)*HA02+HUI(9)*HA03
            HI( 4)=HUI(1)*HA04+HUI(4)*HA05+HUI(7)*HA06
            HI( 5)=HUI(2)*HA04+HUI(5)*HA05+HUI(8)*HA06
            HI( 6)=HUI(3)*HA04+HUI(6)*HA05+HUI(9)*HA06
            HI( 7)=HUI(1)*HA07+HUI(4)*HA08+HUI(7)*HA09
            HI( 8)=HUI(2)*HA07+HUI(5)*HA08+HUI(8)*HA09
            HI( 9)=HUI(3)*HA07+HUI(6)*HA08+HUI(9)*HA09
          END IF
        END IF
        JWAVE=IWAVE
        JRAY=IRAY
      END IF
C
      IF(IUSER.EQ.0) THEN
        CALL AP03A(Y,H)
      ELSE
        CALL AP03A(Y,HA)
        H( 1)=HUI(1)*HA01+HUI(4)*HA02+HUI(7)*HA03
        H( 2)=HUI(2)*HA01+HUI(5)*HA02+HUI(8)*HA03
        H( 3)=HUI(3)*HA01+HUI(6)*HA02+HUI(9)*HA03
        H( 4)=HUI(1)*HA04+HUI(4)*HA05+HUI(7)*HA06
        H( 5)=HUI(2)*HA04+HUI(5)*HA05+HUI(8)*HA06
        H( 6)=HUI(3)*HA04+HUI(6)*HA05+HUI(9)*HA06
        H( 7)=HUI(1)*HA07+HUI(4)*HA08+HUI(7)*HA09
        H( 8)=HUI(2)*HA07+HUI(5)*HA08+HUI(8)*HA09
        H( 9)=HUI(3)*HA07+HUI(6)*HA08+HUI(9)*HA09
        H(10)=HUI(1)*HA10+HUI(4)*HA11+HUI(7)*HA12
        H(11)=HUI(2)*HA10+HUI(5)*HA11+HUI(8)*HA12
        H(12)=HUI(3)*HA10+HUI(6)*HA11+HUI(9)*HA12
        H(13)=HUI(1)*HA13+HUI(4)*HA14+HUI(7)*HA15
        H(14)=HUI(2)*HA13+HUI(5)*HA14+HUI(8)*HA15
        H(15)=HUI(3)*HA13+HUI(6)*HA14+HUI(9)*HA15
        H(16)=HUI(1)*HA16+HUI(4)*HA17+HUI(7)*HA18
        H(17)=HUI(2)*HA16+HUI(5)*HA17+HUI(8)*HA18
        H(18)=HUI(3)*HA16+HUI(6)*HA17+HUI(9)*HA18
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP03A(YA,HA)
      REAL YA(11),HA(18)
C
C AUXILIARY SUBROUTINE TO AP03, EVALUATES THE BASIS OF THE INTRINSIC
C RAY-CENTRED COORDINATE SYSTEM AT THE GIVEN POINT.
C
C INPUT:
C     YA...   QUANTITIES DESCRIBING A POINT OF A RAY
C
C OUTPUT:
C     HA...   COVARIANT (HA(1:9)) AND CONTRAVARIANT (HA(10:18))
C             COMPONENTS OF THE BASIS VECTORS OF THE INTRINSIC
C             RAY-CENTRED COORDINATE SYSTEM AT THE GIVEN POINT.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
*     INTEGER KOOR
*     EXTERNAL KOOR,METRIC,SMVPRD
C     KOOR,METRIC... FILE 'METRIC.FOR' OF THE PACKAGE 'MODEL'.
C     SMVPRD...FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C
C DATE: 1994, JANUARY 15
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR THE METRIC TENSOR ETC.:
C
*     REAL G(12),GAMMA(18),GSQRD
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /APCAUX/ BUT TO SHARE THE MEMORY.
*     COMMON/APCAUX/ G,GAMMA,GSQRD
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C
C.......................................................................
C
      REAL AUX
C
C.......................................................................
C
      HA(01)=YA( 9)
      HA(02)=YA(10)
      HA(03)=YA(11)
*     IF(KOOR().NE.0) THEN
C       CURVILINEAR COORDINATES:
*       CALL METRIC(YA(3),GSQRD,G,GAMMA)
*       CALL SMVPRD(G(7),HA(01),HA(02),HA(03),HA(10),HA(11),HA(12))
*       CALL SMVPRD(G(7),YA(6),YA(7),YA(8),HA(16),HA(17),HA(18))
*       AUX=SQRT(YA(6)*HA(16)+YA(7)*HA(17)+YA(8)*HA(18))
*       HA(07)=YA(6)/AUX
*       HA(08)=YA(7)/AUX
*       HA(09)=YA(8)/AUX
*       HA(16)=HA(16)/AUX
*       HA(17)=HA(17)/AUX
*       HA(18)=HA(18)/AUX
*       HA(04)=(HA(17)*HA(12)-HA(18)*HA(11))*GSQRD
*       HA(05)=(HA(18)*HA(10)-HA(16)*HA(12))*GSQRD
*       HA(06)=(HA(16)*HA(11)-HA(17)*HA(10))*GSQRD
*       HA(13)=(HA(08)*HA(03)-HA(09)*HA(02))/GSQRD
*       HA(14)=(HA(09)*HA(01)-HA(07)*HA(03))/GSQRD
*       HA(15)=(HA(07)*HA(02)-HA(06)*HA(01))/GSQRD
*     ELSE
C       CARTESIAN COORDINATES:
        AUX=SQRT(YA(6)*YA(6)+YA(7)*YA(7)+YA(8)*YA(8))
        HA(07)=YA(6)/AUX
        HA(08)=YA(7)/AUX
        HA(09)=YA(8)/AUX
        HA(10)=HA(01)
        HA(11)=HA(02)
        HA(12)=HA(03)
        HA(16)=HA(07)
        HA(17)=HA(08)
        HA(18)=HA(09)
        HA(04)=HA(17)*HA(12)-HA(18)*HA(11)
        HA(05)=HA(18)*HA(10)-HA(16)*HA(12)
        HA(06)=HA(16)*HA(11)-HA(17)*HA(10)
        HA(13)=HA(04)
        HA(14)=HA(05)
        HA(15)=HA(06)
*     END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP05(IUSER,HUI,Q11,Q21,Q12,Q22)
      INTEGER IUSER
      REAL HUI(5),Q11,Q21,Q12,Q22
C
C THIS SUBROUTINE EVALUATES THE COMPONENTS OF THE MATRIX OF GEOMETRICAL
C SPREADING AT A POINT SITUATED ON A RAY, SEE C.R.T.7.5.
C
C INPUT:
C     IUSER...IUSER=0... INTRINSIC CHOICE OF POLARIZATION VECTORS.
C               ANY OTHER INPUT NEED NOT BE SPECIFIED.
C             OTHERWISE, USER'S CHOICE OF POLARIZATION VECTORS AT THE
C               INITIAL POINT O/O OF THE RAY.
C               HUI(1:9) HAS TO BE SPECIFIED.
C     HUI(1),HUI(2),HUI(4),HUI(5)...  COMPONENTS HUI11,HUI21,HUI12,HUI22
C             OF THE 2*2 TRANSFORMATION MATRIX FROM THE INTRINSIC TO THE
C             USER'S RAY-CENTRED COORDINATE SYSTEM.
C
C OUTPUT:
C     Q11,Q21,Q12,Q22... COMPONENTS OF THE MATRIX OF GEOMETRICAL
C             SPREADING IN THE USER'S RAY-CENTRED COORDINATE SYSTEM AT
C             THE POINT O/S READ INTO THE COMMON BLOCK /POINTC/ BY THE
C             LAST INVOCATION OF THE SUBROUTINE AP00.
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
      REAL AUX
C
      Q11=Y(12)*YI(12)+Y(16)*YI(13)+Y(20)*YI(14)+Y(24)*YI(15)
      Q21=Y(13)*YI(12)+Y(17)*YI(13)+Y(21)*YI(14)+Y(25)*YI(15)
      Q12=Y(12)*YI(16)+Y(16)*YI(17)+Y(20)*YI(18)+Y(24)*YI(19)
      Q22=Y(13)*YI(16)+Y(17)*YI(17)+Y(21)*YI(18)+Y(25)*YI(19)
      IF(IUSER.NE.0) THEN
        AUX=Q11
        Q11=HUI(1)*AUX+HUI(4)*Q21
        Q21=HUI(2)*AUX+HUI(5)*Q21
        AUX=Q12
        Q12=HUI(1)*AUX+HUI(4)*Q22
        Q22=HUI(2)*AUX+HUI(5)*Q22
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP06(IUSER,HUI,P11,P21,P12,P22)
      INTEGER IUSER
      REAL HUI(5),P11,P21,P12,P22
C
C THIS SUBROUTINE EVALUATES THE COMPONENTS OF THE TRANSFORMATION MATRIX
C P AT A POINT SITUATED ON A RAY, SEE C.R.T.7.6.
C
C INPUT:
C     IUSER...IUSER=0... INTRINSIC CHOICE OF POLARIZATION VECTORS.
C               ANY OTHER INPUT NEED NOT BE SPECIFIED.
C             OTHERWISE, USER'S CHOICE OF POLARIZATION VECTORS AT THE
C               INITIAL POINT O/O OF THE RAY.
C               HUI(1:9) HAS TO BE SPECIFIED.
C     HUI(1),HUI(2),HUI(4),HUI(5)...  COMPONENTS HUI11,HUI21,HUI12,HUI22
C             OF THE 2*2 TRANSFORMATION MATRIX FROM THE INTRINSIC TO THE
C             USER'S RAY-CENTRED COORDINATE SYSTEM.
C
C OUTPUT:
C     P11,P21,P12,P22... COMPONENTS OF THE TRAMSFORMATION MATRIX P FROM
C             RAY COORDINATES TO THE USER'S RAY-CENTRED COMPONENTS OF
C             THE SLOWNESS VECTOR AT THE POINT O/S READ INTO THE COMMON
C             BLOCK /POINTC/ BY THE LAST INVOCATION OF THE SUBROUTINE
C             AP00.
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
      REAL AUX
C
      P11=Y(14)*YI(12)+Y(18)*YI(13)+Y(22)*YI(14)+Y(26)*YI(15)
      P21=Y(15)*YI(12)+Y(19)*YI(13)+Y(23)*YI(14)+Y(27)*YI(15)
      P12=Y(14)*YI(16)+Y(18)*YI(17)+Y(22)*YI(18)+Y(26)*YI(19)
      P22=Y(15)*YI(16)+Y(19)*YI(17)+Y(23)*YI(18)+Y(27)*YI(19)
      IF(IUSER.NE.0) THEN
        AUX=P11
        P11=HUI(1)*AUX+HUI(4)*P21
        P21=HUI(2)*AUX+HUI(5)*P21
        AUX=P12
        P12=HUI(1)*AUX+HUI(4)*P22
        P22=HUI(2)*AUX+HUI(5)*P22
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP07(QDET)
      REAL QDET
C
C THIS SUBROUTINE EVALUATES THE GEOMETRICAL SPREADING AT A POINT
C SITUATED ON A RAY, SEE C.R.T.7.7.
C
C NO INPUT.
C
C OUTPUT:
C     QDET... GEOMETRICAL SPREADING AT THE POINT O/S READ INTO THE
C             COMMON BLOCK /POINTC/ BY THE LAST INVOCATION OF THE
C             SUBROUTINE AP00.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL AP05
C     AP05... THIS FILE.
C
C DATE: 1991, MAY 15
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
      REAL DUMMY(5),Q11,Q21,Q12,Q22
C
      CALL AP05(0,DUMMY,Q11,Q21,Q12,Q22)
      QDET=ABS(Q11*Q22-Q12*Q21)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP08(IUSER,H,HUI,RM,RN)
      INTEGER IUSER
      REAL H(9),HUI(9),RM(6),RN(6)
C
C THIS SUBROUTINE EVALUATES THE COMPONENTS OF THE SYMMETRIC 3*3 MATRICES
C M AND N OF SECOND DERIVATIVES OF THE TRAVEL-TIME FIELD AT A POINT
C SITUATED ON A RAY, SEE C.R.T.7.8.  SUBROUTINE AP03 SHOULD BE CALLED
C BEFORE THE INVOCATION OF AP08 TO DEFINE INPUT ARGUMENTS H AND HUI, SEE
C BELOW.
C
C INPUT:
C     IUSER...IUSER=0... INTRINSIC CHOICE OF POLARIZATION VECTORS.
C               HUI(1:9) NEED NOT BE SPECIFIED.
C             OTHERWISE, USER'S CHOICE OF POLARIZATION VECTORS AT THE
C               INITIAL POINT O/O OF THE RAY.
C               HUI(1:9) HAS TO BE SPECIFIED.
C     H...    COVARIANT COMPONENTS OF THE BASIS VECTORS OF THE USER'S
C             RAY-CENTRED COORDINATE SYSTEM AT THE POINT O/S READ INTO
C             THE COMMON BLOCK /POINTC/ BY THE LAST INVOCATION OF THE
C             SUBROUTINE AP00.
C     HUI...  COMPONENTS OF THE 3*3 TRANSFORMATION MATRIX FROM THE
C             INTRINSIC TO THE USER'S RAY-CENTRED COORDINATE SYSTEM.
C
C OUTPUT:
C     RM...   COMPONENTS M11,M12,M22,M13,M23,M33 OF THE SECOND COVARIANT
C             DERIVATIVES OF THE TRAVEL-TIME FIELD IN THE USER'S
C             RAY-CENTRED COORDINATE SYSTEM, AT THE POINT O/S READ INTO
C             THE COMMON BLOCK /POINTC/ BY THE LAST INVOCATION OF THE
C             SUBROUTINE AP00.
C     RN...   COMPONENTS N11,N12,N22,N13,N23,N33 OF THE SECOND PARTIAL
C             DERIVATIVES OF THE TRAVEL-TIME FIELD IN THE GENERAL MODEL
C             COORDINATES, AT THE POINT O/S READ INTO THE COMMON BLOCK
C             /POINTC/ BY THE LAST INVOCATION OF THE SUBROUTINE AP00.
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
*     INTEGER KOOR
*     EXTERNAL KOOR,METRIC
C     KOOR,METRIC... FILE 'METRIC.FOR' OF THE PACKAGE 'MODEL'.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR THE METRIC TENSOR ETC.:
C
*     REAL G(12),GAMMA(18),GSQRD
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /APCAUX/ BUT TO SHARE THE MEMORY.
*     COMMON/APCAUX/ G,GAMMA,GSQRD
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C
C.......................................................................
C
      REAL QDET,Q11,Q21,Q12,Q22,P11,P21,P12,P22
      REAL RM1,RM2,RM3,AUX1,AUX2,AUX3
C
C     MATRIX M IN THE INTRINSIC R.C.C.S.
      CALL AP05(0,HUI,Q11,Q21,Q12,Q22)
      CALL AP06(0,HUI,P11,P21,P12,P22)
      QDET=Q11*Q22-Q12*Q21
      RM(1)=(P11*Q22-P12*Q21)/QDET
      RM(2)=(P21*Q22-P22*Q21)/QDET
      RM(3)=(P21*Q12-P22*Q11)/QDET
      IF(ICB1.GE.0) THEN
        AUX2=YL(1)**2
      ELSE
        AUX2=YL(2)**2
      END IF
C     SLOWNESS GRADIENT IN R.C.C.S.
      RM(4)=-(H(1)*YL(4)+H(2)*YL(5)+H(3)*YL(6))/AUX2
      RM(5)=-(H(4)*YL(4)+H(5)*YL(5)+H(6)*YL(6))/AUX2
      RM(6)=-(H(7)*YL(4)+H(8)*YL(5)+H(9)*YL(6))/AUX2
C
      IF(IUSER.NE.0) THEN
        AUX1=RM(1)*HUI(1)+RM(2)*HUI(4)+RM(4)*HUI(7)
        AUX2=RM(2)*HUI(1)+RM(3)*HUI(4)+RM(5)*HUI(7)
        AUX3=RM(4)*HUI(1)+RM(5)*HUI(4)+RM(6)*HUI(7)
        RM1  =HUI(1)*AUX1+HUI(4)*AUX2+HUI(7)*AUX3
        AUX1=RM(1)*HUI(2)+RM(2)*HUI(5)+RM(4)*HUI(8)
        AUX2=RM(2)*HUI(2)+RM(3)*HUI(5)+RM(5)*HUI(8)
        AUX3=RM(4)*HUI(2)+RM(5)*HUI(5)+RM(6)*HUI(8)
        RM2  =HUI(1)*AUX1+HUI(4)*AUX2+HUI(7)*AUX3
        RM3  =HUI(2)*AUX1+HUI(5)*AUX2+HUI(8)*AUX3
        AUX1=RM(1)*HUI(3)+RM(2)*HUI(6)+RM(4)*HUI(9)
        AUX2=RM(2)*HUI(3)+RM(3)*HUI(6)+RM(5)*HUI(9)
        AUX3=RM(4)*HUI(3)+RM(5)*HUI(6)+RM(6)*HUI(9)
        RM(4)=HUI(1)*AUX1+HUI(4)*AUX2+HUI(7)*AUX3
        RM(5)=HUI(2)*AUX1+HUI(5)*AUX2+HUI(8)*AUX3
        RM(6)=HUI(3)*AUX1+HUI(6)*AUX2+HUI(9)*AUX3
        RM(1)=RM1
        RM(2)=RM2
        RM(3)=RM3
      END IF
      AUX1=RM(1)*H(1)+RM(2)*H(4)+RM(4)*H(7)
      AUX2=RM(2)*H(1)+RM(3)*H(4)+RM(5)*H(7)
      AUX3=RM(4)*H(1)+RM(5)*H(4)+RM(6)*H(7)
      RN(1)=H(1)*AUX1+H(4)*AUX2+H(7)*AUX3
      AUX1=RM(1)*H(2)+RM(2)*H(5)+RM(4)*H(8)
      AUX2=RM(2)*H(2)+RM(3)*H(5)+RM(5)*H(8)
      AUX3=RM(4)*H(2)+RM(5)*H(5)+RM(6)*H(8)
      RN(2)=H(1)*AUX1+H(4)*AUX2+H(7)*AUX3
      RN(3)=H(2)*AUX1+H(5)*AUX2+H(8)*AUX3
      AUX1=RM(1)*H(3)+RM(2)*H(6)+RM(4)*H(9)
      AUX2=RM(2)*H(3)+RM(3)*H(6)+RM(5)*H(9)
      AUX3=RM(4)*H(3)+RM(5)*H(6)+RM(6)*H(9)
      RN(4)=H(1)*AUX1+H(4)*AUX2+H(7)*AUX3
      RN(5)=H(2)*AUX1+H(5)*AUX2+H(8)*AUX3
      RN(6)=H(3)*AUX1+H(6)*AUX2+H(9)*AUX3
*     IF(KOOR().NE.0) THEN
C       CURVILINEAR COORDINATES:
*       CALL METRIC(Y(3),GSQRD,G,GAMMA)
*       RN(1)=RN(1)+GAMMA(1)*Y(6)+GAMMA( 7)*Y(6)+GAMMA(13)*Y(8)
*       RN(2)=RN(2)+GAMMA(2)*Y(6)+GAMMA( 8)*Y(6)+GAMMA(14)*Y(8)
*       RN(3)=RN(3)+GAMMA(3)*Y(6)+GAMMA( 9)*Y(6)+GAMMA(15)*Y(8)
*       RN(4)=RN(4)+GAMMA(4)*Y(6)+GAMMA(10)*Y(6)+GAMMA(16)*Y(8)
*       RN(5)=RN(5)+GAMMA(5)*Y(6)+GAMMA(11)*Y(6)+GAMMA(17)*Y(8)
*       RN(6)=RN(6)+GAMMA(6)*Y(6)+GAMMA(12)*Y(6)+GAMMA(18)*Y(8)
*     END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE AP15(IUSER,HI,H,HUI,ZI,Z,AMPLI,AMPL)
      INTEGER IUSER
      REAL HI(18),H(18),HUI(9),ZI(9),Z(9),AMPLI(6),AMPL(6)
C
C THIS SUBROUTINE EVALUATES THE RAY AMPLITUDES AT A POINT SITUATED ON A
C RAY, SEE C.R.T.7.15.  SUBROUTINE AP03 SHOULD BE CALLED BEFORE THE
C INVOCATION OF AP15 TO DEFINE THE INPUT ARGUMENTS HI AND H, SEE BELOW.
C
C INPUT:
C     IUSER...IUSER=0... INTRINSIC CHOICE OF POLARIZATION VECTORS.
C               HUI(1:9) NEED NOT BE SPECIFIED.
C             OTHERWISE, USER'S CHOICE OF POLARIZATION VECTORS AT THE
C               INITIAL POINT O/O OF THE RAY.
C               HUI(1:9) HAS TO BE SPECIFIED.
C     HI...   COVARIANT (HI(1:9)) AND CONTRAVARIANT (HI(10:18))
C             COMPONENTS OF THE BASIS VECTORS OF THE USER'S RAY-CENTRED
C             COORDINATE SYSTEM AT THE INITIAL POINT O/O OF A RAY.
C     H...    COVARIANT (H(1:9)) AND CONTRAVARIANT (H(10:18)) COMPONENTS
C             OF THE BASIS VECTORS OF THE USER'S RAY-CENTRED COORDINATE
C             SYSTEM AT THE POINT O/S READ INTO THE COMMON BLOCK
C             /POINTC/ BY THE LAST INVOCATION OF THE SUBROUTINE AP00.
C     HUI...  COMPONENTS OF THE 3*3 TRANSFORMATION MATRIX FROM THE
C             INTRINSIC TO THE USER'S RAY-CENTRED COORDINATE SYSTEM.
C     ZI...   CONTRAVARIANT COMPONENTS OF THE BASIS VECTORS OF THE LOCAL
C             COORDINATE SYSTEM AT THE INITIAL POINT O/O OF A RAY.
C             SEE C.R.T.7.15, EQ.(7.44).
C     Z...    CONTRAVARIANT COMPONENTS OF THE BASIS VECTORS OF THE LOCAL
C             RECORDING COORDINATE SYSTEM AT THE POINT O/S READ INTO THE
C             COMMON BLOCK /POINTC/ BY THE LAST INVOCATION OF THE
C             SUBROUTINE AP00.  SEE C.R.T.7.15, EQ.(7.43).
C     AMPLI...COMPONENTS RE(A1),IM(A1),RE(A2),IM(A2),RE(A3),IM(A3) OF
C             THE RAY AMPLITUDE IN THE LOCAL COORDINATE SYSTEM
C             MULTIPLIED BY THE SQUARE ROOT OF THE GEOMETRICAL
C             SPREADING, AT THE INITIAL POINT O/O OF A RAY.
C             SEE C.R.T.7.15, EQ.(7.47).
C
C OUTPUT:
C     AMPL... COMPONENTS RE(A1),IM(A1),RE(A2),IM(A2),RE(A3),IM(A3) OF
C             THE RAY AMPLITUDE IN THE LOCAL RECORDING COORDINATE SYSTEM
C             AT THE POINT O/S READ INTO THE COMMON BLOCK /POINTC/ BY
C             THE LAST INVOCATION OF THE SUBROUTINE AP00.
C             SEE C.R.T.7.15, EQ.(7.45).
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL AP07
C     AP05... THIS FILE.
C     AP07... THIS FILE.
C
C ERROR MESSAGES:
C     715...  USER'S CHOICE OF R.C.C.S. NOT ALLOWED:
C             NONZERO INPUT PARAMETER IUSER OF THE SUBROUTINE AP15
C             INDICATES USER'S CHOICE OF POLARIZATION VECTORS.  THIS
C             OPTION HAS NOT BEEN INCLUDED IN THIS SUBROUTINE YET.
C             SUBROUTINE AP03 SHOULD BE CALLED BEFORE THE INVOCATION OF
C             AP15 WITH IUSER=0.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
      INTEGER I
      REAL AUX01,AUX02,AUX03,AUX04,AUX05,AUX06,AUX07,AUX08,AUX09,AUX10
      REAL AUX11,AUX12,AUX13,AUX14,AUX15,AUX16,AUX17,AUX18,AUX19,SCALAR
C
      IF(IUSER.NE.0) THEN
        PAUSE'ERROR 715 IN AP15: USER''S CHOICE OF R.C.C.S. NOT ALLOWED'
        STOP
      END IF
C
C     SCALAR MULTIPLICATION FACTOR IN EQ.(7.45)
      CALL AP07(SCALAR)
      SCALAR=YLI(3)/(YL(3)*SCALAR)
      IF(ICB1.GT.0) THEN
C       P-WAVE
        SCALAR=SCALAR/YL(1)
      ELSE
C       S-WAVE
        SCALAR=SCALAR/YL(2)
      END IF
C     NOTE: THE SQUARE ROOT AND VELOCITY AT THE INITIAL POINT OF THE RAY
C     WILL BE APPLIED LATER ON.
C
C     TRANSPOSED INVERSE MATRIX TO Z MULTIPLIED BY ITS DETERMINANT
      AUX01=Z(5)*Z(9)-Z(6)*Z(8)
      AUX02=Z(6)*Z(7)-Z(4)*Z(9)
      AUX03=Z(4)*Z(8)-Z(5)*Z(7)
      AUX04=Z(3)*Z(8)-Z(2)*Z(9)
      AUX05=Z(1)*Z(9)-Z(3)*Z(7)
      AUX06=Z(2)*Z(7)-Z(1)*Z(8)
      AUX07=Z(2)*Z(6)-Z(3)*Z(5)
      AUX08=Z(3)*Z(4)-Z(1)*Z(6)
      AUX09=Z(1)*Z(5)-Z(2)*Z(4)
      AUX10=Z(1)*AUX01+Z(2)*AUX02+Z(3)*AUX03
C
C     TRANSFORMATION MATRIX FROM R.C.C.S. TO THE LOCAL RECORDING
C     COORDINATE SYSTEM Z
      IF(ICB1.GT.0.OR.NY.EQ.33.OR.NY.EQ.39) THEN
C       P-WAVE OR INTERFACE CONVERSION COEFFICIENTS APPLIED
        AUX17=AUX01*H(16)+AUX04*H(17)+AUX07*H(18)
        AUX18=AUX02*H(16)+AUX05*H(17)+AUX08*H(18)
        AUX19=AUX03*H(16)+AUX06*H(17)+AUX09*H(18)
      END IF
      IF(ICB1.LT.0.OR.NY.EQ.33.OR.NY.EQ.39) THEN
C       S-WAVE OR INTERFACE CONVERSION COEFFICIENTS APPLIED
        AUX11=AUX01*H(10)+AUX04*H(11)+AUX07*H(12)
        AUX12=AUX02*H(10)+AUX05*H(11)+AUX08*H(12)
        AUX13=AUX03*H(10)+AUX06*H(11)+AUX09*H(12)
        AUX14=AUX01*H(13)+AUX04*H(14)+AUX07*H(15)
        AUX15=AUX02*H(13)+AUX05*H(14)+AUX08*H(15)
        AUX16=AUX03*H(13)+AUX06*H(14)+AUX09*H(15)
      END IF
C
      IF(ICB1I.GT.0) THEN
C       P-WAVE AT THE INITIAL POINT:
        SCALAR=SQRT(SCALAR*YLI(1))/AUX10
C       MATRIX ZI IN R.C.C.S.
        AUX01=HI(16)*ZI(1)+HI(17)*ZI(2)+HI(18)*ZI(3)
        AUX02=HI(16)*ZI(4)+HI(17)*ZI(5)+HI(18)*ZI(6)
        AUX03=HI(16)*ZI(7)+HI(17)*ZI(8)+HI(18)*ZI(9)
C       AMPLI IN R.C.C.S. MULTIPLIED BY SCALAR
        AUX07=SCALAR*(AUX01*AMPLI(1)+AUX02*AMPLI(3)+AUX03*AMPLI(5))
        AUX08=SCALAR*(AUX01*AMPLI(2)+AUX02*AMPLI(4)+AUX03*AMPLI(6))
C       AMPL IN R.C.C.S.
        AUX01=Y(28)*AUX07-Y(29)*AUX08
        AUX02=Y(29)*AUX07+Y(28)*AUX08
        IF(NY.GT.29) THEN
C         P-WAVE TO S-WAVE OR INTERFACE CONVERSION COEFFICIENTS APPLIED
          AUX03=Y(30)*AUX07-Y(31)*AUX08
          AUX04=Y(31)*AUX07+Y(30)*AUX08
          IF(NY.GT.31) THEN
C           INTERFACE CONVERSION COEFFICIENTS APPLIED
            AUX05=Y(32)*AUX07-Y(33)*AUX08
            AUX06=Y(33)*AUX07+Y(32)*AUX08
          END IF
        END IF
      ELSE
C       S-WAVE AT THE INITIAL POINT:
        SCALAR=SQRT(SCALAR*YLI(2))/AUX10
C       MATRIX ZI IN R.C.C.S.
        AUX01=HI(10)*ZI(1)+HI(11)*ZI(2)+HI(12)*ZI(3)
        AUX02=HI(10)*ZI(4)+HI(11)*ZI(5)+HI(12)*ZI(6)
        AUX03=HI(10)*ZI(7)+HI(11)*ZI(8)+HI(12)*ZI(9)
        AUX04=HI(13)*ZI(1)+HI(14)*ZI(2)+HI(15)*ZI(3)
        AUX05=HI(13)*ZI(4)+HI(14)*ZI(5)+HI(15)*ZI(6)
        AUX06=HI(13)*ZI(7)+HI(14)*ZI(8)+HI(15)*ZI(9)
C       AMPLI IN R.C.C.S. MULTIPLIED BY SCALAR
        AUX07=SCALAR*(AUX01*AMPLI(1)+AUX02*AMPLI(3)+AUX03*AMPLI(5))
        AUX08=SCALAR*(AUX01*AMPLI(2)+AUX02*AMPLI(4)+AUX03*AMPLI(6))
        AUX09=SCALAR*(AUX04*AMPLI(1)+AUX05*AMPLI(3)+AUX06*AMPLI(5))
        AUX10=SCALAR*(AUX04*AMPLI(2)+AUX05*AMPLI(4)+AUX06*AMPLI(6))
C       AMPL IN R.C.C.S.
        I=(NY+27)/2
        AUX01=Y(28)*AUX07-Y(29)*AUX08+Y(I+1)*AUX09-Y(I+2)*AUX10
        AUX02=Y(29)*AUX07+Y(28)*AUX08+Y(I+2)*AUX09+Y(I+1)*AUX10
        IF(NY.GT.31) THEN
C         S-WAVE TO S-WAVE OR INTERFACE CONVERSION COEFFICIENTS APPLIED
          AUX03=Y(30)*AUX07-Y(31)*AUX08+Y(I+3)*AUX09-Y(I+4)*AUX10
          AUX04=Y(31)*AUX07+Y(30)*AUX08+Y(I+4)*AUX09-Y(I+3)*AUX10
          IF(NY.GT.35) THEN
C           INTERFACE CONVERSION COEFFICIENTS APPLIED
            AUX05=Y(32)*AUX07-Y(33)*AUX08+Y(38)*AUX09-Y(39)*AUX10
            AUX06=Y(33)*AUX07+Y(32)*AUX08+Y(39)*AUX09-Y(38)*AUX10
          END IF
        END IF
      END IF
C
C     AMPL IN THE LOCAL RECORDING COORDINATE SYSTEM Z
      IF(NY.EQ.33.OR.NY.EQ.39) THEN
C       INTERFACE CONVERSION COEFFICIENTS APPLIED
        AMPL(1)=AUX11*AUX01+AUX14*AUX03+AUX17*AUX05
        AMPL(2)=AUX11*AUX02+AUX14*AUX04+AUX17*AUX06
        AMPL(3)=AUX12*AUX01+AUX15*AUX03+AUX18*AUX05
        AMPL(4)=AUX12*AUX02+AUX15*AUX04+AUX18*AUX06
        AMPL(5)=AUX13*AUX01+AUX16*AUX03+AUX19*AUX05
        AMPL(6)=AUX13*AUX02+AUX16*AUX04+AUX19*AUX06
      ELSE IF(ICB1.GT.0) THEN
C       P-WAVE
        AMPL(1)=AUX17*AUX01
        AMPL(2)=AUX17*AUX02
        AMPL(3)=AUX18*AUX01
        AMPL(4)=AUX18*AUX02
        AMPL(5)=AUX19*AUX01
        AMPL(6)=AUX19*AUX02
      ELSE
C       S-WAVE
        AMPL(1)=AUX11*AUX01+AUX14*AUX03
        AMPL(2)=AUX11*AUX02+AUX14*AUX04
        AMPL(3)=AUX12*AUX01+AUX15*AUX03
        AMPL(4)=AUX12*AUX02+AUX15*AUX04
        AMPL(5)=AUX13*AUX01+AUX16*AUX03
        AMPL(6)=AUX13*AUX02+AUX16*AUX04
      END IF
      RETURN
      END
C
C=======================================================================
C