      SUBROUTINE RAYB(W,T0,DT)
C
C  INITIAL-VALUE RAY TRACING BY THE RUNGE-KUTTA METHOD
C
      COMPLEX AY,W(2)
      COMMON /FORCE/ F(3)
      COMMON /RAYW/ AY(3,2000),E(3,3,2000),OMEGA,N,NTOT,IND,IND1
      COMMON /ISOTR/T(2000),XX(2000),Y(2000),Z(2000),VP(2000),VS(2000),
     1RHO(2000)
C
      X=T0
      PI=3.14159632
      RHO0=RHO(1)
      VS0=VS(1)
      AUX=4*PI*SQRT(RHO0)*RHO0*VS0*VS0
      W(1)=(E(1,1,1)*F(1)+E(1,2,1)*F(2)+E(1,3,1)*F(3))/AUX
      W(2)=(E(2,1,1)*F(1)+E(2,2,1)*F(2)+E(2,3,1)*F(3))/AUX
      AY(1,1)=T0
      AY(2,1)=W(1)
      AY(3,1)=W(2)
C
C  SOLVING OF THE SYSTEM OF COUPLED EQUATIONS
C
      CALL RK(X,W)
      RETURN
      END
C
C     ***********************************************************
C
      SUBROUTINE RK(X,Y)
C
C     MODIFIED EULER'S METHOD TO SOLVE A SYSTEM OF LINEAR
C     ORDINARY DIFFERENTIAL EQUATIONS OF FIRST ORDER
C
      COMPLEX AY,Y(2),DERY(2),Y1(2),Y2(2)
      COMMON /RAYW/ AY(3,2000),HHH(3,3,2000),OMEGA,N,NTOT,IND,IND1
      COMMON /ISOTR/T(2000),XX(2000),YY(2000),ZZ(2000),VP(2000),
     1VS(2000),RHO(2000)
C
      NDIM=2
      N=1
    1 H=T(N+1)-T(N)
      X=X+H
      CALL FCTA(Y,DERY,NDIM)
      DO 2 I=1,NDIM
      Y1(I)=DERY(I)
    2 Y2(I)=Y(I)+H*Y1(I)
      N=N+1
      CALL FCTA(Y2,DERY,NDIM)
      DO 3 I=1,NDIM
    3 Y(I)=Y(I)+.5*H*(Y1(I)+DERY(I))
      CALL OUT(X,Y,DERY,NDIM)
      IF(N.EQ.NTOT) GO TO 4
      GO TO 1
    4 RETURN
      END
C
C********************************************************
C
      SUBROUTINE FCTA(W,DERY,NDIM)
C
C  COMPUTATION OF THE RIGHT HAND SIDES OF THE QI EQUATIONS
C
      SAVE VS2
      COMPLEX AY,DERY(NDIM),W(NDIM),IMAG
      DIMENSION YY(3)
      COMMON /RAYW/ AY(3,2000),E(3,3,2000),OMEGA,N,NTOT,IND,IND1
      COMMON /ISOTR/T(2000),X(2000),Y(2000),Z(2000),VP(2000),VS(2000),
     1RHO(2000)
C
      IMAG=(0.,1.)
      VS2=VS(N)*VS(N)
      YY(1)=X(N)
      YY(2)=Y(N)
      YY(3)=Z(N)
      CALL PARD(YY)
      CALL BIJ(B11,1,1)
      CALL BIJ(B12,1,2)
      CALL BIJ(B22,2,2)
C     
C     MODIFICATION OF WEAK ANISOTROPY MATRIX BY ROUTINE NEWB
C
      IF(IND1.EQ.1)THEN
        VSS=VS(N)
        CALL NEWB(VSS,B11,B12,B22)
      END IF
C
      DERY(1)=-0.5*IMAG*(OMEGA/VS2)*(W(1)*B11+W(2)*B12)
      DERY(2)=-0.5*IMAG*(OMEGA/VS2)*(W(1)*B12+W(2)*B22)
      RETURN
      END
C
C****************************************************
C
      SUBROUTINE OUT(X,Y,DERY,NDIM)
C
      COMPLEX AY,Y(2),DERY(2)
      COMMON /RAYW/ AY(3,2000),E(3,3,2000),OMEGA,N,NTOT,IND,IND1
C
      AY(1,N)=X
      AY(2,N)=Y(1)
      AY(3,N)=Y(2)
      RETURN
      END
C
C*********************************************************
C
      SUBROUTINE BIJ(B,M,N)
C
      COMPLEX AY
      COMMON /RAYW/   AY(3,2000),E(3,3,2000),OMEGA,NN,NTOT,IND,IND1
      COMMON /ELPAR/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66
      COMMON /ISOTR/T(2000),X(2000),Y(2000),Z(2000),VP(2000),VS(2000),
     1RHO(2000)
C
      A2344=A23+A44
      A3645=A36+A45
      A2546=A25+A46
      A1355=A13+A55
      A1456=A14+A56
      A1266=A12+A66
      T1=E(3,1,NN)
      T2=E(3,2,NN)
      T3=E(3,3,NN)
      EM1=E(M,1,NN)
      EM2=E(M,2,NN)
      EM3=E(M,3,NN)
      EN1=E(N,1,NN)
      EN2=E(N,2,NN)
      EN3=E(N,3,NN)
      T2T3=T2*T3
      T1T2=T1*T2
      T1T3=T1*T3
      T1T1=T1*T1
      T2T2=T2*T2
      T3T3=T3*T3
      C11=T1T1*A11+T2T2*A66+T3T3*A55
     1+2.*(T2T3*A56+T1T3*A15+T1T2*A16)
      C22=T1T1*A66+T2T2*A22+T3T3*A44
     1+2.*(T2T3*A24+T1T3*A46+T1T2*A26)
      C33=T1T1*A55+T2T2*A44+T3T3*A33
     1+2.*(T2T3*A34+T1T3*A35+T1T2*A45)
      C23=T1T1*A56+T2T2*A24+T3T3*A34
     1   +T2T3*A2344+T1T3*A3645+T1T2*A2546
      C13=T1T1*A15+T2T2*A46+T3T3*A35
     1   +T2T3*A3645+T1T3*A1355+T1T2*A1456
      C12=T1T1*A16+T2T2*A26+T3T3*A45
     1   +T2T3*A2546+T1T3*A1456+T1T2*A1266
      B=C11*EM1*EN1+C22*EM2*EN2+C33*EM3*EN3+C12*(EM1*EN2+EM2*EN1)
     1   +C13*(EM1*EN3+EM3*EN1)+C23*(EM2*EN3+EM3*EN2)
      IF(M.EQ.N)THEN
        V2=VS(NN)*VS(NN)
        IF(M.EQ.3)V2=VP(NN)*VP(NN)
        B=B-V2
      END IF
      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(6),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
      SUBROUTINE DISC(Y,DEP)
C
C     DETERMINATION OF DEPTH OF 3D INTERFACES AND ITS DERIVATIVES
C     FOR BICUBIC POLYNOMIAL APPROXIMATION
C
      SAVE B,DX,DY,llay,ibb,iu,il,ju,jl
      DIMENSION Y(18),DEP(6),B(16,2),DX(2),DY(2)
C
      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 /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT,
     1   XINTA
      COMMON/ZCOEF/ A02(1000),A20(1000),A22(1000)
C
      IBB=0
      LB=2
      IF(INTR.EQ.LAY)LB=1
      MX=NX(INTR)
      MY=NY(INTR)
      DO 1 I=2,MX
      K=MMX(INTR)+I-1
      IF(Y(1).LT.SX(K))GO TO 2
    1 CONTINUE
    2 I1=K
      DO 3 I=2,MY
      K=MMY(INTR)+I-1
      IF(Y(2).LT.SY(K))GO TO 4
    3 CONTINUE
    4 J1=K
      IF(MAUX.EQ.0) GOTO 8
      IF(LB.EQ.2) GOTO 5
      IF(I1.EQ.IU.AND.J1.EQ.JU.AND.LLAY.EQ.LAY) GOTO 10
      GOTO 8
    5 IF(I1.EQ.IL.AND.J1.EQ.JL.AND.LLAY.EQ.LAY) GOTO 10
      IL=I1
      JL=J1
      GOTO 9
    8 IU=I1
      JU=J1
    9 LLAY=LAY
      I=I1-MMX(INTR)
      J=J1-MMY(INTR)
      DX(LB)=SX(I1-1)
      DY(LB)=SY(J1-1)
      MM=MMXY(INTR)-1
      K=MM+(I-1)*MY+J
      B20=A20(K)
      B02=A02(K)
      B22=A22(K)
      B00=Z(K)
      K=MM+I*MY+J
      C20=A20(K)
      C02=A02(K)
      C22=A22(K)
      C00=Z(K)
      K=MM+(I-1)*MY+J+1
      D20=A20(K)
      D02=A02(K)
      D22=A22(K)
      D00=Z(K)
      K=MM+I*MY+J+1
      E20=A20(K)
      E02=A02(K)
      E22=A22(K)
      E00=Z(K)
      HX=SX(I1)-DX(LB)
      HY=SY(J1)-DY(LB)
      XA=3.*HX
      YA=3.*HY
      XB=HX/3.
      YB=HY/3.
      D32=(E22-D22)/XA
      D30=(E20-D20)/XA
      B30=(C20-B20)/XA
      B32=(C22-B22)/XA
      D12=(E02-D02)/HX-XB*(E22+2.*D22)
      D10=(E00-D00)/HX-XB*(E20+2.*D20)
      B10=(C00-B00)/HX-XB*(C20+2.*B20)
      B12=(C02-B02)/HX-XB*(C22+2.*B22)
      B03=(D02-B02)/YA
      B13=(D12-B12)/YA
      B23=(D22-B22)/YA
      B33=(D32-B32)/YA
      B01=(D00-B00)/HY-YB*(D02+2.*B02)
      B11=(D10-B10)/HY-YB*(D12+2.*B12)
      B21=(D20-B20)/HY-YB*(D22+2.*B22)
      B31=(D30-B30)/HY-YB*(D32+2.*B32)
      MAUX=1
      B(1,LB)=B00
      B(2,LB)=B01
      B(3,LB)=B02
      B(4,LB)=B03
      B(5,LB)=B10
      B(6,LB)=B11
      B(7,LB)=B12
      B(8,LB)=B13
      B(9,LB)=B20
      B(10,LB)=B21
      B(11,LB)=B22
      B(12,LB)=B23
      B(13,LB)=B30
      B(14,LB)=B31
      B(15,LB)=B32
      B(16,LB)=B33
      IBB=1
   10 AX=Y(1)-DX(LB)
      AZ=Y(2)-DY(LB)
      IF(IBB.EQ.1) GOTO 11
      B00=B(1,LB)
      B01=B(2,LB)
      B02=B(3,LB)
      B03=B(4,LB)
      B10=B(5,LB)
      B11=B(6,LB)
      B12=B(7,LB)
      B13=B(8,LB)
      B20=B(9,LB)
      B21=B(10,LB)
      B22=B(11,LB)
      B23=B(12,LB)
      B30=B(13,LB)
      B31=B(14,LB)
      B32=B(15,LB)
      B33=B(16,LB)
   11 AUX1=((B33*AZ+B32)*AZ+B31)*AZ+B30
      AUX2=((B23*AZ+B22)*AZ+B21)*AZ+B20
      AUX3=((B13*AZ+B12)*AZ+B11)*AZ+B10
      AUX4=((B03*AZ+B02)*AZ+B01)*AZ+B00
      DEP(1)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4
      IF(NDER.EQ.0) RETURN
      DEP(2)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3
      IF(NDER.EQ.1)GO TO 7
      DEP(4)=6.*AUX1*AX+2.*AUX2
    7 AUX1=(3.*B33*AZ+2.*B32)*AZ+B31
      AUX2=(3.*B23*AZ+2.*B22)*AZ+B21
      AUX3=(3.*B13*AZ+2.*B12)*AZ+B11
      AUX4=(3.*B03*AZ+2.*B02)*AZ+B01
      DEP(3)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4
      IF(NDER.EQ.1) RETURN
      DEP(5)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3
      AUX1=3.*B33*AZ+B32
      AUX2=3.*B23*AZ+B22
      AUX3=3.*B13*AZ+B12
      AUX4=3.*B03*AZ+B02
      DEP(6)=2.*(((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4)
      RETURN
      END
C
C*********************************************************
C
      SUBROUTINE SPLIN(X,FX,NMIN,NMAX)
C
      DIMENSION A(200),B(200),H(200),F(200),X(200),FX(200)
C
C     CUBIC SPLINE INTERPOLATION WITH ZERO CURVATURES AT
C     THE END POINTS
C
      IF((NMAX-NMIN).EQ.1)GO TO 4
      NMIN1=NMIN+1
      NMAX1=NMAX-1
      H(NMIN1)=X(NMIN1)-X(NMIN)
      D2=(FX(NMIN1)-FX(NMIN))/H(NMIN1)
      DO 1 I=NMIN1,NMAX1
      H(I+1)=X(I+1)-X(I)
      D1=D2
      D2=(FX(I+1)-FX(I))/H(I+1)
      B(I)=H(I)+H(I+1)
      FX(I)=3.*(D2-D1)/B(I)
      A(I)=H(I)/B(I)
    1 B(I)=H(I+1)/B(I)
    4 FX(NMIN)=0.
      FX(NMAX)=0.
      IF((NMAX-NMIN).EQ.1)RETURN
      H(NMIN)=0.
      F(NMIN)=0.
      DO 2 I=NMIN1,NMAX1
      XPOM=2.+A(I)*H(I-1)
      H(I)=-B(I)/XPOM
    2 F(I)=(FX(I)-A(I)*F(I-1))/XPOM
      DO 3 I=NMIN,NMAX1
      J=NMAX1-(I-NMIN)
    3 FX(J)=H(J)*FX(J+1)+F(J)
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE FACETS(N1,N2,NSRF)
      INTEGER LU,N1,N2,NSRF
C
C Subroutine FACETS writes the index file listing the vertices of each
C tetragon covering the structural interface.  The vertices are assumed
C to be stored in a separate file, with inner loop over N1 points along
C the first horizontal axis, middle loop over N2 points along the second
C horizontal axis and outer loop over the surfaces.  The vertices are
C indexed by positive integers according to their order in the vertex
C file.
C
C Input:
C     LU...   Logical unit number connected to the output file to be
C             written by this subroutine.
C     N1...   Number of points along the first horizontal axis.
C     N2...   Number of points along the second horizontal axis.
C     NSRF... Number of interfaces.
C The input parameters are not altered.
C
C No output.
C
C Output index file with the tetragons:
C For each tetragon, a line containing I1,I2,I3,I4,/
C     I1,I2,I3,I4... Indices of the vertices of the tetragon.
C             The vertices are indexed by positive integers according to
C             their order in the respective vertex file.
C     /...    List of vertices is terminated by a slash.
C
C Date: 1999, October 4
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      CHARACTER*9 FORMAT
      INTEGER I1,I2,ISRF
      COMMON/VRML/LUBRD,LUGRD,LU,LURAY
C
      IF(LU.EQ.0)RETURN
C     Setting output format:
      FORMAT='(4(I0,A))'
      I1=INT(ALOG10(FLOAT(N1*N2*NSRF)+0.5))+1
      FORMAT(5:5)=CHAR(ICHAR('0')+I1)
C
C     Writing the file:
      DO 33 ISRF=0,N1*N2*(NSRF-1),N1*N2
        DO 32 I2=ISRF,ISRF+N1*(N2-2),N1
          DO 31 I1=I2+1,I2+N1-1
            WRITE(LU,FORMAT) I1,' ',I1+1,' ',I1+1+N1,' ',I1+N1,' /'
   31     CONTINUE
   32   CONTINUE
   33 CONTINUE
C
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE NEWB(VS,B11,B12,B22)
      REAL VS,B11,B12,B22
C
C Subroutine to slightly modify matrix (7) of Psencik (1998):
C Quasi-shear waves in the zero-order approximation of the
C quasi-isotropic approach.  Preliminary results.  In: Seismic Waves
C in Complex 3-D Structures, Report 7, pp.225-266.
C The modification makes the coupling equations dependent on the
C reference isotropic model just in terms of the reference ray path.
C
C Input:
C     VS...   S-wave velocity in the reference isotropic model.
C     B11,B12,B22... Components of symmetric 2*2 matrix B,
C             B(M,N) = DeltaA_ijkl E(M)_i E(3)_j E(3)_k E(N)_l
C
C Output:
C     B11,B12,B22... Components of symmetric 2*2 matrix
C             2*{E-[E+B*VS**(-2)]**(-1/2)}*VS**2,
C             where E is the unit matrix and B is the input matrix.
C
C Date: 1999, March 13
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL VS2,A11,A12,A22,AD,AT,AS,BD,BT,BS
C
C     B=E+B*VS**(-2):
      VS2=VS**2
      B11=B11/VS2+1.
      B12=B12/VS2
      B22=B22/VS2+1.
C
C     A=SQRT(B):
      BD=B11*B22-B12*B12
      BT=B11+B22
      BS=B11-B22
      AD=SQRT(BD)
      AT=SQRT(BT+AD+AD)
      AS=BS/AT
      A11=0.5*(AT+AS)
      A12=B12/AT
      A22=0.5*(AT-AS)
C
C     B=2*[A**(-1)-E]*VS**2:
      B11=2.*(1.-A22/AD)*VS2
      B12=2.*(   A12/AD)*VS2
      B22=2.*(1.-A11/AD)*VS2
C
      RETURN
      END
