      SUBROUTINE RAYA (X0,Y0,Z0,T0,FI0,PSI0,PX,PY,PZ,XX,YY,ZZ,T,DT,AC)
C
C  INITIAL-VALUE RAY TRACING BY THE RUNGE-KUTTA METHOD
C
      DIMENSION Y(18),DEP(6),PRM(5),DERY(18),AUX(8,18),DIN(18),VSQ(3)
      dimension pn(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
      COMMON /AUXI/  IANI(20),INTR,INT1,IOUT,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 /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),brd(4),NINT
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),DINC,DTOLD,N,IREF,IND,IND1
      COMMON /ZERO/ RNULL
      common/dyn/xdyn(3,3),ydyn(3,3)
      common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18)
      COMMON/KM/KMAH,SPROLD,TSTOLD
      EXTERNAL FCT,OUT
C
      IAC=0
      kmah=0
      Y(1)=X0
      Y(2)=Y0
      Y(3)=Z0
      IREFR=0
      KRE=KREF
      IF(KC.EQ.0) KRE=0
      N=0
      IREF=1
      IOUT=0
      IF(IND.GT.0) GOTO 6
C
C  FOR IND=-1 DETERMINATION OF THE NUMBER OF THE LAYER IN WHICH THE
C  INITIAL POINT OF THE RAY IS SITIUATED
C
      IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2).OR.Y(2).LT.BRD(3).OR.Y(2).GT.
     1   BRD(4)) GOTO 4
      INTR=1
   1  CALL DISC (Y,DEP)
      ZINT=DEP(1)
      IF(ABS(Y(3)).GT.0.00001) GOTO 2
      ISOUR=1
C
C  REDEFINITION OF Z-COORDINATES FOR A SOURCE ON THE EARTH SURFACE
C
      Z0=ZINT
      GOTO 5
   2  IF(Y(3).LT.ZINT.AND.INTR.EQ.1) GOTO 4
      IF(Y(3).LT.ZINT) GOTO 5
      IF(ABS(Y(3)-ZINT).LT.RNULL.AND.INTR.EQ.NINT) GOTO 5
      IF(INTR.EQ.NINT) GOTO 4
      ISOUR=INTR
      INTR=INTR+1
      GOTO 1
   4  WRITE(lou,'(A,/,6F10.5,/,4F10.5)')' Y,BRD',Y,BRD
      IND=50
      RETURN
C
C  DETERMINATION OF INITIAL CONDITIONS FOR THE RUNGE-KUTTA PROCEDURE
C
   5  IF(IND.GE.0) GOTO 6
      IND=ISOUR
      RETURN
   6  LAY=ISOUR
      is(3,1)=lay
      INT1=ISOUR
      IF(ISOUR.NE.CODE(1,1)) IND=14
      IF(ISOUR.NE.CODE(1,1)) RETURN
      ITYPE=CODE(1,2)
      IF(IANI(ISOUR).EQ.0.AND.ITYPE.EQ.2.and.kc.ne.0) IND=38
      IF(IANI(ISOUR).EQ.0.AND.ITYPE.EQ.2.and.kc.ne.0) RETURN
      CALL PARDIS(Y,0)
C
C  DETERMINATION OF INITIAL VALUES FOR RAY TRACING 
C  AND DYNAMIC RAY TRACING
C
      csp=cos(psi0)
      snp=sin(psi0)
      csf=cos(fi0)
      snf=sin(fi0)
      if(mori.eq.0)then
        Y(4)=CSP*CSF
        Y(5)=SNP*CSF
        Y(6)=SNF
      else
        y(4)=csf*csp
        y(5)=snp
        y(6)=snf*csp
      end if
      do 8 k=1,3
        pn(k)=y(k+3)
    8 continue
      if(nder.gt.1)then
        do 3 k=7,12
          y(k)=0.
    3   continue
        if(mori.eq.0)then
 	    y(13)=-snp*csf
          y(14)=csp*csf
	    y(15)=0.
	    y(16)=-csp*snf
	    y(17)=-snp*snf
	    y(18)=csf
        else
          y(13)=-snf*csp
   	    y(14)=0.
          y(15)=csf*csp
	    y(16)=-csf*snp
	    y(17)=csp
	    y(18)=-snf*snp
        end if
      end if
      IF(IANI(ISOUR).ne.0)then
C
C  SOURCE LOCATED IN AN ANISOTROPIC LAYER
C
        CALL INIT(pn,VSQ)
        IF(IPRINT.GT.2)WRITE(lou,'(a,3F14.6)')' V1,V2,V3=', VSQ
        VP=AMAX1(VSQ(1),VSQ(2),VSQ(3))
        VS1=AMIN1(VSQ(1),VSQ(2),VSQ(3))
        VS2=VSQ(1)+VSQ(2)+VSQ(3)-VP-VS1
        IF(ITYPE.EQ.3)V=SQRT(VP)
        IF(ITYPE.EQ.1)V=SQRT(VS1)
        IF(ITYPE.EQ.2)V=SQRT(VS2)
        do 7 i=4,6
          y(i)=y(i)/v
    7   continue
        if(nder.gt.1)then
          nder=1
          call fct(0.,y,dery)
          nder=2
          vg=sqrt(dery(1)*dery(1)+dery(2)*dery(2)+dery(3)*dery(3))
          if(mori.eq.0)then
            auxp=-dery(1)*snp*csf+dery(2)*csp*csf
            auxf=-dery(1)*csp*snf-dery(2)*snp*snf+dery(3)*csf
            y(13)=y(13)-auxp*csp*csf/v
            y(14)=y(14)-auxp*snp*csf/v
            y(15)=y(15)-auxp*snf/v
            y(16)=y(16)-auxf*csp*csf/v
            y(17)=y(17)-auxf*snp*csf/v
            y(18)=y(18)-auxf*snf/v
          else
            auxp=-dery(1)*snf*csp+dery(3)*csf*csp
            auxf=-dery(1)*csf*snp+dery(2)*csp-dery(3)*snf*snp
            y(13)=y(13)-auxp*csf*csp/v
            y(14)=y(14)-auxp*snp/v
            y(15)=y(15)-auxp*snf*csp/v
            y(16)=y(16)-auxf*csf*csp/v
            y(17)=y(17)-auxf*snp/v
            y(18)=y(18)-auxf*snf*csp/v
          end if
          do 11 i=13,18
            y(i)=y(i)/v
   11     continue
C
C     DETERMINATION OF KMAH INDEX AT SOURCE IN ANISTROPIC MEDIUM
C
          call fct(0.,y,dery)
          if(mori.eq.0)then
            el=-(-dery(7)*snp*csf+dery(8)*csp*csf)/vg
            em=-(-dery(10)*snp*csf+dery(11)*csp*csf)/vg
            en=-(-dery(10)*csp*snf-dery(11)*snp*snf+dery(12)*csf)/vg
          else
            el=-(-dery(7)*snf*csp+dery(9)*csf*csp)/vg
            em=-(-dery(10)*snf*csp+dery(12)*csf*csp)/vg
            en=-(-dery(10)*csf*snp+dery(11)*csp-dery(12)*snf*snp)/vg
          end if
          ee=y(13)*y(13)+y(14)*y(14)+y(15)*y(15)        
          ff=y(13)*y(16)+y(14)*y(17)+y(15)*y(18)        
          gg=y(16)*y(16)+y(17)*y(17)+y(18)*y(18)        
          egf=(ee*gg-ff*ff)/v/v
          elnm=el*en-em*em
          be=el*gg+en*ee-2.*em*ff
          if(abs(egf).lt..000001)egf=sign(1.,egf)*.00001
          gk=elnm/egf
          ak=.5*be/egf
          if(gk.gt.0.)then
            if(ak.lt.0.)kmah=0
            if(ak.gt.0.)kmah=2
          end if
          if(gk.lt.0.)kmah=1
        end if
      end if
C
C  SOURCE LOCATED IN AN ISOTROPIC LAYER
C
      IF(IANI(ISOUR).eq.0)then
        IF(ITYPE.EQ.3)V=SQRT(A11)
        IF(ITYPE.NE.3)V=SQRT(A44)
        kdim=6
        if(nder.gt.1)kdim=18
        do 9 i=4,kdim
            y(i)=y(i)/v
    9   continue
      end if
C
      IND=0
      IND1=0
      PRM(1)=T0
      PRM(2)=TMAX
      PRM(3)=DT
      IF(ITYPE.NE.3) PRM(3)=DT*1.7
      PRM(4)=AC
      T=PRM(1)
   20 DO 10 I=1,3
        auxx=y(4)*y(4)+y(5)*y(5)+y(6)*y(6)
        auxx=sqrt(auxx)
        DIN(I)=auxx
        din(i+3)=prm(3)/auxx
   10 CONTINUE
      kdim=6
      if(nder.gt.1)kdim=18
      do 25 i=4,kdim
        din(i)=0.
   25 continue
      DO 30 I=1,kdim
        DERY(I)=DIN(I)
   30 CONTINUE
C
C  COMPUTATION OF THE RAY
C
      CALL RKGS(PRM,Y,DERY,kdim,IHLF,FCT,OUT,AUX)
      IF(IHLF.EQ.11) IND=5
      IF(IHLF.EQ.12) IND=6
      IF(IHLF.EQ.13) IND=7
      IF(IND.GE.5.AND.IND.LE.7) RETURN
      IF(ABS(PRM(5)).GT.0.0001) GOTO 35
      IF(IND.eq.12) GOTO 70
      GOTO 60
   35 IF(ABS(PRM(5)-1.1).GT.0.0001) GOTO 50
C
C  LOOP FOR THE ITERATIVE DETERMINATION OF THE POINT OF INCIDENCE
C
      PRM(1)=AY(1,N)
      DO 40 I=1,6
        Y(I)=Yold(I)
  40  CONTINUE
      if(nder.gt.1)then
	do 37 i=1,12
	  y(i+6)=yold(i+6)
  37    continue
      end if
      T=AY(1,N)
      N=N-1
      GOTO 20
  50  XX=Y(1)
      YY=Y(2)
      ZZ=Y(3)
      T=AY(1,N)
      IF(ABS(PRM(5)-2.).GT.0.0001) GOTO 80
C
C  INTEGRATION FROM THE POINT OF REFLECTION/TRANSMISSION TO THE CLOSEST
C  POINT OF THE RAY CORRESPONDING TO REGULAR TIME STEP
C
      PRM(1)=AY(1,N)
      I=INT((PRM(1)-T0)/DT)
      PRM(2)=FLOAT(I+1)*DT+T0
      PRM(3)=DT
      GOTO 20
   60 PRM(1)=PRM(2)
      PRM(2)=TMAX
      PRM(3)=DT
      IF(ITYPE.NE.3) PRM(3)=1.7*DT
      N=N-1
      GOTO 20
   70 CONTINUE
      XX=Y(1)
      YY=Y(2)
      ZZ=Y(3)
   80 continue
c      if(kmah.ne.0)ind1=ind1+50
      IF(IREFR.EQ.1) IND1=-IND1
      PX=Y(4)
      PY=Y(5)
      PZ=Y(6)
      if(nder.gt.1)then
        do 90 i=1,3
          xdyn(i,1)=y(i+6)
          xdyn(i,2)=y(i+9)
          xdyn(i,3)=dery(i)
          ydyn(i,1)=y(i+12)
          ydyn(i,2)=y(i+15)
          ydyn(i,3)=dery(i+3)
   90   continue
      end if
      RETURN
      END
