C 
C Subroutine file 'trans.for' to transform the computed quantities at an
C interface
C
C Date: 1997, June 27
C Coded by: Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     TRANS...Subroutine transforming the computed quantities across a
C             curved interface, see
C             C.R.T.5.9.
C             TRANS
C     CONV... Subroutine designed to replace the amplitudes Y(28) to
C             Y(IY(1)) expressed in the ray-centred coordinate system by
C             the amplitudes involving appropriate conversion
C             coefficients, see
C             C.R.T.5.5.4.
C             CONV
C
C=======================================================================
C
C     
C
      SUBROUTINE TRANS(YL,Y,YY,IY,KODNEW,ICBNEW,IEND)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12),KODNEW,ICBNEW,IEND
C
C This subroutine transforms the computed quantities across a curved
C interface, see
C C.R.T.5.9.
C
C Input:
C     YL...   Array containing local quantities at the point of
C             incidence.
C     Y...    Array containing basic quantities computed along the ray,
C             corresponding to the incident wave.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray, corresponding to the incident wave.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray, corresponding to the incident wave.
C     KODNEW..Position in the code (index in the array KODE)
C             corresponding to the next element of the ray.  Output from
C             subroutine CODE.
C     ICBNEW..The index of the complex block in which the next element
C             of the ray is to be situated, supplemented by the sign '+'
C             for P wave or '-' for S wave. Output from subroutine CODE.
C
C Output:
C     YL...   Array containing local quantities at the R/T point,
C             corresponding to the generated wave.
C     Y...    Array containing basic quantities computed along the ray,
C             corresponding to the generated wave.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray, corresponding to the generated wave.  Unchanged
C             input values.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray, corresponding to the generated wave.
C     IEND... Indication of the termination of the ray computation, see
C             C.R.T.5.4.
C             IEND.EQ.0... Computation of the ray continues.
C             IEND.NE.0... The computation of the ray terminates,
C               different values of IEND may specify the reason of the
C               termination:
C             24... Transmission is required by the code at a free
C               surface.
C             25... Ray of the required reflected or transmitted wave
C               is not real-valued (overcritical reflection or
C               transmission).
C             26... S wave in a liquid block is required by the code.
C             32... Reflection or type conversion at the fictitious
C               part of the interface.  Here, the interface is
C               considered to be fictitious if P and S wave velocities
C               and the density are the same at both sides of the
C               interface at the point of incidence.
C             33... Zero reflection/transmission coefficient.
C             37... The next element of the ray is required by the code
C               to be situated in a complex block that does not touch
C               the point of incidence.
C
C Subroutines and external functions required:
      EXTERNAL ISIDE,VELOC,KOOR,METRIC,SRFC2,PARM2,SMVPRD,COEF50,COEFSH
      INTEGER ISIDE,KOOR
C     ISIDE,VELOC... File 'model.for' of the package 'MODEL'.
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C              files (like 'val.for' and 'fit.for') of the package
C              'MODEL'.
C     PARM2 and subsequent routines... File 'parm.for' and subsequent
C              files (like 'val.for' and 'fit.for') of the package
C              'MODEL'.
C     SMVPRD... File 'means.for' of the package 'MODEL'.
C     COEF50,COEFSH... File 'coef.for'.
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
C     Auxiliary storage locations:
C
      INTEGER NAMPL1,NAMPL2,ISB1,ICB1,ICB2,I,J,NCODE,ND
      REAL VP1,VS1,RO1,VP2,VS2,RO2,V1,VD1(4),V2,VD2(10)
      REAL H11,H13,H43,HHH,Z11,Z12,Z13,Z41,Z42,Z43,ZZZ
      REAL H21,H23,H53,    Z21,Z22,Z23,Z51,Z52,Z53
      REAL H31,H33,H63,    Z31,Z32,Z33,Z61,Z62,Z63
      REAL E11,E21,E31,            C1,S1,C2,S2,C,S
      REAL V11,V21,V31,V12,V22,V32,D11,D12,D22,AUX1,AUX2,AUX3
      REAL P,AUX,E12,CFC11,CFC12,CFC22,Q1,Q2,P1,P2
      REAL COEF,RMOD,RPH,RMODSH,RPHSH,RR,YR
      REAL R1A1,Y1A1,R2A1,Y2A1,R1B1,Y1B1,R2B1,Y2B1
      REAL R1A2,Y1A2,R2A2,Y2A2,R1B2,Y1B2,R2B2,Y2B2
C
C     NAMPL1..Number of the amplitude coefficients of the incident wave
C             (5.9.1a).
C     NAMPL2..Number of the amplitude coefficients of the R/T wave
C             (5.9.1b).
C     ISB1... Index of a simple block of the incident wave (5.9.1e),
C             (5.9.3).
C     ICB1... Index of a complex block of the incident wave (5.9.1e).
C     ICB2... Index of a complex block at the other side (5.9.1e).
C     I,J...  Temporary storage locations (5.9.3), (5.9.5).
C     NCODE,ND... Type of the R/T coefficient (5.9.6).
C     VP1,VS1,RO1... Parameters of the complex block IABS(ICB1) (5.9.2).
C     VP2,VS2,RO2... Parameters of the complex block ICB2 (5.9.2).
C     V1,VD1(2:4)... Incident wave velocity and its derivatives (5.9.2).
C     V2,VD2(2:4)...R/T wave velocity and its derivatives (5.9.2).
C     VD2(1),VD2(5:10)... Temporary storage locations (5.9.2).
C     H11,H21,H31,H13,H23,H33... Covariant components of the basis
C             vectors of ray-centred coordinate system of the incident
C             wave (5.9.3).
C     H43,H53,H63... Contravariant components of the slowness vector of
C             the incident wave (5.9.3).  Not used in Cartesian
C             coordinates.
C     HHH...  Norm of the slowness vector (5.9.3).
C     Z11,Z21,Z31,Z12,Z22,Z32,Z13,Z23,Z33... Covariant components of the
C             local Cartesian coordinate system (5.9.3).
C     Z41,Z51,Z61,Z42,Z52,Z62,Z43,Z53,Z63... Contravariant components of
C             the local Cartesian coordinate system (5.9.3).
C     ZZZ...  Norm of the gradient of the function describing the
C             interface (5.9.3,5.9.4).
C     E11,E21,E31... Covariant components of the unit vector
C             perpendicular to the R/T ray in the plane of incidence
C             (5.9.3).
C             E11 is also a temporary storage location in (5.9.5).
C     E13,E23,E33... Covariant components of the unit vector tangent to
C             the R/T ray (5.9.3), not used in this version.
C     C1,S1...Cosine and sine of the angle of incidence (5.9.3,5.9.5).
C     C2,S2...Cosine and sine of the R/T angle (5.9.3,5.9.5).
C     C,S...  Cosine and sine of the revolution angle
C             (5.9.3,5.9.5,5.9.6).
C     V11,V21,V31... Velocity gradient in local Cartesian coordinates
C             corresponding to the incident wave (5.9.4).
C     V12,V22,V32... Velocity gradient in local Cartesian coordinates
C             corresponding to the incident wave (5.9.4).
C     D11,D12,D22... Matrix of the curvature of the interface (5.9.4).
C     AUX1,AUX2,AUX3... Temporary storage locations (5.9.4).
C     P...    Snell ray parameter(5.9.5,5.9.6)
C     AUX,E12,CFC11,CFC12,CFC22,Q1,Q2,P1,P2... Temporary storage
C             locations (5.9.5).
C     COEF... (5.9.6).
C     RMOD,RPH... P-SV R/T coefficient (5.9.6).
C     RMODSH,RPHSH... SH R/T coefficient (5.9.6).
C     RR,YR...Real and imaginary part of a complex-valued R/T
C             coefficient (5.9.6).
C     ( R1A1,Y1A1,R2A1,Y2A1 )
C     ( R1B1,Y1B1,R2B1,Y2B1 )... Amplitude coefficients (5.9.6):
C     ( R1A2,Y1A2,R2A2,Y2A2 )
C     ( R1B2,Y1B2,R2B2,Y2B2 )               Initial
C             Part:          Component:     polarization:  Wave:
C             R...Real       1...P-SV       A...First      1...Incident
C             Y...Imaginary  2...SH         B...Second     2...R/T
C
C.......................................................................
C
C
C     (5.9.1) Transformation of auxiliary quantities, travel time and
C     coordinates:
C     (a) Number of real-valued quantities describing the reduced
C     amplitudes of the incident wave:
      NAMPL1=IY(1)-27
C     (b) Number of real-valued quantities describing the reduced
C     amplitudes of the generated wave:
      IF(IY(5).GT.0.AND.ICBNEW.LT.0) THEN
        NAMPL2=NAMPL1*2
      ELSE IF(IY(5).LT.0.AND.ICBNEW.GT.0) THEN
        NAMPL2=NAMPL1/2
      ELSE
        NAMPL2=NAMPL1
      END IF
C     (c) Output number of basic quantities:
      IY(1)=27+NAMPL2
C     (d) New position in the code:
      IY(2)=KODNEW
C     (e) Indices of simple and complex blocks:
      ISB1=IY(4)
      ICB1=IY(5)
      ICB2=IABS(IY(8))
      IF(IABS(ICBNEW).EQ.ICB2) THEN
C       Transmission:
        IY(3)=IABS(IY(5))
        IY(4)=IY(7)
        IY(5)=ICBNEW
        IY(6)=-IY(6)
      ELSE IF(IABS(ICBNEW).EQ.IABS(IY(5))) THEN
C       Reflection:
        IY(5)=ICBNEW
      ELSE
C       Required complex block is not attainable:
        IEND=37
        RETURN
      END IF
C     (f) IY(7), IY(8) may be undefined
      IY(7)=0
      IY(8)=0
C     (g) Number of transformations at an interface:
      IY(11)=IY(11)+1
C     (h) IY(9), IY(10), IY(12), YY(1:5) are unchanged
C     (i) Travel time and coordinates Y(1:5) remain unchanged
C
C     (5.9.2) Velocities (metric tensor is determined later on)
C     (b) Material parameters corresponding to the incident wave:
      VP1=YL(1)
      VS1=YL(2)
      RO1=YL(3)
      VD1(2)=YL(4)
      VD1(3)=YL(5)
      VD1(4)=YL(6)
      IF(ICB1.GE.0) THEN
        V1=VP1
      ELSE
        V1=VS1
      END IF
C     (c) Material parameters on the other side of the interface:
      IF(ICB2.NE.0) THEN
        CALL PARM2(ICB2,Y(3),UP,US,RO2,QP,QS)
        CALL VELOC(ISIGN(ICB2,ICBNEW),UP,US,QP,QS,VP2,VS2,VD2,QL)
      ELSE
C       Free space:
        RO2=0.
      END IF
      IF(VP1.EQ.VP2) THEN
        IF(ICBNEW.NE.ISIGN(ICB2,ICB1)) THEN
C         Transmission without conversion is excluded
          IF(VS1.EQ.VS2.AND.RO1.EQ.RO2) THEN
C           Reflection or wave type conversion at a fictitious interface
            IEND=32
            RETURN
          END IF
        END IF
      END IF
C     (d) Material parameters corresponding to the generated wave:
      IF(ICBNEW.EQ.ICB1) THEN
C       Reflection without conversion:
        V2=V1
        VD2(2)=VD1(2)
        VD2(3)=VD1(3)
        VD2(4)=VD1(4)
      ELSE IF(ICBNEW.EQ.-ICB1) THEN
C       Reflection with conversion:
        CALL PARM2(IABS(ICBNEW),Y(3),UP,US,AUX,QP,QS)
        CALL VELOC(ICBNEW,UP,US,QP,QS,AUX,AUX,VD2,QL)
        V2=VD2(1)
      ELSE
C       Transmission:
        V2=VD2(1)
        IF(RO2.EQ.0.) THEN
C         Transmission into free space:
          IEND=24
          RETURN
        END IF
        IF(V2.EQ.0.) THEN
C         Zero velocity (S-wave in liquid):
          IEND=26
          RETURN
        END IF
        YL(1)=VP2
        YL(2)=VS2
        YL(3)=RO2
      END IF
      YL(4)=VD2(2)
      YL(5)=VD2(3)
      YL(6)=VD2(4)
C
C     (5.9.3) Transformation of the slowness vector and of the
C     basis of the ray-centred coordinate system:
      CALL SRFC2(IABS(IY(6)),Y(3),FAUX)
C     Non-unit vectors - covariant components
      H13=Y(6)
      H23=Y(7)
      H33=Y(8)
      H11=Y(9)
      H21=Y(10)
      H31=Y(11)
      Z13=FAUX(2)
      Z23=FAUX(3)
      Z33=FAUX(4)
C     (c) Non-unit vectors - contravariant components
      IF(KOOR().EQ.0) THEN
C       Cartesian coordinates:
C       Norm of the slowness vector
        HHH=SQRT(H13*H13+H23*H23+H33*H33)
C       Norm of the gradient of the function describing the interface
        ZZZ=SQRT(Z13*Z13+Z23*Z23+Z33*Z33)
C       Unit normal to the interface:
C       Covariant components
        Z13=Z13/ZZZ
        Z23=Z23/ZZZ
        Z33=Z33/ZZZ
C       Contravariant components
        Z43=Z13
        Z53=Z23
        Z63=Z33
      ELSE
C       Curvilinear coordinates:
        CALL METRIC(Y(3),GSQRD,G,GAMMA)
C       Second covariant derivatives of the interface function
        FAUX(5)=FAUX(5)
     *             -GAMMA(1)*FAUX(2)-GAMMA(7)*FAUX(3)-GAMMA(13)*FAUX(4)
        FAUX(6)=FAUX(6)
     *             -GAMMA(2)*FAUX(2)-GAMMA(8)*FAUX(3)-GAMMA(14)*FAUX(4)
        FAUX(7)=FAUX(7)
     *             -GAMMA(3)*FAUX(2)-GAMMA(9)*FAUX(3)-GAMMA(15)*FAUX(4)
        FAUX(8)=FAUX(8)
     *             -GAMMA(4)*FAUX(2)-GAMMA(10)*FAUX(3)-GAMMA(16)*FAUX(4)
        FAUX(9)=FAUX(9)
     *             -GAMMA(5)*FAUX(2)-GAMMA(11)*FAUX(3)-GAMMA(17)*FAUX(4)
        FAUX(10)=FAUX(10)
     *             -GAMMA(6)*FAUX(2)-GAMMA(12)*FAUX(3)-GAMMA(18)*FAUX(4)
C       Slowness vector - contravariant components
        CALL SMVPRD(G(7),H13,H23,H33,H43,H53,H63)
C       Norm of the slowness vector
        HHH=SQRT(H13*H43+H23*H53+H33*H63)
C       Contravariant gradient of the function describing the interface
        CALL SMVPRD(G(7),Z13,Z23,Z33,Z43,Z53,Z63)
C       Norm of the gradient of the function describing the interface
        ZZZ=SQRT(Z13*Z43+Z23*Z53+Z33*Z63)
C       Unit normal to the interface:
C       Covariant components
        Z13=Z13/ZZZ
        Z23=Z23/ZZZ
        Z33=Z33/ZZZ
C       Contravariant components
        Z43=Z43/ZZZ
        Z53=Z53/ZZZ
        Z63=Z63/ZZZ
      END IF
C     Unit vector tangent to the incident ray
      H13=H13/HHH
      H23=H23/HHH
      H33=H33/HHH
C     Cosine of the angle of incidence
      C1=H13*Z43+H23*Z53+H33*Z63
C     Vector tangent to the interface in the plane of incidence
      Z11=H13-Z13*C1
      Z21=H23-Z23*C1
      Z31=H33-Z33*C1
      IF(KOOR().EQ.0) THEN
C       Cartesian coordinates:
        AUX=Z11*Z11+Z21*Z21+Z31*Z31
      ELSE
C       Curvilinear coordinates:
        Z41=H43/HHH-Z43*C1
        Z51=H53/HHH-Z53*C1
        Z61=H63/HHH-Z63*C1
        AUX=Z11*Z41+Z21*Z51+Z31*Z61
      END IF
C     Sine of the angle of incidence
      IF(AUX.LT.0.5) THEN
        S1=SQRT(AUX)
      ELSE
        S1=SQRT(1.-C1*C1)
      END IF
      IF(S1.LT.0.0001) THEN
C       Correction for nearly normal incidence
        AUX1=Z43*H11+Z53*H21+Z63*H31
        AUX2=SQRT(1.-AUX1*AUX1)
        Z11=(H11-Z13*AUX1)/AUX2
        Z21=(H21-Z23*AUX1)/AUX2
        Z31=(H31-Z33*AUX1)/AUX2
      ELSE
        Z11=Z11/S1
        Z21=Z21/S1
        Z31=Z31/S1
      END IF
C     Angle of reflection/transmission
      S2=S1*V2/V1
      C2=1.-S2*S2
      IF(C2.LE.0.) THEN
C       Overcritical ray
        IEND=25
        RETURN
      ELSE
        C2=SIGN(SQRT(C2),C1)
      END IF
      IF(IABS(ICBNEW).EQ.IABS(ICB1)) C2=-C2
C     Correction for nearly grazing rays (not included in the paper)
      I=ISIDE(IY(6),ISB1)
      IF(I.EQ.2) THEN
C       591
        CALL ERROR('591 in TRANS: Wrong interface')
C       This error should not appear.  Contact the authors.
      ELSE IF(FLOAT(I)*C1.GT.0.) THEN
C       Incident ray goes from the interface - interchanging R/T
        C2=-C2
      END IF
      IF(KOOR().EQ.0) THEN
C       Cartesian coordinates:
C       Vectors tangent to the interface:
C       Vector in the plane of incidence - contravariant components
        Z41=Z11
        Z51=Z21
        Z61=Z31
C       Vector perpendicular to the plane of incidence - covariant comp.
        Z12=Z23*Z31-Z33*Z21
        Z22=Z33*Z11-Z13*Z31
        Z32=Z13*Z21-Z23*Z11
C       Vector perpendicular to the plane of incidence - contravariant
        Z42=Z12
        Z52=Z22
        Z62=Z32
      ELSE
C       Curvilinear coordinates:
C       Vectors tangent to the interface:
C       Vector in the plane of incidence - contravariant components
        CALL SMVPRD(G(7),Z11,Z21,Z31,Z41,Z51,Z61)
C       Vector perpendicular to the plane of incidence - covariant comp.
        Z12=(Z53*Z61-Z63*Z51)*GSQRD
        Z22=(Z63*Z41-Z43*Z61)*GSQRD
        Z32=(Z43*Z51-Z53*Z41)*GSQRD
C       Vector perpendicular to the plane of incidence - contravariant
        Z42=(Z23*Z31-Z33*Z21)/GSQRD
        Z52=(Z33*Z11-Z13*Z31)/GSQRD
        Z62=(Z13*Z21-Z23*Z11)/GSQRD
      END IF
C     Revolution angle of the ray centred coordinate system
      C=(H11*Z41+H21*Z51+H31*Z61)/C1
      S=-H11*Z42-H21*Z52-H31*Z62
C     Unit vector perpendicular to the R/T ray in the plane of incidence
      E11=Z11*C2-Z13*S2
      E21=Z21*C2-Z23*S2
      E31=Z31*C2-Z33*S2
C     Unit vector tangent to the R/T ray
C1    E13=Z11*S2+Z13*C2
C2    E23=Z21*S2+Z23*C2
C3    E33=Z31*S2+Z33*C2
C     Slowness vector
C4    Y(6)=E13/V2
C5    Y(7)=E23/V2
C6    Y(8)=E33/V2
C     For a nearly normal incidence, the vector (Z11,Z21,Z31) is not
C     situated in the plane of incidence and the above six equations
C     cannot be used.  Thus, they are replaced by the following three
C     ones
      Y(6)=Y(6)+Z13*(C2/V2-C1/V1)
      Y(7)=Y(7)+Z23*(C2/V2-C1/V1)
      Y(8)=Y(8)+Z33*(C2/V2-C1/V1)
C     First basis vector of the ray centred coordinate system
      Y(9) =E11*C-Z12*S
      Y(10)=E21*C-Z22*S
      Y(11)=E31*C-Z32*S
C
C     (5.9.4) Curvature of the interface and velocity gradients in the
C     local Cartesian coordinate system:
      CALL SMVPRD(FAUX(5),Z41,Z51,Z61,AUX1,AUX2,AUX3)
C     Note: In the original paper on C.R.T., the curvature matrix D is
C     erroneously defined with opposite sign.
      D11=-(AUX1*Z41+AUX2*Z51+AUX3*Z61)/ZZZ
      D12=-(AUX1*Z42+AUX2*Z52+AUX3*Z62)/ZZZ
      CALL SMVPRD(FAUX(5),Z42,Z52,Z62,AUX1,AUX2,AUX3)
      D22=-(AUX1*Z42+AUX2*Z52+AUX3*Z62)/ZZZ
      V11=VD1(2)*Z41+VD1(3)*Z51+VD1(4)*Z61
      V21=VD1(2)*Z42+VD1(3)*Z52+VD1(4)*Z62
      V31=VD1(2)*Z43+VD1(3)*Z53+VD1(4)*Z63
      V12=VD2(2)*Z41+VD2(3)*Z51+VD2(4)*Z61
      V22=VD2(2)*Z42+VD2(3)*Z52+VD2(4)*Z62
      V32=VD2(2)*Z43+VD2(3)*Z53+VD2(4)*Z63
C
C     (5.9.5) Dynamic ray tracing across a curved interface
      P=S1/V1
      AUX=C1/V1-C2/V2
      E11=(S1*C1*V31-(1.+C1*C1)*V11)/V1-(S2*C2*V32-(1.+C2*C2)*V12)/V2
      E12=V22/V2-V21/V1
      CFC11=(AUX*D11+P*E11)/C2/C2
      CFC12=(AUX*D12+P*E12)/C2
      CFC22= AUX*D22
      AUX=C1/C2
      DO 51 I=1,4
        J=4*I
        Q1= C*Y(8+J)+S*Y(9+J)
        Q2=-S*Y(8+J)+C*Y(9+J)
        P1= C*Y(10+J)+S*Y(11+J)
        P2=-S*Y(10+J)+C*Y(11+J)
        Q1=Q1/AUX
        P1=P1*AUX+CFC11*Q1+CFC12*Q2
        P2=P2    +CFC12*Q1+CFC22*Q2
        Y(8+J)=C*Q1-S*Q2
        Y(9+J)=S*Q1+C*Q2
        Y(10+J)=C*P1-S*P2
        Y(11+J)=S*P1+C*P2
   51 CONTINUE
C
C     (5.9.6) Transformation of reduced amplitudes
      IF(IABS(ICBNEW).EQ.IABS(ICB1)) THEN
C       Reflection:
        NCODE=1
        COEF=1.
      ELSE
C       Transmission:
        NCODE=3
        COEF=RO2/RO1
      END IF
      COEF=SQRT(COEF*ABS(C2/C1)*V2/V1)
      IF(C1.GE.0) THEN
C       Epsilon=+1.
C       (see C.R.T.5.9.7)
        ND=0
      ELSE
C       Epsilon=-1.
C       (see C.R.T.5.9.7)
        ND=1
      END IF
      RMODSH=0.
      R2A2=0.
      Y2A2=0.
      R2B2=0.
      Y2B2=0.
      IF(ICB1.GE.0) THEN
C       Incident P wave:
        IF(ICBNEW.LT.0) THEN
C         Outgoing S wave:
          NCODE=NCODE+1
        END IF
        R1A1=Y(28)
        Y1A1=Y(29)
        IF(NAMPL1.GT.2) THEN
          R1B1= Y(30)
          Y1B1= Y(31)
        END IF
      ELSE
C       Incident S wave:
        IF(ICB2.EQ.0) THEN
          NCODE=NCODE+2
        ELSE
          NCODE=NCODE+4
        END IF
        R1A1= C*Y(28)+S*Y(30)
        Y1A1= C*Y(29)+S*Y(31)
        IF(NAMPL1.GT.4) THEN
          R1B1= C*Y(32)+S*Y(34)
          Y1B1= C*Y(33)+S*Y(35)
        END IF
        IF(ICBNEW.LT.0) THEN
C         Outgoing S wave:
          NCODE=NCODE+1
          IF(ICBNEW.EQ.ICB1) THEN
            CALL COEFSH(P,VS1,RO1,VS2,RO2,1,RMODSH,RPHSH)
          ELSE
            CALL COEFSH(P,VS1,RO1,VS2,RO2,2,RMODSH,RPHSH)
          END IF
          RMODSH=COEF*RMODSH
          RR=RMODSH*COS(RPHSH)
          YR=RMODSH*SIN(RPHSH)
          R2A1=-S*Y(28)+C*Y(30)
          Y2A1=-S*Y(29)+C*Y(31)
          R2A2=RR*R2A1-YR*Y2A1
          Y2A2=YR*R2A1+RR*Y2A1
          IF(NAMPL1.GT.4) THEN
            R2B1=-S*Y(32)+C*Y(34)
            Y2B1=-S*Y(33)+C*Y(35)
            R2B2=RR*R2B1-YR*Y2B1
            Y2B2=YR*R2B1+RR*Y2B1
          END IF
        END IF
      END IF
      CALL COEF50(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH)
      RMOD=COEF*RMOD
      IF(RMOD.EQ.0..AND.RMODSH.EQ.0.) THEN
C       Zero reflection/transmission coefficient
        IEND=33
        RETURN
      END IF
      RR=RMOD*COS(RPH)
      YR=RMOD*SIN(RPH)
      R1A2=RR*R1A1-YR*Y1A1
      Y1A2=YR*R1A1+RR*Y1A1
      IF(ICBNEW.GE.0) THEN
C       Outgoing P wave:
        Y(28)=R1A2
        Y(29)=Y1A2
        IF(NAMPL2.GT.2) THEN
          R1B2=RR*R1B1-YR*Y1B1
          Y1B2=YR*R1B1+RR*Y1B1
          Y(30)=R1B2
          Y(31)=Y1B2
        END IF
      ELSE
C       Outgoing S wave:
        Y(28)= C*R1A2-S*R2A2
        Y(29)= C*Y1A2-S*Y2A2
        Y(30)= S*R1A2+C*R2A2
        Y(31)= S*Y1A2+C*Y2A2
        IF(NAMPL2.GT.4) THEN
          R1B2=RR*R1B1-YR*Y1B1
          Y1B2=YR*R1B1+RR*Y1B1
          Y(32)= C*R1B2-S*R2B2
          Y(33)= C*Y1B2-S*Y2B2
          Y(34)= S*R1B2+C*R2B2
          Y(35)= S*Y1B2+C*Y2B2
        END IF
      END IF
C
      IEND=0
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE CONV(KSTORE,YL,Y,IY,NAMPL,YC)
      INTEGER KSTORE,IY(12),NAMPL
      REAL YL(6),Y(35),YC(12)
C
C This subroutine replaces the amplitudes Y(28) to Y(IY(1)) expressed in
C the ray-centred coordinate system by the amplitudes involving
C appropriate conversion coefficients, see
C C.R.T.5.5.4.
C
C Input:
C     KSTORE..Specifies whether the conversion coefficients are to be
C             considered, see
C             C.R.T.5.5.4
C             and 5.6e.
C             KSTORE.LE.1... No amplitude conversion coefficients are
C               applied at the point of intersection with an interface.
C             KSTORE.GE.2... Amplitude conversion coefficients are
C               applied at the point of intersection with an interface.
C     YL...   Array containing local quantities at the point of
C             incidence.
C     Y...    Array of the dimension at least 39.  The locations Y(1) to
C             Y(IY(1)) contain basic quantities computed along the ray,
C             corresponding to the incident wave.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray, corresponding to the incident wave.
C     YC...   Array of the dimension at least 12.
C
C Output:
C     NAMPL...Number of real quantities representing complex-valued
C             vectorial reduced amplitudes.  If no conversion
C             coefficients are applied NAMPL=IY(1)-27,  otherwise
C             NAMPL=6 or 12, see
C             C.R.T.5.5.4.
C     YC...   Array containing real quantities representing complex-
C             -valued vectorial reduced amplitudes.  If no conversion
C             coefficients are applied, YC is a copy of Y(28) to
C             Y(IY(1)).  Otherwise,  YC represents the vectorial reduced
C             amplitudes involving appropriate conversion coefficients,
C             expressed in ray-centred coordinate system, see
C             C.R.T.5.5.4.
C             P wave at the initial point of the ray:
C               NAMPL=6,
C               YC(1)=REAL(A13),  YC(2)=AIMAG(A13),
C               YC(3)=REAL(A23),  YC(4)=AIMAG(A23),
C               YC(5)=REAL(A33),  YC(6)=AIMAG(A33).
C             S wave at the initial point of the ray:
C               NAMPL=12,
C               YC(1)=REAL(A11),  YC(2)=AIMAG(A11),
C               YC(3)=REAL(A21),  YC(4)=AIMAG(A21),
C               YC(5)=REAL(A31),  YC(6)=AIMAG(A31),
C               YC(7)=REAL(A12),  YC(8)=AIMAG(A12),
C               YC(9)=REAL(A22),  YC(10)=AIMAG(A22),
C               YC(11)=REAL(A32), YC(12)=AIMAG(A32).
C
C Subroutines and external functions required:
      EXTERNAL NSRFC,VELOC,KOOR,METRIC,SRFC2,PARM2,SMVPRD,COEF50
      INTEGER NSRFC,KOOR
C     NSRFC,VELOC... File 'model.for' of the package 'MODEL'.
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C              files (like 'val.for' and 'fit.for') of the package
C              'MODEL'.
C     PARM2 and subsequent routines... File 'parm.for' and subsequent
C              files (like 'val.for' and 'fit.for') of the package
C              'MODEL'.
C     SMVPRD... File 'means.for' of the package 'MODEL'.
C     COEF50... File 'coef.for'.
C
C Date: 1997, April 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
C     Auxiliary storage locations:
C
      INTEGER ICB2,NCODE,ND,NW,IW,I
      REAL VP1,VS1,RO1,V1,VP2,VS2,RO2,V2
      REAL H11,H13,H43,HHH,Z13,Z42,Z43,ZZZ
      REAL H21,H23,H53,    Z23,Z52,Z53
      REAL H31,H33,H63,    Z33,Z62,Z63
      REAL P,C1,S1,C2,S2,C,S,CR,CI,SR,SI
      REAL RMOD,RPH,RR,YR
      REAL R1A1,Y1A1,R2A1,Y2A1,R1B1,Y1B1,R2B1,Y2B1
      REAL R1A2,Y1A2,R2A2,Y2A2,R1B2,Y1B2,R2B2,Y2B2,RA2,YA2,RB2,YB2
C
C     ICB2... Index of a complex block at the other side of the
C             interface.  Auxiliary variable.
C     NCODE,ND... Type of the R/T coefficient.
C     NW...   Total number of R/T waves.
C     IW...   Loop variable. Index of an R/T wave.
C     I...    Auxiliary loop variable.
C     VP1,VS1,RO1... Parameters of the complex block IABS(IY(5)).
C     V1...   Incident wave velocity.
C     VP2,VS2,RO2... Parameters of the complex block ICB2.
C     V2...   R/T wave velocity.
C     H11,H21,H31,H13,H23,H33... Covariant components of the basis
C             vectors of ray-centred coordinate system of the incident
C             wave.
C     H43,H53,H63... Contravariant components of the slowness vector of
C             the incident wave.  Not used in Cartesian coordinates.
C     HHH...  Norm of the slowness vector.
C     Z13,Z23,Z33... Covariant components of the local Cartesian
C             coordinate system.
C     Z42,Z52,Z62,Z43,Z53,Z63... Contravariant components of the local
C             Cartesian coordinate system.
C     ZZZ...  Norm of the gradient of the function describing the
C             interface.
C     C1,S1...Cosine and sine of the angle of incidence.
C     C2,S2...Cosine and sine of the R/T angle.  Auxiliary variables.
C     C,S...  Cosine and sine of the revolution angle.
C     P...    Snell ray parameter.
C     CR,CI,SR,SI...Real and imaginary parts of the cosine and sine of
C             the rotation angle in the plane of incidence.
C     RMOD,RPH... Any R/T coefficient.
C     RR,YR...Real and imaginary part of a complex-valued R/T
C             coefficient.
C     ( R1A1,Y1A1,R2A1,Y2A1 )... Amplitude coefficients in ray-centred
C     ( R1B1,Y1B1,R2B1,Y2B1 )    coordinate system of the incident wave
C     ( R1A2,Y1A2,R2A2,Y2A2 )
C     ( R1B2,Y1B2,R2B2,Y2B2 )               Initial
C             Part:          Component:     polarization:  Wave:
C             R...Real       1...P-SV, SV   A...First      1...Incident
C             Y...Imaginary  2...SH I  I    B...Second     2...R/T
C                                   I  I...R/T wave
C                                   I...Incident wave
C     ( RA2,YA2 )... P-SV amplitude coefficients in ray-centred
C     ( RB2,YB2 )    coordinate system of the R/T wave
C             Part:          Polarization:  Wave:
C             R...Real       A...First      2...R/T
C             Y...Imaginary  B...Second
C
C.......................................................................
C
C     Check whether the amplitude conversion coefficients are applied
      IF(IABS(IY(6)).GT.NSRFC().OR.KSTORE.LE.1) THEN
        NAMPL=IY(1)-27
        DO 10 I=1,NAMPL
          YC(I)=Y(27+I)
   10   CONTINUE
        RETURN
      END IF
C
C     Material parameters corresponding to the incident wave:
      VP1=YL(1)
      VS1=YL(2)
      RO1=YL(3)
      IF(IY(5).GE.0) THEN
        V1=VP1
      ELSE
        V1=VS1
      END IF
C     Material parameters on the other side of the interface:
      ICB2=IABS(IY(8))
      IF(ICB2.NE.0) THEN
        CALL PARM2(ICB2,Y(3),UP,US,RO2,QP,QS)
        CALL VELOC(ICB2,UP,US,QP,QS,VP2,VS2,VD,QL)
        NW=4
      ELSE
C       Free space:
        RO2=0.
        NW=2
      END IF
C     VP1, VS1, RO1, V1, VP2, VS2, RO2 and NW are defined.
C
C     Basis of the ray-centred coordinate system and the local Cartesian
C     coordinate system connected with the interface:
      CALL SRFC2(IABS(IY(6)),Y(3),FAUX)
C     Non-unit vectors - covariant components
      H13=Y(6)
      H23=Y(7)
      H33=Y(8)
      H11=Y(9)
      H21=Y(10)
      H31=Y(11)
      Z13=FAUX(2)
      Z23=FAUX(3)
      Z33=FAUX(4)
C     (c) Non-unit vectors - contravariant components
      IF(KOOR().EQ.0) THEN
C       Cartesian coordinates:
        GSQRD=1.
C       Norm of the slowness vector
        HHH=SQRT(H13*H13+H23*H23+H33*H33)
C       Norm of the gradient of the function describing the interface
        ZZZ=SQRT(Z13*Z13+Z23*Z23+Z33*Z33)
C       Unit normal to the interface:
C       Covariant components
        Z13=Z13/ZZZ
        Z23=Z23/ZZZ
        Z33=Z33/ZZZ
C       Contravariant components
        Z43=Z13
        Z53=Z23
        Z63=Z33
      ELSE
C       Curvilinear coordinates:
        CALL METRIC(Y(3),GSQRD,G,GAMMA)
C       Slowness vector - contravariant components
        CALL SMVPRD(G(7),H13,H23,H33,H43,H53,H63)
C       Norm of the slowness vector
        HHH=SQRT(H13*H43+H23*H53+H33*H63)
C       Contravariant gradient of the function describing the interface
        CALL SMVPRD(G(7),Z13,Z23,Z33,Z43,Z53,Z63)
C       Norm of the gradient of the function describing the interface
        ZZZ=SQRT(Z13*Z43+Z23*Z53+Z33*Z63)
C       Unit normal to the interface:
C       Covariant components
        Z13=Z13/ZZZ
        Z23=Z23/ZZZ
        Z33=Z33/ZZZ
C       Contravariant components
        Z43=Z43/ZZZ
        Z53=Z53/ZZZ
        Z63=Z63/ZZZ
      END IF
C     Unit vector tangent to the incident ray
      H13=H13/HHH
      H23=H23/HHH
      H33=H33/HHH
C     Cosine of the angle of incidence
      C1=H13*Z43+H23*Z53+H33*Z63
C     Sine of the angle of incidence
      S1=SQRT(1.-C1*C1)
C     C1 and S1 are the cosine and sine of the angle of incidence.
C
C     Revolution angle of the basis of the ray-centred coordinate system
C     around ray:
      IF(S1.LT.0.0001) THEN
C       Nearly normal incidence
        C=1.
        S=0.
      ELSE
C       Vector perpendicular to the plane of incidence - contravariant
        Z42=(Z23*H33-Z33*H23)
        Z52=(Z33*H13-Z13*H33)
        Z62=(Z13*H23-Z23*H13)
C       Revolution angle of the ray centred coordinate system
        C=-(H11*Z43+H23*Z53+H31*Z63)/S1
        S=-(H11*Z42+H21*Z52+H31*Z62)/S1/GSQRD
      END IF
C     C and S are the cosine and sine of the revolution angle.
C
C     Amplitudes of the incident wave in the rotated ray-centred
C     coordinate system of the incident wave:
      NAMPL=6
      IF(IY(5).GE.0) THEN
C       Incident P wave:
        NCODE=0
        R1A1=Y(28)
        Y1A1=Y(29)
        R1A2=0.
        Y1A2=0.
        R2A2=0.
        Y2A2=0.
        YC(5)=R1A1
        YC(6)=Y1A1
        IF(IY(1).GT.29) THEN
          NAMPL=12
          R1B1= Y(30)
          Y1B1= Y(31)
          R1B2=0.
          Y1B2=0.
          R2B2=0.
          Y2B2=0.
          YC(11)=R1B1
          YC(12)=Y1B1
        END IF
      ELSE
C       Incident S wave:
        NCODE=NW
        R1A1= C*Y(28)+S*Y(30)
        Y1A1= C*Y(29)+S*Y(31)
        R2A1=-S*Y(28)+C*Y(30)
        Y2A1=-S*Y(29)+C*Y(31)
        R1A2=R1A1
        Y1A2=Y1A1
        R2A2=R2A1
        Y2A2=Y2A1
        YC(5)=0.
        YC(6)=0.
        IF(IY(1).GT.31) THEN
          NAMPL=12
          R1B1= C*Y(32)+S*Y(34)
          Y1B1= C*Y(33)+S*Y(35)
          R2B1=-S*Y(32)+C*Y(34)
          Y2B1=-S*Y(33)+C*Y(35)
          R1B2=R1B1
          Y1B2=Y1B1
          R2B2=R2B1
          Y2B2=Y2B1
          YC(11)=0.
          YC(12)=0.
        END IF
      END IF
C
C     Snell parameter
      P=S1/V1
C
      IF(C1.GE.0) THEN
C       Epsilon=+1.
C       (see C.R.T.5.9.7)
        ND=0
      ELSE
C       Epsilon=-1.
C       (see C.R.T.5.9.7)
        ND=1
      END IF
C
C     Loop over the generated waves:
      DO 90 IW=1,2
        NCODE=NCODE+1
C
C       Propagation velocity of the generated wave
        IF(IW.EQ.1) THEN
          V2=VP1
        ELSE IF(IW.EQ.2) THEN
          V2=VS1
        ELSE IF(IW.EQ.3) THEN
          V2=VP2
        ELSE
          V2=VS2
        END IF
C
C       Rotation angle of the basis of the ray-centred coordinate system
C       in the plane of incidence:
        S2=P*V2
        C2=1.-S2*S2
        CR=S1*S2
        SR=C1*S2
        IF(C2.GE.0.) THEN
          C2=SIGN(SQRT(C2),C1)
          IF(IW.LE.2) C2=-C2
          CR=CR+C1*C2
          SR=SR-S1*C2
          CI=0.
          SI=0.
        ELSE
          C2=SIGN(SQRT(-C2),C1)
          IF(IW.LE.2) C2=-C2
          CI= C1*C2
          SI=-S1*C2
        END IF
C       Cosine of the rotation angle is CR+i*CI, where i=SQRT(-1),
C       sine of the rotation angle is   SR+i*SI.
C
C       Transformation of reduced amplitudes (5.9.6)
        CALL COEF50(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH)
        RR=RMOD*COS(RPH)
        YR=RMOD*SIN(RPH)
        RA2=RR*R1A1-YR*Y1A1
        YA2=YR*R1A1+RR*Y1A1
        IF(NAMPL.GT.6) THEN
          RB2=RR*R1B1-YR*Y1B1
          YB2=YR*R1B1+RR*Y1B1
        END IF
        IF(IW.EQ.1.OR.IW.EQ.3) THEN
C         Outgoing P wave:
          R1A2=R1A2+SR*RA2-SI*YA2
          Y1A2=Y1A2+SI*RA2+SR*YA2
          YC(5)=YC(5)+CR*RA2-CI*YA2
          YC(6)=YC(6)+CI*RA2+CR*YA2
          IF(NAMPL.GT.6) THEN
            R1B2=R1B2+SR*RB2-SI*YB2
            Y1B2=Y1B2+SI*RB2+SR*YB2
            YC(11)=YC(11)+CR*RB2-CI*YB2
            YC(12)=YC(12)+CI*RB2+CR*YB2
          END IF
        ELSE
C         Outgoing S wave:
          R1A2=R1A2+CR*RA2-CI*YA2
          Y1A2=Y1A2+CI*RA2+CR*YA2
          YC(5)=YC(5)-SR*RA2+SI*YA2
          YC(6)=YC(6)-SI*RA2-SR*YA2
          IF(NAMPL.GT.6) THEN
            R1B2=R1B2+CR*RB2-CI*YB2
            Y1B2=Y1B2+CI*RB2+CR*YB2
            YC(11)=YC(11)-SR*RB2+SI*YB2
            YC(12)=YC(12)-SI*RB2-SR*YB2
          END IF
          IF(IY(5).LT.0) THEN
C           SH wave coefficients:
            IF(IW.EQ.2) THEN
              CALL COEFSH(P,VS1,RO1,VS2,RO2,1,RMOD,RPH)
            ELSE
              CALL COEFSH(P,VS1,RO1,VS2,RO2,2,RMOD,RPH)
            END IF
            RR=RMOD*COS(RPH)
            YR=RMOD*SIN(RPH)
            R2A2=R2A2+RR*R2A1-YR*Y2A1
            Y2A2=Y2A2+YR*R2A1+RR*Y2A1
            IF(NAMPL.GT.6) THEN
              R2B2=R2B2+RR*R2B1-YR*Y2B1
              Y2B2=Y2B2+YR*R2B1+RR*Y2B1
            END IF
          END IF
        END IF
   90 CONTINUE
C
C     Rotation of the basis of the ray-centred coordinate system around
C     ray:
      YC(1)= C*R1A2-S*R2A2
      YC(2)= C*Y1A2-S*Y2A2
      YC(3)= S*R1A2+C*R2A2
      YC(4)= S*Y1A2+C*Y2A2
      IF(NAMPL.GT.6) THEN
        YC(7)= C*R1B2-S*R2B2
        YC(8)= C*Y1B2-S*Y2B2
        YC(9)= S*R1B2+C*R2B2
        YC(10)=S*Y1B2+C*Y2B2
      END IF
C
      RETURN
      END
C
C=======================================================================
C