C SUBROUTINE FILE 'TRANS.FOR' TO TRANSFORM THE COMPUTED QUANTITIES AT AN
C INTERFACE
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C THIS FILE CONSISTS OF THE FOLLOWING SUBROUTINES:
C     TRANS...SUBROUTINE TRANSFORMING THE COMPUTED QUANTITIES ACROSS A
C             CURVED INTERFACE, SEE C.R.T.5.9.
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.R.T.5.5.4.
C
C=======================================================================
C
      SUBROUTINE TRANS(YL,Y,YY,IY,KODNEW,ICBNEW,IEND)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12),ICBNEW,IEND
C
C THIS SUBROUTINE TRANSFORMS THE COMPUTED QUANTITIES ACROSS A CURVED
C INTERFACE, SEE 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 ERROR MESSAGES:
C     591...  WRONG INTERFACE:
C             THIS ERROR SHOULD NOT APPEAR.  CONTACT THE AUTHORS.
C
C DATE: 1992, DECEMBER 19
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS:
      REAL G(12),GAMMA(18),GSQRD,FAUX(10)
      REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL
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 /AUXMOD/ BUT TO SHARE THE MEMORY.
      COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C     FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS
C             DERIVATIVES.
C     UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'.
C     UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE
C             'MODEL.FOR'.
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     NAMPL1..NUMBER OF THE APMLITUDE COEFFICIENTS OF THE INCIDENT WAVE
C             (5.9.1A).
C     NAMPL2..NUMBER OF THE APMLITUDE 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 )... AMPLIDUDE COEFFICIENTS (5.9.6):
C     ( R1A2,Y1A2,R2A2,Y2A2 )
C     ( R1B2,Y1B2,R2B2,Y2B2 )               INITIAL
C             PART:          COMPONENT:     POLARISATION:  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
        PAUSE 'ERROR 591 IN TRANS: WRONG INTERFACE'
        STOP
      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. (SEE C.R.T.5.9.7)
        ND=0
      ELSE
C       EPSILON=-1. (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(NAMPL1.GT.4) THEN
        R1B2=RR*R1B1-YR*Y1B1
        Y1B2=YR*R1B1+RR*Y1B1
      END IF
      IF(ICBNEW.GE.0) THEN
C       OUTGOING P WAVE:
        Y(28)=R1A2
        Y(29)=Y1A2
        IF(NAMPL2.GT.2) THEN
          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(NAMPL1.GT.4) THEN
          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
      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.R.T.5.5.4.
C
C INPUT:
C     KSTORE..SPECIFIES WHETHER THE CONVERSION COEFFICIENTS ARE TO BE
C             CONSIDERED (SEE C.R.T.5.5.4 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.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 VELOC,KOOR,METRIC,SRFC2,PARM2,SMVPRD,COEF50
      INTEGER KOOR
C     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: 1990, OCTOBER 1
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS:
      REAL G(12),GAMMA(18),GSQRD,FAUX(10)
      REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL
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 /AUXMOD/ BUT TO SHARE THE MEMORY.
      COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE PACKAGE 'METRIC'.
C     FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS
C             DERIVATIVES.
C     UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'.
C     UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE
C             'MODEL.FOR'.
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     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 )... AMPLIDUDE 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:     POLARISATION:  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 AMPLIDUDE COEFFICIENTS IN RAY-CENTRED
C     ( RB2,YB2 )    COORDINATE SYSTEM OF THE R/T WAVE
C             PART:          POLARISATION:  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(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. (SEE C.R.T.5.9.7)
        ND=0
      ELSE
C       EPSILON=-1. (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
