C                                                            <HTML><PRE>
C     PROGRAM  A N R A Y, VERSION 4.01  (PRAHA, SEPTEMBER 1997)
C
C*******************************************************************
C
C     PROGRAM ANRAY IS DESIGNED FOR RAY, TRAVEL TIME AND
C     AMPLITUDE COMPUTATIONS IN 3D GENERAL ANISOTROPIC AND ISOTROPIC
C     LATERALLY VARYING LAYERED MEDIA. THE PROGRAM MAKES POSSIBLE
C     COMPUTATION OF RAYS SPECIFIED BY INITIAL ANGLES AT THE SOURCE,
C     I.E., INITIAL-VALUE RAY TRACING, OR RAYS STARTING FROM THE
C     SOURCE AND TERMINATING ON A VERTICAL OR SURFACE PROFILE, I.E.
C     BOUNDARY-VALUE RAY TRACING. RAY AMPLITUDES CAN BE COMPUTED
C     ALONG RAYS.
C
C*******************************************************************
C
C
      CHARACTER*80 MTEXT,FILEIN,FILEOU,FILE1,FILE2,FILE3
      DIMENSION Y(18)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      COMMON /AUXX/  MMX(20),MMY(20),MMXY(20)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST,
     1xprf,yprf
      COMMON /DENS/ RHO(20)
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /ZERO/ RNULL
      COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP
C
C**************************************************
C
      LIN=5
      LOU=6
      LU1=1
      LU2=2
      LU3=3
      FILEIN='anray.dat'
      FILEOU='anray.out'
      FILE1='lu1.out'
      FILE2='lu2.out'
      FILE3='lu3.out'
      WRITE(*,'(2A)') ' Specify names of input and output files',
     1' LIN, LOU, LU1, LU2, LU3: '
      READ(*,*) FILEIN,FILEOU,FILE1,FILE2,FILE3
      IF(FILE1.EQ.' ') LU1=0
      IF(FILE2.EQ.' ') LU2=0
      IF(FILE3.EQ.' ') LU3=0
      LOUT=LOU
      OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LOU,FILE=FILEOU,FORM='FORMATTED')
      IF(LU1.NE.0)OPEN(LU1,FILE=FILE1,FORM='FORMATTED')
      IF(LU2.NE.0)OPEN(LU2,FILE=FILE2,FORM='FORMATTED')
      IF(LU3.NE.0)OPEN(LU3,FILE=FILE3,FORM='FORMATTED')
C
C**************************************************
C
      WRITE(LOU,777)
 777  FORMAT(///,'***********************'
     1,//,'  PROGRAM   A N R A Y      ',//,
     2'***********************',//)
      NCODE=1
    1 READ(LIN,115)MTEXT
      WRITE(LOU,115)MTEXT
      READ(LIN,106)INULL,ISURF
      IF(INULL.EQ.0)INULL=4
      RNULL=10.**(-INULL)
      WRITE(LOU,106)INULL,ISURF
C
C
C     SPECIFICATION OF THE MODEL
C
      CALL MODEL(MTEXT,LIN)
C
C     GENERATE FILE FOR PLOTTING VARIOUS CHARACTERISTIC SURFACES
C
      IF(ISURF.EQ.1)THEN
        CALL SURFPL(LIN,LU3)
      END IF
C
C     SPECIFICATION OF SYNTHETIC SEISMOGRAMS
C
    2 READ(LIN,100)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX,
     1IPOL,IPREC,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI
      WRITE(LOU,102)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX,
     1IPOL,IPREC,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI
      IF(ICONT.EQ.(-1))GO TO 1
      IF(ICONT.EQ.0)GO TO 99
C
C
      if(mep.ne.0.and.mdim.eq.0)mdim=1
      IF(ITMAX.EQ.0)ITMAX=10
      IVSP=0
      IF(ILOC.EQ.0)ITPR=3
      IF(ILOC.EQ.1)THEN
        IVSP=1
        ITPR=43
        MREG=1
      END IF
      IF(ILOC.GT.1)THEN
        ITPR=ILOC+100
        MREG=1
      END IF
C
      IF(MEP.EQ.0)THEN
        NDST=0
      END IF
C
      if(mep.eq.1)then
        ndst=1
        read(lin,104)xrec,yrec
        write(lou,104)xrec,yrec
        go to 4
      end if
      IF(MEP.LT.0)THEN
        NDST=-MEP
        READ(LIN,104)PROF(1),(DST(I),I=1,NDST),xprf,yprf
        WRITE(LOU,104)PROF(1),(DST(I),I=1,NDST),xprf,yprf
        IF(NDST.EQ.1)RSTEP=1.
        IF(NDST.EQ.1)DST(2)=DST(1)+1.
        IF(NDST.EQ.1)GO TO 4
        RSTEP=(DST(NDST)-DST(1))/FLOAT(NDST-1)
      END IF
C
      IF(MEP.GT.0)THEN
        NDST=MEP
        READ(LIN,104)PROF(1),RMIN,RSTEP,xprf,yprf
        WRITE(LOU,104)PROF(1),RMIN,RSTEP,xprf,yprf
        DO 13 I=1,MEP
   13   DST(I)=RMIN+(I-1)*RSTEP
        if(ndst.eq.1)dst(2)=rmin+rstep
      END IF
      PROF(2)=PROF(1)+1.
      NDSTP=1
C
      IF(IVSP.EQ.1.AND.NDST.NE.0)THEN
        READ(LIN,104)XVSP,YVSP
        WRITE(LOU,104)XVSP,YVSP
      END IF
C
    4 READ(LIN,104)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS
      WRITE(LOU,104)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS
C
      if(abs(xprf).lt..000001.and.abs(yprf).lt..000001)then
        xprf=xsour
        yprf=ysour
      end if
      if(mep.eq.1)then
        xe=xrec-xprf
        ye=yrec-yprf
        rprf=sqrt(xe*xe+ye*ye)
        xatan=atan2(ye,xe)
        prof(1)=xatan
        rmin=rprf
        dst(1)=rmin
        write(lou,104)rprf,xatan
        RSTEP=100.
        DST(2)=DST(1)+100.
        prof(2)=prof(1)+1.
        ndstp=1
      end if
C
      IF(IVSP.EQ.1.AND.NDST.NE.0)THEN
        XNRM=XVSP-XSOUR
        YNRM=YVSP-YSOUR
        aux=sqrt(xnrm*xnrm+ynrm*ynrm)
        xnrm=xnrm/aux
        ynrm=ynrm/aux
        PROF(1)=ATAN2(YNRM,XNRM)
        PROF(2)=PROF(1)+1.
        xprf=xsour
        yprf=ysour
      END IF
      IF(MCOD.EQ.0)THEN
        READ(LIN,104)AMIN,ASTEP,AMAX
        WRITE(LOU,104)AMIN,ASTEP,AMAX
        READ(LIN,104)BMIN,BSTEP,BMAX
        if(abs(bstep).lt..000001)then
          bmin=prof(1)-.3
          bmax=prof(1)+.4
          bstep=.6
        end if
        WRITE(LOU,104)BMIN,BSTEP,BMAX
      END IF
      IF(MREG.EQ.0.AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)')
     1 ' COEFFICIENTS OF CONVERSION ARE APPLIED'
      IF(MREG.NE.0.AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)')
     1 ' COEFFICIENTS OF CONVERSION ARE *** NOT *** APPLIED'
      IF(DT.LT.0.0000001)DT=1.
      IF(AC.LT.0.0000001)AC=0.0001
      TMAX=10000.
      IF(REPS.LT..0000001)REPS=.05
      IF(PREPS.LT..0000001)PREPS=.05
      IND=-1
      NDER=1
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,AMIN1,BMIN,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      Y(1)=XSOUR
      Y(2)=YSOUR
      Y(3)=ZSOUR
      IF(IND.EQ.50)WRITE(LOU,111)IND
      IF(IND.EQ.50)GO TO 99
      LAY=IND
      ISOUR=IND
      itype=3
      CALL PARDIS(Y,0)
      VP=SQRT(A11)
      IF(IRHO.EQ.0)RO=1.7+.2*VP
      IF(IRHO.EQ.1)RO=RHO(IND)
C
C     GENERATE FILE LU2 FOR SYNTHETIC SEISMOGRAM COMPUTATIONS
C
      IF(LU2.NE.0.AND.NDST.NE.0)THEN
        WRITE(LU2,115)MTEXT
        KSH=2
        WRITE(LU2,100)NDST,KSH,ILOC
        WRITE(LU2,104)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,RO
        if(mep.ne.1)WRITE(LU2,104)(DST(I),I=1,NDST)
        if(mep.eq.1)WRITE(LU2,104)xrec,yrec,rprf,xatan
      END IF
C
C    LOOP FOR ELEMENTARY WAVES
C
   20 READ(LIN,100)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF)
      WRITE(LOU,100)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF)
      IF(KREF.EQ.0)GOTO 2
      IF(MOUT.NE.0)WRITE(LOU,107)
      WRITE(LOU,103)NCODE,KC,KREF,((CODE(I,K),K=1,2),I=1,KREF)
C
      IF(MCOD.NE.0)THEN
        READ(LIN,104)AMIN,ASTEP,AMAX
        WRITE(LOU,104)AMIN,ASTEP,AMAX
        READ(LIN,104)BMIN,BSTEP,BMAX
        if(abs(bstep).lt..000001)then
          bmin=prof(1)-.3
          bmax=prof(1)+.4
          bstep=.6
        end if
        WRITE(LOU,104)BMIN,BSTEP,BMAX
      END IF
C
C     GENERATE FILE LU1 FOR PLOTTING OF RAY DIAGRAMS,
C     TIME-DISTANCE AND AMPLITUDE-DISTANCE CURVES
C
      IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 21
      WRITE(LU1,100)ICONT,NDST,ILOC
      WRITE(LU1,104)RO
      NPN=2
      APN=0.
      WRITE(LU1,100)NPN,NPN,NPN
      WRITE(LU1,101)APN,APN,APN,APN,APN
      WRITE(LU1,101)APN,APN,APN,APN,APN
      WRITE(LU1,104)Xprf,Yprf,0.0,PROF(1)
      WRITE(LU1,104)(DST(I),I=1,NDST)
   21 CONTINUE
C
C
C     SEARCH FOR THE NUMBER OF THE ELEMENT OF THE RAY, STARTING FROM
C     WHICH THE WAVE DOES UNDERTAKE NEITHER REFLECTION NOR CONVERSION
C
      ICOD=0
      IF(IVSP.EQ.0)GO TO 35
      DO 34 I=1,KREF
      ICOD=KREF-I+1
      IF(ICOD.EQ.1) GO TO 34
      IC1=CODE(ICOD,1)
      IC2=CODE(ICOD-1,1)
      IF((IC1-IC2).EQ.0)GO TO 35
      IC1=CODE(ICOD,2)
      IC2=CODE(ICOD-1,2)
      IF((IC1-IC2).NE.0)GO TO 35
   34 CONTINUE
   35 CONTINUE
      IF(MOUT.NE.0)WRITE(LOU,108)
C
C
      CALL RECEIV(XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,ITMAX,AMIN,ASTEP,
     1AMAX,BMIN,BSTEP,BMAX,MOUT,LU1,LU2,METHOD,ITPR,NCODE)
      IF(IND.EQ.20)WRITE(LOU,111)IND
      IF(IND.EQ.20)GO TO 99
      IF(IND.EQ.14) WRITE(LOU,111) IND
      NCODE=NCODE+1
      GOTO 20
C
C     END OF LOOP FOR ELEMENTARY WAVES
C
C
  100 FORMAT(26I3)
  101 FORMAT(5E15.5)
  102 FORMAT(1H0,////,2X,26I3)
  104 FORMAT(8F10.5)
  106 FORMAT(17I5)
  103 FORMAT(4X,I4,9X,36I3)
  107 FORMAT(//2X,'INT.CODE',5X,'E X T E R N A L   C O D E')
  108 FORMAT(//)
  111 FORMAT(/2X,'IND=',I5,/)
  115 FORMAT(A)
C
   99 CONTINUE
      IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)ICONT,ICONT
      IF(LU1.NE.0)REWIND LU1
      IF(LU2.NE.0)REWIND LU2
C
      STOP
      END
C
C     *********************************************************
C
      SUBROUTINE AMPL (AMPX,AMPY,AMPZ,UU)
C
C     ROUTINE FOR COMPUTING COMPLEX VECTORIAL RAY AMPLITUDES
C
C     OUTPUT PARAMETERS
C     AMPX(2),AMPY(2),AMPZ(2) - X,Y AND Z COMPONENTS OF COMPLEX
C     VECTORIAL RAY AMPLITUDES IN THE MODEL COORDINATES. FOR P WAVE
C     IN ANY MEDIUM AND FOR S WAVES IN AN ANISOTROPIC MEDIUM, I=1.
C     FOR S WAVE GENERATED IN AN ISOTROPIC MEDIUM, I=1,2. I=1 AND 2
C     CORRESPOND TO S WAVES SPECIFIED AT THE SOURCE BY VECTORS E1
C     E2. VECTORS E1 AND E2 TOGETHER WITH UNIT VECTOR TANGENT TO
C     THE RAY FORM A BASIS OF RAY CENTRED COORDINATE SYSTEM.
C     UU - PRODUCT OF RATIOS OF DENSITIES AND COSINES OF INCIDENCE
C     AND OF REFLECTION/TRANSMISSION AT POINTS WHERE THE RAY CROSSES
c     INTERFACES.
C
C     CALLED FROM: RECEIV
C     ROUTINES CALLED: POLAR,TRANSL,COEF
C
      DIMENSION Y(18),UN(3),POLD(3),PNEW(3)
      COMPLEX  AMPX(2),AMPY(2),AMPZ(2),CR(3),UC(3),STU(6)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DENS/ RHO(20)
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /RAY2/  DRY(3,400)
C
      KSS=1
      ISHEAR=0
      ITYPE=CODE(1,2)
      IF(IANI(ISOUR).EQ.0.AND.ITYPE.NE.3)then
        ISHEAR=1
        itype=1
      end if
      ITP=ITYPE
      do 1 i=1,2
      AMPX(i)=CMPLX(0.,0.)
      AMPY(i)=CMPLX(0.,0.)
      AMPZ(i)=CMPLX(0.,0.)
    1 continue
C
 3000 NN=N
      IDD=0
      N2=0
      N1=1
      IRE=IREF
      AV=1.
C
C  SPECIFICATION OF DISPLACEMENT VECTOR AT SOURCE
C  IN RAY CENTERED COORDINATES
C
      DO 5 I=1,3
      CR(I)=(0.,0.)
    5 CONTINUE
      CR(ITP)=(1.,0.)
      IREF1=IREF-1
      IF(IREF1.EQ.0) GOTO 100
C
C  LOOP OVER INTERFACES
C
      DO 10 I=1,IREF1
      IREF=I
      IF(KC.NE.0) ITYPE=CODE(IREF,2)
      N=KINT(IREF)
      IF(N.EQ.0) THEN
        IDD=1
        GO TO 10
      ELSE
        N1=N2+1
        N2=N
        IF(IDD.NE.0) N2=-N2
        IDD=0
C
C  COMPUTATION OF POLARIZATION VECTORS
C  CONSIDERED POLARIZATION VECTOR(S) ARE STORED IN CORRESPONDING
C  COLUMNS OF THE MATRIX HHH. OTHER COLUMNS ARE ZERO.
C
        CALL POLAR(N1,N2,NN,IREF)
      END IF
      DO 20 K=1,6
      Y(K)=AY(K+1,N)
  20  CONTINUE
      IF(IAMP.GT.0)WRITE(LOU,'(a,2i5,6f10.5)')' AMPL:I,N,Y',I,N,
     1(Y(L),L=1,6)
      DO 30 K=1,3
      POLD(K)=Y(K+3)
      PS(K,7,IREF)=Y(K+3)
  30  CONTINUE
      DO 40 K=1,3
      UN(K)=DS(K,IREF)
  40  CONTINUE
      LAY=IS(1,IREF)
      ITRANS=IS(2,IREF)
      ITR1=ITRANS
      IF(UN(3).GT.0.0) GOTO 50
C
C  RAY STRIKING THE INTERFACE FROM ABOVE
C
      IF(ITRANS.EQ.0) THEN
        LAY=LAY+1
        ITRANS=1
        GOTO 70
      END IF
      IF(ITRANS.GT.0) THEN
        LAY=LAY-1
        ITRANS=0
        GOTO 70
      END IF
C
C  RAY STRIKING THE INTERFACE FROM BELOW
C
  50  IF(ITRANS.EQ.0) THEN
        LAY=LAY-1
        ITRANS=1
        GOTO 70
      END IF
      IF(ITRANS.GT.0) THEN
        LAY=LAY+1
        ITRANS=0
        GOTO 70
      END IF
C
C  SLOWNESS VECTORS ON THE SIDE OF THE INTERFACE WHERE GENERATED
C  WAVE PROPAGATES WERE DETERMINED DURING THE CALL OF TRANSL IN THE
C  ROUTINE OUT. HERE REMAINING SLOWNESS VECTORS ON THE OTHER SIDE
C  OF THE INTERFACE ARE DETERMINED
C
C  REDEFINITION OF IREF FOR CALL OF ROUTINE TRANSL
C
  70  IF(LAY.EQ.0) THEN
        DO 71 K=4,6
        DO 71 L=1,3
        PS(L,K,IREF)=CMPLX(0.,0.)
  71    CONTINUE
        GO TO 75
      END IF
      IREF=IREF+1
      CALL TRANSL(Y,POLD,PNEW,UN,ITRANS,0)
      IF(IND.EQ.10)RETURN
      IREF=IREF-1
  75  IF(IAMP.ne.0)then
        WRITE(LOU,'(a)')' REFLECTED/TRANSMITTED SLOWNESS VECTORS'
        WRITE(LOU,'(6f12.6)')((PS(L,K,IREF),L=1,3),k=1,6)
      end if
      AV1=(DS(11,IREF)*ds(10,iref))/(DS(8,IREF)*ds(7,iref))
      AV=AV*AV1
      IF(IAMP.GT.0) THEN
        WRITE(LOU,'(A)') 'ROI,ROG,UNVGI,UNVGG,AV1,AV'
        WRITE(LOU,'(6F10.5)') DS(8,IREF),
     1  DS(11,IREF),ds(7,iref),ds(10,iref),AV1,AV
        WRITE(LOU,'(A,/,6F12.5,/,3(3F12.5/))') ' CR,HHH',
     2  CR,((HHH(J,K),J=1,3),K=1,3)
      END IF
C
C  COMPUTATION OF AMPLITUDE COEFFICIENTS OF REFLECTED/TRANSMITTED WAVES
C
C
C  COMPUTATION OF CARTESIAN COMPONENTS OF INCIDENT DISPLACEMENT VECTOR
C
      DO 87 K=1,3
      STU(K)=CMPLX(0.,0.)
      DO 87 J=1,3
      STU(K)=HHH(J,K)*CR(J)+STU(K)
  87  CONTINUE
      IF(IAMP.GT.0)WRITE(LOU,'(A,6F10.5)') ' STU',(STU(K),K=1,3)
      IF(KC.NE.0)ITYPE=CODE(IREF+1,2)
      CALL COEF(STU,CR,ITR1)
      IF(IND.EQ.11)RETURN
      BCR=SQRT(REAL(CR(1)*CONJG(CR(1))+CR(2)*CONJG(CR(2))
     1             +CR(3)*CONJG(CR(3))))
      IF(BCR.LT.1.E-10) THEN
        DO 88 K=1,3
        UC(K)=(0.,0.)
  88    CONTINUE
        GOTO 130
      END IF
  10  CONTINUE
C
C  END OF LOOP OVER INTERFACES
C
C  TERMINATION POINT
C
  100 N1=N2+1
      N2=NN
      IF(KC.NE.0)ITYPE=CODE(IRE,2)
      CALL POLAR(N1,N2,NN,IRE)
C
C COMPUTATION OF CARTESIAN COMPONENTS OF INCIDENT DISPLACEMENT VECTOR
C
      DO 107 K=1,3
      STU(K)=CMPLX(0.,0.)
      DO 107 J=1,3
      STU(K)=HHH(J,K)*CR(J)+STU(K)
 107  CONTINUE
      IF(IAMP.GT.0)WRITE(LOU,'(A,6F10.5)') ' STU',(STU(K),K=1,3)
      DO 105 K=1,6
      Y(K)=AY(K+1,Nn)
      IF(K.LE.3)GO TO 105
      PS(K-3,7,IRE)=Y(K)
  105 CONTINUE
      n=nn
      iref=ire
      IF(MREG.NE.0) THEN
        UC(1)=STU(1)
        UC(2)=STU(2)
        UC(3)=STU(3)
        GOTO 110
      END IF
      IF(IAMP.GT.0)then
        WRITE(LOU,'(a)')
     1  ' REFLECTED SLOWNESS VECTORS AT TERMINATION POINT'
        WRITE(LOU,'(6f12.6)')((PS(L,K,IRE),L=1,3),K=1,3)
      end if
C
C COMPUTATION OF CONVERSION COEFFICIENTS
      KTR=999
      CALL COEF(STU,UC,KTR)
      IF(IND.EQ.11)RETURN
  110 VpEND=1./SQRT(Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6))
      do 120 k=1,3
      y(k)=ay(k+4,1)
  120 continue
      Vp0=1./SQRT(Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3))
      RHO0=0.2*SQRT(AY(8,1))+1.7
      IF(IRHO.NE.0) RHO0=RHO(ISOUR)
      RHEND=0.2*SQRT(AY(8,NN))+1.7
      IF(IRHO.NE.0) RHEND=RHO(LAY)
      AV=AV*Vp0*RHO0
      AV=AV/(VpEND*RHEND)
      UU=SQRT(ABS(AV))
      IF(IAMP.GT.0)
     1WRITE(LOU,'(A,4F12.6)')'VP0,RH0,VPEND,RHEND',Vp0,RHO0,VpEND,RHEND
 130  CONTINUE
      n=nn
      iref=ire
      AMPX(KSS)=UC(1)
      AMPY(KSS)=UC(2)
      AMPZ(KSS)=UC(3)
      IF(ISHEAR.NE.0.AND.KSS.NE.2) THEN
        KSS=2
        ITP=2
        GOTO 3000
      END IF
      RETURN
      END
C
C     *********************************************************
C
      subroutine approx(x,y,yd,kdim)
c
c     The routine performs third-order interpolation between points
c     yold and ynew parameterized by an independent variable x.
c     dold, dnew are the first derivatives of y with respect
c     to x at the points yold and ynew.
c
      dimension y(18),yd(18)
      common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18)
c
      a=(x-xnew)/(xnew-xold)
      aux=a+1.
      a1=(2.*a+3.)*a*a
      a2=1.-a1
      b1=aux*a*(x-xnew)
      b2=aux*a*(x-xold)
      ad1=6.*a*aux/(xnew-xold)
      ad2=-ad1
      bd1=a*(3.*a+2.)
      bd2=aux*(3.*a+1.)
      do 1 i=1,kdim
      y(i)=a1*yold(i)+a2*ynew(i)+b1*dold(i)+b2*dnew(i)
      yd(i)=ad1*yold(i)+ad2*ynew(i)+bd1*dold(i)+bd2*dnew(i)
    1 continue
      return
      end
C
C     *********************************************************
C
      SUBROUTINE BIAP(MX1,MX,MY1,MY,MXY1)
C
      DIMENSION X(200),FX(200),V(1000)
      COMMON/ZCOEF/ A02(1000),A20(1000),A22(1000)
      COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(4),NINT,
     1   XINTA
      EQUIVALENCE(Z(1),V(1))
C
C     ROUTINE DETERMINING THE COEFFICIENTS
C     OF BICUBIC SPLINE INTERPOLATION
C
      DO 1 J=1,MX
      L=MX1+J-1
    1 X(J)=SX(L)
      DO 3 I=1,MY
      DO 2 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    2 FX(J)=V(K)
      CALL SPLIN(X,FX,1,MX)
      DO 3 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    3 A20(K)=FX(J)
C
      DO 4 I=1,MY
      L=MY1+I-1
    4 X(I)=SY(L)
      DO 6 J=1,MX
      DO 5 I=1,MY
      K=MXY1+(J-1)*MY+I-1
    5 FX(I)=V(K)
      CALL SPLIN(X,FX,1,MY)
      DO 6 I=1,MY
      K=MXY1+(J-1)*MY+I-1
    6 A02(K)=FX(I)
C
      DO 7 J=1,MX
      L=MX1+J-1
    7 X(J)=SX(L)
      DO 9 I=1,MY
      DO 8 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    8 FX(J)=A02(K)
      CALL SPLIN(X,FX,1,MX)
      DO 9 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    9 A22(K)=FX(J)
C
      RETURN
      END
C
C
C     *********************************************************
C
      SUBROUTINE CHRM(Y)
C
C  ROUTINE FOR THE COMPUTATION OF THE ELEMENTS OF THE CHRISTOFFEL
C  MATRIX FOR AN ARBITRARY ANISOTROPIC MEDIUM
C
      DIMENSION Y(18)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DJK/  D11,D12,D13,D22,D23,D33,DTR
      COMMON /GAM/  C11,C12,C13,C22,C23,C33
C
      P1=Y(4)
      P2=Y(5)
      P3=Y(6)
      P2P3=P2*P3
      P1P2=P1*P2
      P1P3=P1*P3
      P1P1=P1*P1
      P2P2=P2*P2
      P3P3=P3*P3
      C11=P1P1*A11+P2P2*A66+P3P3*A55
     1+2.*(P2P3*A56+P1P3*A15+P1P2*A16)
      C22=P1P1*A66+P2P2*A22+P3P3*A44
     1+2.*(P2P3*A24+P1P3*A46+P1P2*A26)
      C33=P1P1*A55+P2P2*A44+P3P3*A33
     1+2.*(P2P3*A34+P1P3*A35+P1P2*A45)
      C23=P1P1*A56+P2P2*A24+P3P3*A34
     1   +P2P3*A2344+P1P3*A3645+P1P2*A2546
      C13=P1P1*A15+P2P2*A46+P3P3*A35
     1   +P2P3*A3645+P1P3*A1355+P1P2*A1456
      C12=P1P1*A16+P2P2*A26+P3P3*A45
     1   +P2P3*A2546+P1P3*A1456+P1P2*A1266
      C11N=C11-1.
      C22N=C22-1.
      C33N=C33-1.
      C23SQ=C23*C23
      C13SQ=C13*C13
      C12SQ=C12*C12
      D11=C22N*C33N-C23SQ
      D22=C11N*C33N-C13SQ
      D33=C11N*C22N-C12SQ
      D12=C13*C23-C12*C33N
      D13=C12*C23-C13*C22N
      D23=C12*C13-C23*C11N
      DTR=D11+D22+D33
      IF(ABS(DTR).LT.0.00001)THEN
        IND=10
      END IF
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE CHRM1(C,PN,UN)
C
C  ROUTINE FOR THE COMPUTATION OF THE ELEMENTS OF THE CHRISTOFFEL
C  MATRIX FOR AN ARBITRARY ANISOTROPIC MEDIUM
C
      DIMENSION C(3,3),PN(3),UN(3)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
C
      P1=PN(1)
      P2=PN(2)
      P3=PN(3)
      U1=UN(1)
      U2=UN(2)
      U3=UN(3)
      P2U3=P2*U3
      P3U2=P3*U2
      P1U2=P1*U2
      P2U1=P2*U1
      P1U3=P1*U3
      P3U1=P3*U1
      P1U1=P1*U1
      P2U2=P2*U2
      P3U3=P3*U3
      C(1,1)=P1U1*A11+P2U2*A66+P3U3*A55
     1+(P2U3+P3U2)*A56+(P1U3+P3U1)*A15+(P1U2+P2U1)*A16
      C(2,2)=P1U1*A66+P2U2*A22+P3U3*A44
     1+(P2U3+P3U2)*A24+(P1U3+P3U1)*A46+(P1U2+P2U1)*A26
      C(3,3)=P1U1*A55+P2U2*A44+P3U3*A33
     1+(P2U3+P3U2)*A34+(P1U3+P3U1)*A35+(P1U2+P2U1)*A45
      C(2,3)=P1U1*A56+P2U2*A24+P3U3*A34
     1+0.5*((P2U3+P3U2)*A2344+(P1U3+P3U1)*A3645+(P1U2+P2U1)*A2546)
      C(1,3)=P1U1*A15+P2U2*A46+P3U3*A35
     1+0.5*((P2U3+P3U2)*A3645+(P1U3+P3U1)*A1355+(P1U2+P2U1)*A1456)
      C(1,2)=P1U1*A16+P2U2*A26+P3U3*A45
     1+0.5*((P2U3+P3U2)*A2546+(P1U3+P3U1)*A1456+(P1U2+P2U1)*A1266)
      C(2,1)=C(1,2)
      C(3,2)=C(2,3)
      C(3,1)=C(1,3)
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE CHRM2(Y,G,i)
C
C     EVALUATES ELEMENTS OF THE CHRISTOFFEL MATRIX
C
      DIMENSION a(21),Y(18),G(3,3)
      COMMON/GAM/G11,G12,G13,G22,G23,G33
      COMMON /APROX1/ e(21,10)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
C
      do 1 j=1,21
      a(j)=e(j,i)
    1 continue
      P1=Y(4)
      P2=Y(5)
      P3=Y(6)
      P11=P1*P1
      P12=P1*P2
      P13=P1*P3
      P22=P2*P2
      P23=P2*P3
      P33=P3*P3
      G11=A(1)*P11+A(21)*P22+A(19)*P33+
     1    2.*(A(6)*P12+A(5)*P13+A(20)*P23)
      G22=A(21)*P11+A(7)*P22+A(16)*P33+
     1    2.*(A(11)*P12+A(18)*P13+A(9)*P23)
      G33=A(19)*P11+A(16)*P22+A(12)*P33+
     1    2.*(A(17)*P12+A(14)*P13+A(13)*P23)
      G12=A(6)*P11+A(11)*P22+A(17)*P33+
     1    (A(21)+A(2))*P12+(A(20)+A(4))*P13+(A(10)+A(18))*P23
      G13=A(5)*P11+A(18)*P22+A(14)*P33+
     1    (A(20)+A(4))*P12+(A(19)+A(3))*P13+(A(17)+A(15))*P23
      G23=A(20)*P11+A(9)*P22+A(13)*P33+
     1    (A(10)+A(18))*P12+(A(17)+A(15))*P13+(A(16)+A(8))*P23
      G(1,1)=G11
      G(1,2)=G12
      G(1,3)=G13
      G(2,1)=G12
      G(2,2)=G22
      G(2,3)=G23
      G(3,1)=G13
      G(3,2)=G23
      G(3,3)=G33
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE PCHRM(Y,G,L,i)
C
C     EVALUATES FIRST DERIVATIVES OF ELEMENTS OF CHRISTOFFEL MATRIX
C     WITH RESPECT TO THE L-TH COMPONENT OF THE SLOWNESS VECTOR
C
      DIMENSION a(21),Y(18),G(3,3)
      COMMON /APROX1/ e(21,10)
C
      do 1 j=1,21
      a(j)=e(j,i)
    1 continue
      P1=Y(4)
      P2=Y(5)
      P3=Y(6)
      IF(L.EQ.1)THEN
        G(1,1)=2.*(A(1)*P1+A(6)*P2+A(5)*P3)
        G(2,2)=2.*(A(21)*P1+A(11)*P2+A(18)*P3)
        G(3,3)=2.*(A(19)*P1+A(17)*P2+A(14)*P3)
        AUX=2.*A(6)*P1+(A(21)+A(2))*P2+(A(20)+A(4))*P3
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(5)*P1+(A(20)+A(4))*P2+(A(19)+A(3))*P3
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(20)*P1+(A(10)+A(18))*P2+(A(17)+A(15))*P3
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.2)THEN
        G(1,1)=2.*(A(6)*P1+A(21)*P2+A(20)*P3)
        G(2,2)=2.*(A(11)*P1+A(7)*P2+A(9)*P3)
        G(3,3)=2.*(A(17)*P1+A(16)*P2+A(13)*P3)
        AUX=2.*A(11)*P2+(A(21)+A(2))*P1+(A(10)+A(18))*P3
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(18)*P2+(A(20)+A(4))*P1+(A(17)+A(15))*P3
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(9)*P2+(A(10)+A(18))*P1+(A(16)+A(8))*P3
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.3)THEN
        G(1,1)=2.*(A(5)*P1+A(20)*P2+A(19)*P3)
        G(2,2)=2.*(A(18)*P1+A(9)*P2+A(16)*P3)
        G(3,3)=2.*(A(14)*P1+A(13)*P2+A(12)*P3)
        AUX=2.*A(17)*P3+(A(20)+A(4))*P1+(A(10)+A(18))*P2
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(14)*P3+(A(19)+A(3))*P1+(A(17)+A(15))*P2
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(13)*P3+(A(17)+A(15))*P1+(A(16)+A(8))*P2
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE PPCHRM(G,L,M,i)
C
C     EVALUATES SECOND DERIVATIVES OF ELEMENTS OF CHRISTOFFEL MATRIX
C     WITH RESPECT TO THE L-TH AND M-TH COMPONENTS OF THE SLOWNESS
C     VECTOR
C
      DIMENSION a(21),G(3,3)
      COMMON /APROX1/ e(21,10)
C
      do 1 j=1,21
      a(j)=e(j,i)
    1 continue
      IF(L.EQ.1.AND.M.EQ.1)THEN
        G(1,1)=2.*A(1)
        G(2,2)=2.*A(21)
        G(3,3)=2.*A(19)
        AUX=2.*A(6)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(5)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(20)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.2.AND.M.EQ.2)THEN
        G(1,1)=2.*A(21)
        G(2,2)=2.*A(7)
        G(3,3)=2.*A(16)
        AUX=2.*A(11)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(18)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(9)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.3.AND.M.EQ.3)THEN
        G(1,1)=2.*A(19)
        G(2,2)=2.*A(16)
        G(3,3)=2.*A(12)
        AUX=2.*A(17)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(14)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(13)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF((L.EQ.1.AND.M.EQ.2).OR.(L.EQ.2.AND.M.EQ.1))THEN
        G(1,1)=2.*A(6)
        G(2,2)=2.*A(11)
        G(3,3)=2.*A(17)
        AUX=A(21)+A(2)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=A(20)+A(4)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=A(10)+A(18)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF((L.EQ.1.AND.M.EQ.3).OR.(L.EQ.3.AND.M.EQ.1))THEN
        G(1,1)=2.*A(5)
        G(2,2)=2.*A(18)
        G(3,3)=2.*A(14)
        AUX=A(20)+A(4)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=A(19)+A(3)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=A(17)+A(15)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF((L.EQ.2.AND.M.EQ.3).OR.(L.EQ.3.AND.M.EQ.2))THEN
        G(1,1)=2.*A(20)
        G(2,2)=2.*A(9)
        G(3,3)=2.*A(13)
        AUX=A(10)+A(18)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=A(17)+A(15)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=A(16)+A(8)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'a2.for'
C     <A HREF="a2.for">a2.for</A>
      INCLUDE 'a3.for'
C     <A HREF="a3.for">a3.for</A>
      INCLUDE 'a4.for'
C     <A HREF="a4.for">a4.for</A>
      INCLUDE 'a5.for'
C     <A HREF="a5.for">a5.for</A>
C
C <A NAME="MOD"></A><FONT COLOR="RED">Interpolation method:</FONT>
C Include just one of the following files 'mod*.for':
C (a) Isosurface interpolation:
*     INCLUDE 'modis.for'
C     <A HREF="modis.for">modis.for</A>
C (b) (Bi-)(tri-)cubic B-spline interpolation:
*     INCLUDE 'modbs.for'
C     <A HREF="modbs.for">modbs.for</A>
C
C=======================================================================
C                                                                 </PRE>
