C
C *** RICKER WAVE IN 2-D IN-PLANE LAYERED HALF-SPACE
C     USING EXACT GREEN FUNCTION AND RECIPROCITY RELATION
C
      PROGRAM MSVHASK
C
      PARAMETER(MOM=100,ML=1,MOB=6)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP(ML),VS(ML),FW0(MOM,MOB),FU0(MOM,MOB),CPX
     *          ,FSZ1,FCZ1,FSX1,FCX1,FSZ2,FCZ2,FSX2,FCX2,AA,CC
      DIMENSION DNS(ML),THK(ML),XOB(MOB)
      COMMON / CN / C1,N1,C2,N2 /FDAT/ OM,NL,QX,QZ
      EXTERNAL FNCA, FNCB
      CHARACTER FNAME*15
C
      WRITE(*,1000)
 1000 FORMAT(' MSVHASK: P-SV WAVE IN 2-D IN-PLANE HOMOGENEOUS'/
     *       '   HALF-SPACE USING EXACT GREEN FUNCTION AND'/
     *       '   RECIPLOCITY RELATION (SP, OB IN ANY LAYER),')
C
      CALL DATIN2(DOM,NOM,ML,NL1,DNS,VP,VS,THK,LS,SDEP,LO,ODEP)
      NL=NL1
      IF(NL.NE.1)THEN
        WRITE(*,*)' NL.NE.1 !!!!!!!!!!'
        STOP
      END IF
C
   30 WRITE(*,1410)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : ENTER NOB : MOB=',MOB
        GOTO 30
      END IF
      IF(NOB.EQ.0)GOTO 30
C
      DO 100 I=1,NOB
        WRITE(*,1420)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
 1410 FORMAT(/' ENTER NUMBER OF OB : ',$)
 1420 FORMAT( ' ENTER X-AXIS OF OB (',I2,') : ',$) 
C
      TAU=3.1415927D0
      VR=1000.D0
      DV=1.D0
      FW=3141.5927D0
C
      WRITE(*,1521)TAU,DV,VR,FW
 1521 FORMAT(/' TAU, DV, VR, FW =',2F10.5,2F10.3) 
C
      WRITE(*,1550)
      READ(*,*)DA
      DA=DA*0.0174532D0
C
 1550 FORMAT(/' ENTER DELTA (DEG) : ',$)
C
C *** GREEN FUNC ***
C
      CPX=DCMPLX(0.D0,1.D0)
      AA=DNS(1)*VS(1)*VS(1)
      TW=FW/VR
      SID=DSIN(2.D0*DA)
      COD=DCOS(2.D0*DA)
C
      DO 200 IOB=1,NOB
C
        XX= -XOB(IOB)
C
        DO 210 IOM=1,NOM
C
          OM=DOM*DFLOAT(IOM)
          XW=OM*TW/2.D0
          XT=OM*TAU/2.D0
          SW=DSIN(XW)/XW
          ST=DSIN(XT)/XT
          CC=AA*DV*FW*SW*ST*CDEXP(-CPX*(XW+XT))/(CPX*OM)
C
          WNB=OM/C2
          WNA=WNB/(1000.D0*N1)
          WNC=OM/C1
C
          IF((LO.LT.LS).OR.((LO.EQ.LS).AND.(ODEP.LE.SDEP)))THEN
C
            CALL NEWTNN(FNCA,WNA,WNB,N1,FSZ1,FCZ1,FSX1,FCX1
     *                 ,XX,DNS,VP,VS,THK,ML,IER,LS,SDEP,LO,ODEP)
            IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
            CALL NEWTNN(FNCA,WNB,WNC,N2,FSZ2,FCZ2,FSX2,FCX2
     *                 ,XX,DNS,VP,VS,THK,ML,IER,LS,SDEP,LO,ODEP)
            IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
C
          ELSE
C
            CALL NEWTNN(FNCB,WNA,WNB,N1,FSZ1,FCZ1,FSX1,FCX1
     *                 ,XX,DNS,VP,VS,THK,ML,IER,LO,ODEP,LS,SDEP)
            IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
            CALL NEWTNN(FNCB,WNB,WNC,N2,FSZ2,FCZ2,FSX2,FCX2
     *                 ,XX,DNS,VP,VS,THK,ML,IER,LO,ODEP,LS,SDEP)
            IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
C
          END IF
C
          FW0(IOM,IOB)=CC*(SID*(FSZ1+FSZ2)-COD*(FCZ1+FCZ2))
          FU0(IOM,IOB)=CC*(SID*(FSX1+FSX2)-COD*(FCX1+FCX2))
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** OUTPUT OF DATA ***
C
        WRITE(*,2400)DOM,NOM
        DO 300 IOB=1,NOB
        WRITE(*,2410)XOB(IOB)
          DO 310 K=1,NOM,4
            WRITE(*,2420)(FW0(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
        DO 320 IOB=1,NOB
          WRITE(*,2415)XOB(IOB)
          DO 330 K=1,NOM,4
            WRITE(*,2420)(FU0(K+J,IOB),J=0,3)
  330     CONTINUE
  320   CONTINUE
C
 2400 FORMAT(/' DOM,NOM=',F6.2,I5)
 2410 FORMAT(/' X=',F10.2,'; FW0*4')
 2415 FORMAT(/' X=',F10.2,'; FU0*4')
 2420 FORMAT(1X,8D10.3)
C
C *** OUTPUT OF DATA FILE ***
C
      WRITE(*,*)
   50 WRITE(*,2310)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 50
      WRITE(*,2320)
 2310 FORMAT( ' Enter file name : ',$)
 2320 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VV roop')
C      
      OPEN(8,FILE=FNAME,STATUS='NEW')
        WRITE(8,2500)DOM,NOM
        WRITE(8,2510)
        DO 500 IOB=1,NOB
          DO 510 K=1,NOM,5
            WRITE(8,2520)(FW0(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
        DO 520 IOB=1,NOB
          DO 530 K=1,NOM,5
            WRITE(8,2520)(FU0(K+J,IOB),J=0,4)
  530     CONTINUE
  520   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F15.8,I15)
 2510 FORMAT(///)
 2520 FORMAT(1X,10F8.5)
C
      STOP
      END
C
      SUBROUTINE DATIN2(DOM,NOM,ML,NL,DNS,VP,VS,THK,LS,SDEP,LO,ODEP)
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30
      COMPLEX*16 VP(ML),VS(ML)
      DIMENSION DNS(ML),THK(ML)
      COMMON / CN / C1,N1,C2,N2
C
C *** SOIL DATA INPUT FROM FILE DATA ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1000 FORMAT( ' From SUBROUTINE DATAIN :')
 1010 FORMAT(/' Enter Soil Data File Name ',$)

      OPEN (7,FILE=FNAME,STATUS='OLD')
      READ(7,1200)NL
      DO 100 IL=1,NL
        READ(7,1210)DNS(IL),VPP,VSS,THK(IL),QF
        VP(IL)=VPP*DCMPLX(1.D0,0.5D0/QF)
        VS(IL)=VSS*DCMPLX(1.D0,0.5D0/QF)
  100 CONTINUE
      CLOSE(7)
 1200 FORMAT(I10)
 1210 FORMAT(5F10.0)
C
C *** DATA ECHO ***
C
      WRITE(*,2000)NL
      DO 200 IL=1,NL
        WRITE(*,2010)IL,DNS(IL),VP(IL),VS(IL),THK(IL)
  200 CONTINUE
 2000 FORMAT(/' Total Layers Number            =',I5/
     *       /' Layer Num., Density , Vp(real), Vp(imag)'
     *        ', Vs(real), Vs(imag),Thickness')
 2010 FORMAT(I5,5X,6F10.3)
C
   20 WRITE(*,3000)
      READ(*,*)DOM,NOM
 3000 FORMAT(/' Enter Delta Omega and Omega Number        ? :',$)
C
      WRITE(*,4000)
      READ(*,*)C1
      WRITE(*,4010)
      READ(*,*)C2
      WRITE(*,4020)
      READ(*,*)N1
      WRITE(*,4030)
      READ(*,*)N2
C
 4000 FORMAT(/' Enter the lowest phase vel. C1        : ',$)
 4010 FORMAT( ' Enter the second phase vel. C2(>C1)   : ',$)
 4020 FORMAT( ' Enter div. num. for int. (0-OM/C2)    : ',$)
 4030 FORMAT( ' Enter div. num. for int. (OM/C2-OM/C1): ',$)
C
C *** ENTER OF SOURCE DEPTH IN ANY LAYER ***
C
   30 WRITE(*,5000)NL
C      READ(*,*)LS
C      IF(LS.EQ.0)GOTO 30
C      IF(LS.GT.NL)GOTO 30
C      WRITE(*,5010)THK(LS)
C      READ(*,*)SDEP
      LO=1
      WRITE(*,5010)THK(LO)
      READ(*,*)ODEP
 5000 FORMAT(/' Total layers number =',I5,
     *       /' Enter source layer number (from top layer) : ',$)
 5010 FORMAT( ' Source layer thickness =',F10.4,
     *       /' Enter source depth from upper boundary : ',$)
C
C *** ENTER OF OP DEPTH IN ANY LAYER ***
C
   40 WRITE(*,5100)NL
C      READ(*,*)LO
C      IF(LO.EQ.0)GOTO 40
C      IF(LO.GT.NL)GOTO 40
C      WRITE(*,5110)THK(LO)
C      READ(*,*)ODEP
      LS=1
      WRITE(*,5110)THK(LS)
      READ(*,*)SDEP
 5100 FORMAT(/' Total layers number =',I5,
     *       /' Enter OP layer number (from top layer) : ',$)
 5110 FORMAT( ' Source layer thickness =',F10.4,
     *       /' Enter OP depth from upper boundary : ',$)
C
      WRITE(*,6000)
      READ(*,*)IY
      IF(IY.NE.0)GOTO 20
 6000 FORMAT(/' All data OK? (Yes=0 enter) : ',$)
C
      WRITE(*,*)
C
      RETURN
      END
C
      SUBROUTINE NEWTNN(FNC,A,B,N,SZ,CZ,SX,CX,XX,DNS,VP,VS
     *                 ,THK,ML,IER,L1,DEP1,L2,DEP2)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP(ML),VS(ML),SZ,CZ,SX,CX
     *          ,FSZ,FCZ,FSX,FCX,YSZ,YCZ,YSX,YCX,QSZ,QCZ,QSX,QCX
      DIMENSION DNS(ML),THK(ML)
C
      IF(N.LT.1)THEN
        IER=1000
        SZ=(0.D0,0.D0)
        CZ=(0.D0,0.D0)
        SX=(0.D0,0.D0)
        CX=(0.D0,0.D0)
        GO TO 9999
      END IF
      IF(A.EQ.B)THEN
        IER=2000
        SZ=(0.D0,0.D0)
        CZ=(0.D0,0.D0)
        SX=(0.D0,0.D0)
        CX=(0.D0,0.D0)
        GO TO 9999
      END IF
C
      H=(B-A)/N
      CALL FNC(A,YSZ,YCZ,YSX,YCX,XX,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2)
      QSZ=(0.D0,0.D0)
      QCZ=(0.D0,0.D0)
      QSX=(0.D0,0.D0)
      QCX=(0.D0,0.D0)
C
      DO  100 I=1,N
      X=A+H*I
      CALL FNC(X,FSZ,FCZ,FSX,FCX,XX,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2)
      QSZ  = QSZ +H*(YSZ+FSZ)/2.D0
      QCZ  = QCZ +H*(YCZ+FCZ)/2.D0
      QSX  = QSX +H*(YSX+FSX)/2.D0
      QCX  = QCX +H*(YCX+FCX)/2.D0
      YSZ  = FSZ
      YCZ  = FCZ
      YSX  = FSX
      YCX  = FCX
  100 CONTINUE
C
      SZ = QSZ
      CZ = QCZ
      SX = QSX
      CX = QCX
      IER=0
C
 9999 RETURN
C
      END
C
      SUBROUTINE FNCA(WN,FSZ,FCZ,FSX,FCX,XX,DNS,VP,VS
     *               ,THK,ML,L1,DEP1,L2,DEP2)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),THK(ML)
      COMPLEX*16 FSZ,FCZ,FSX,FCX,VP(ML),VS(ML)
     *          ,AA(4,4),SA(4,4),SB(4,4),CPX,CCC
     *          ,WP2,WS2,CMP,CMS,BM2,PM,QM,RBM,RAM
     *          ,CQ,SQ,STQ,SDQ,CP,SP,STP,SDP,GM
     *          ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
     *          ,GN,RAN,RBN,RL,RK,RN,RM,GS2,RS2,HS2,SS2
     *          ,FWZ0,FWX0,FUZ0,FUX0,FR,EPD,ESD,CMPN,CMSN
     *          ,B11,B12,B13,B14,B21,B22,B23,B24
     *          ,AS2(4,4),AR1(4,4),FWZR,FWXR,FUZR,FUXR
     *          ,FWZRZ,FWXRZ,FUZRZ,FUXRZ
     *          ,UXXX,UZXX,UXZX,UZZX,UXXZ,UZXZ,UXZZ,UZZZ
      COMMON /FDAT/ OM,NL,QX,QZ
C
      WN2=WN*WN
      OM2=OM*OM
      CV=OM/WN
      CV2=CV*CV
      CPX=DCMPLX(0.D0,1.D0)
C
C *** AA MATRIX ***
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J) =(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
      DO 150 IL=NL-1,1,-1
C
        DO 160 I=1,4
          DO 170 J=1,4
            SA(I,J) =AA(I,J)
  170     CONTINUE
  160   CONTINUE
C
        BM2=VS(IL)*VS(IL)
        WP2=OM2/(VP(IL)*VP(IL))
        WS2=OM2/BM2
C
        CMP=CDSQRT(WP2-WN2)
        CMS=CDSQRT(WS2-WN2)
        IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
        IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
        RAM=CMP/WN
        PM=CMP*THK(IL)
        CP=CDCOS(PM)
        SP=CDSIN(PM)
        STP=RAM*SP
        SDP=SP/RAM
C
        RBM=CMS/WN
        QM=CMS*THK(IL)
        CQ=CDCOS(QM)
        SQ=CDSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
        GM=2.D0*BM2/CV2
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CV2
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= C11*CP-C12*CQ
        SB(1,2)= (C12*SDP+C11*STQ)*CPX
        SB(1,3)=-C15*D11
        SB(1,4)= C15*D12*CPX
        SB(2,1)=-(C11*STP+C12*SDQ)*CPX
        SB(2,2)=-C12*CP+C11*CQ
        SB(2,3)= C15*D13*CPX
        SB(2,4)= SB(1,3)
        SB(3,1)= C16*C11*C12*D11
        SB(3,2)= C16*D14*CPX
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= C16*D15*CPX
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
        DO 180 I=1,4
          DO 190 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 195 K=1,4
              AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
  195       CONTINUE
  190     CONTINUE
  180   CONTINUE
C
  150 CONTINUE
C
C *** RL-RM AND FR ***
C
      BM2=VS(NL)*VS(NL)
      WP2=OM2/(VP(NL)*VP(NL))
      WS2=OM2/BM2
      CMPN=CDSQRT(WP2-WN2)
      CMSN=CDSQRT(WS2-WN2)
      IF(DIMAG(CMPN).GT.0.D0)CMPN=-CMPN
      IF(DIMAG(CMSN).GT.0.D0)CMSN=-CMSN
C
      GN=2.D0*BM2/CV2
      RAN=CMPN/WN
      RBN=CMSN/WN
C
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DNS(NL)*CV2)
      B13=-RAN*B14
      B21=-B12
      B22=GN*RBN
      B23=B14
      B24=RBN*B14
C
      RL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
      RK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
      RN=B21*AA(1,1)+B22*AA(2,1)+B23*AA(3,1)+B24*AA(4,1)
      RM=B21*AA(1,2)+B22*AA(2,2)+B23*AA(3,2)+B24*AA(4,2)
C
      FR=RN*RK-RL*RM
C
C *** SP AND OB ***
C
      IF(L1.EQ.NL)THEN
C
        EPD=CDEXP(-CPX*CMPN*DEP1)
        ESD=CDEXP(-CPX*CMSN*DEP1)
C
        RS2= B14*EPD
        GS2=-RAN*RS2
        HS2= B14*ESD
        SS2= RBN*HS2
C
      ELSE
C
        CALL FNC22(AS2,L1,DEP1,OM2,CV2,WN,DNS,VP,VS,THK,ML,NL)
C
        GS2=B11*AS2(1,3)+B12*AS2(2,3)+B13*AS2(3,3)+B14*AS2(4,3)
        RS2=B11*AS2(1,4)+B12*AS2(2,4)+B13*AS2(3,4)+B14*AS2(4,4)
        HS2=B21*AS2(1,3)+B22*AS2(2,3)+B23*AS2(3,3)+B24*AS2(4,3)
        SS2=B21*AS2(1,4)+B22*AS2(2,4)+B23*AS2(3,4)+B24*AS2(4,4)
C
      END IF

      FWZ0=GS2*RN-HS2*RL
      FWX0=RS2*RN-SS2*RL
      FUZ0=GS2*RM-HS2*RK
      FUX0=RS2*RM-SS2*RK
C
      CALL FNC21(AR1,L2,DEP2,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
      FWZR=AR1(2,2)*FWZ0-AR1(2,1)*FUZ0
      FWXR=AR1(2,2)*FWX0-AR1(2,1)*FUX0
      FUZR=AR1(1,1)*FUZ0-AR1(1,2)*FWZ0
      FUXR=AR1(1,1)*FUX0-AR1(1,2)*FWX0
C
      CALL FNC21Z(AR1,L2,DEP2,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
      FWZRZ=AR1(2,2)*FWZ0-AR1(2,1)*FUZ0
      FWXRZ=AR1(2,2)*FWX0-AR1(2,1)*FUX0
      FUZRZ=AR1(1,1)*FUZ0-AR1(1,2)*FWZ0
      FUXRZ=AR1(1,1)*FUX0-AR1(1,2)*FWX0
C
      WNX=WN*XX
      CKX=DCOS(WNX)
      SKX=DSIN(WNX)
      PAI=3.1415927D0
C
      CCC=1.D0/(PAI*FR)
      UZXX=-CCC*FWXR*CKX
      UXZX= CCC*FUZR*CKX
      UXXX=-CPX*CCC*FUXR*SKX
      UZZX= CPX*CCC*FWZR*SKX
      CCC=CCC/WN
      UZXZ=-CCC*FWXRZ*SKX
      UXZZ= CCC*FUZRZ*SKX
      UXXZ= CPX*CCC*FUXRZ*CKX
      UZZZ=-CPX*CCC*FWZRZ*CKX
C
      FSX=UXXX-UZXZ
      FCX=UXXZ+UZXX
      FSZ=UXZX-UZZZ
      FCZ=UXZZ+UZZX
C
      RETURN
      END
C
      SUBROUTINE FNCB(WN,FSZ,FCZ,FSX,FCX,XX,DNS,VP,VS
     *               ,THK,ML,L1,DEP1,L2,DEP2)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),THK(ML)
      COMPLEX*16 FSZ,FCZ,FSX,FCX,VP(ML),VS(ML)
     *          ,AA(4,4),SA(4,4),SB(4,4),CPX,CCC
     *          ,WP2,WS2,CMP,CMS,BM2,PM,QM,RBM,RAM
     *          ,CQ,SQ,STQ,SDQ,CP,SP,STP,SDP,GM
     *          ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
     *          ,GN,RAN,RBN,RL,RK,RN,RM,GS2,RS2,HS2,SS2
     *          ,FWZ0,FWX0,FUZ0,FUX0,FR,EPD,ESD,CMPN,CMSN
     *          ,B11,B12,B13,B14,B21,B22,B23,B24
     *          ,AS2(4,4),AR1(4,4),FWZR,FWXR,FUZR,FUXR
     *          ,FWZRZ,FWXRZ,FUZRZ,FUXRZ
     *          ,UXXX,UZXX,UXZX,UZZX,UXXZ,UZXZ,UXZZ,UZZZ
      COMMON /FDAT/ OM,NL,QX,QZ
C
      WN2=WN*WN
      OM2=OM*OM
      CV=OM/WN
      CV2=CV*CV
      CPX=DCMPLX(0.D0,1.D0)
C
C *** AA MATRIX ***
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J) =(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
      DO 150 IL=NL-1,1,-1
C
        DO 160 I=1,4
          DO 170 J=1,4
            SA(I,J) =AA(I,J)
  170     CONTINUE
  160   CONTINUE
C
        BM2=VS(IL)*VS(IL)
        WP2=OM2/(VP(IL)*VP(IL))
        WS2=OM2/BM2
C
        CMP=CDSQRT(WP2-WN2)
        CMS=CDSQRT(WS2-WN2)
        IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
        IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
        RAM=CMP/WN
        PM=CMP*THK(IL)
        CP=CDCOS(PM)
        SP=CDSIN(PM)
        STP=RAM*SP
        SDP=SP/RAM
C
        RBM=CMS/WN
        QM=CMS*THK(IL)
        CQ=CDCOS(QM)
        SQ=CDSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
        GM=2.D0*BM2/CV2
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CV2
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= C11*CP-C12*CQ
        SB(1,2)= (C12*SDP+C11*STQ)*CPX
        SB(1,3)=-C15*D11
        SB(1,4)= C15*D12*CPX
        SB(2,1)=-(C11*STP+C12*SDQ)*CPX
        SB(2,2)=-C12*CP+C11*CQ
        SB(2,3)= C15*D13*CPX
        SB(2,4)= SB(1,3)
        SB(3,1)= C16*C11*C12*D11
        SB(3,2)= C16*D14*CPX
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= C16*D15*CPX
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
        DO 180 I=1,4
          DO 190 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 195 K=1,4
              AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
  195       CONTINUE
  190     CONTINUE
  180   CONTINUE
C
  150 CONTINUE
C
C *** RL-RM AND FR ***
C
      BM2=VS(NL)*VS(NL)
      WP2=OM2/(VP(NL)*VP(NL))
      WS2=OM2/BM2
      CMPN=CDSQRT(WP2-WN2)
      CMSN=CDSQRT(WS2-WN2)
      IF(DIMAG(CMPN).GT.0.D0)CMPN=-CMPN
      IF(DIMAG(CMSN).GT.0.D0)CMSN=-CMSN
C
      GN=2.D0*BM2/CV2
      RAN=CMPN/WN
      RBN=CMSN/WN
C
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DNS(NL)*CV2)
      B13=-RAN*B14
      B21=-B12
      B22=GN*RBN
      B23=B14
      B24=RBN*B14
C
      RL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
      RK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
      RN=B21*AA(1,1)+B22*AA(2,1)+B23*AA(3,1)+B24*AA(4,1)
      RM=B21*AA(1,2)+B22*AA(2,2)+B23*AA(3,2)+B24*AA(4,2)
C
      FR=RN*RK-RL*RM
C
C *** SP AND OB ***
C
      IF(L1.EQ.NL)THEN
C
        EPD=CDEXP(-CPX*CMPN*DEP1)
        ESD=CDEXP(-CPX*CMSN*DEP1)
C
        RS2= B14*EPD
        GS2=-RAN*RS2
        HS2= B14*ESD
        SS2= RBN*HS2
C
      ELSE
C
        CALL FNC22(AS2,L1,DEP1,OM2,CV2,WN,DNS,VP,VS,THK,ML,NL)
C
        GS2=B11*AS2(1,3)+B12*AS2(2,3)+B13*AS2(3,3)+B14*AS2(4,3)
        RS2=B11*AS2(1,4)+B12*AS2(2,4)+B13*AS2(3,4)+B14*AS2(4,4)
        HS2=B21*AS2(1,3)+B22*AS2(2,3)+B23*AS2(3,3)+B24*AS2(4,3)
        SS2=B21*AS2(1,4)+B22*AS2(2,4)+B23*AS2(3,4)+B24*AS2(4,4)
C
      END IF

      FWZ0=GS2*RN-HS2*RL
      FWX0=RS2*RN-SS2*RL
      FUZ0=GS2*RM-HS2*RK
      FUX0=RS2*RM-SS2*RK
C
      CALL FNC21(AR1,L2,DEP2,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
      FWZR=AR1(2,2)*FWZ0-AR1(2,1)*FUZ0
      FUZR=AR1(2,2)*FWX0-AR1(2,1)*FUX0
      FWXR=AR1(1,1)*FUZ0-AR1(1,2)*FWZ0
      FUXR=AR1(1,1)*FUX0-AR1(1,2)*FWX0
C
      IF(L1.EQ.NL)THEN
C
        EPD=-CPX*CMPN*EPD
        ESD=-CPX*CMSN*ESD
C
        RS2= B14*EPD
        GS2=-RAN*RS2
        HS2= B14*ESD
        SS2= RBN*HS2
C
      ELSE
C
        CALL FNC22(AS2,L1,DEP1,OM2,CV2,WN,DNS,VP,VS,THK,ML,NL)
C
        GS2=B11*AS2(1,3)+B12*AS2(2,3)+B13*AS2(3,3)+B14*AS2(4,3)
        RS2=B11*AS2(1,4)+B12*AS2(2,4)+B13*AS2(3,4)+B14*AS2(4,4)
        HS2=B21*AS2(1,3)+B22*AS2(2,3)+B23*AS2(3,3)+B24*AS2(4,3)
        SS2=B21*AS2(1,4)+B22*AS2(2,4)+B23*AS2(3,4)+B24*AS2(4,4)
C
      END IF

      FWZ0=GS2*RN-HS2*RL
      FWX0=RS2*RN-SS2*RL
      FUZ0=GS2*RM-HS2*RK
      FUX0=RS2*RM-SS2*RK
C
      CALL FNC21(AR1,L2,DEP2,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
      FWZRZ=AR1(2,2)*FWZ0-AR1(2,1)*FUZ0
      FUZRZ=AR1(2,2)*FWX0-AR1(2,1)*FUX0
      FWXRZ=AR1(1,1)*FUZ0-AR1(1,2)*FWZ0
      FUXRZ=AR1(1,1)*FUX0-AR1(1,2)*FWX0
C
      WNX=WN*XX
      CKX=DCOS(WNX)
      SKX=DSIN(WNX)
      PAI=3.1415927D0
C
      CCC=1.D0/(PAI*FR)
      UZXX=-CCC*FWXR*CKX
      UXZX= CCC*FUZR*CKX
      UXXX=-CPX*CCC*FUXR*SKX
      UZZX= CPX*CCC*FWZR*SKX
      CCC=CCC/WN
      UZXZ=-CCC*FWXRZ*SKX
      UXZZ= CCC*FUZRZ*SKX
      UXXZ= CPX*CCC*FUXRZ*CKX
      UZZZ=-CPX*CCC*FWZRZ*CKX
C
      FSX=UXXX-UZXZ
      FCX=UXXZ+UZXX
      FSZ=UXZX-UZZZ
      FCZ=UXZZ+UZZX
C
      RETURN
      END
C
      SUBROUTINE FNC21(AA,LL,DEP,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),THK(ML)
      COMPLEX*16 VP(ML),VS(ML)
     *          ,AA(4,4),SA(4,4),SB(4,4),CPX
     *          ,BM2,WP2,WS2,CMP,CMS,PM,QM,RBM,RAM
     *          ,CQ,SQ,STQ,SDQ,CP,SP,STP,SDP,GM
     *          ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
C
      CPX=DCMPLX(0.D0,1.D0)
      WN2=WN*WN
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J) =(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
      DO 150 IL=LL,1,-1
C
        DO 160 I=1,4
          DO 170 J=1,4
            SA(I,J) =AA(I,J)
  170     CONTINUE
  160   CONTINUE
C
        BM2=VS(IL)*VS(IL)
        WP2=OM2/(VP(IL)*VP(IL))
        WS2=OM2/BM2
C
        CMP=CDSQRT(WP2-WN2)
        CMS=CDSQRT(WS2-WN2)
        IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
        IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
        RAM=CMP/WN
        PM=CMP*THK(IL)
        IF(IL.EQ.LL)PM=CMP*DEP
        CP=CDCOS(PM)
        SP=CDSIN(PM)
        STP=RAM*SP
        SDP=SP/RAM
C
        RBM=CMS/WN
        QM=CMS*THK(IL)
        IF(IL.EQ.LL)QM=CMS*DEP
        CQ=CDCOS(QM)
        SQ=CDSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
        GM=2.D0*BM2/CV2
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CV2
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= C11*CP-C12*CQ
        SB(1,2)= (C12*SDP+C11*STQ)*CPX
        SB(1,3)=-C15*D11
        SB(1,4)= C15*D12*CPX
        SB(2,1)=-(C11*STP+C12*SDQ)*CPX
        SB(2,2)=-C12*CP+C11*CQ
        SB(2,3)= C15*D13*CPX
        SB(2,4)= SB(1,3)
        SB(3,1)= C16*C11*C12*D11
        SB(3,2)= C16*D14*CPX
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= C16*D15*CPX
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
        DO 180 I=1,4
          DO 190 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 195 K=1,4
              AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
  195       CONTINUE
  190     CONTINUE
  180   CONTINUE
C
  150 CONTINUE
C
      RETURN
      END
C
C
      SUBROUTINE FNC21Z(AA,LL,DEP,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),THK(ML)
      COMPLEX*16 VP(ML),VS(ML)
     *          ,AA(4,4),SA(4,4),SB(4,4),CPX
     *          ,BM2,WP2,WS2,CMP,CMS,PM,QM,RBM,RAM
     *          ,CQ,SQ,STQ,SDQ,CP,SP,STP,SDP,GM
     *          ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
C
      CPX=DCMPLX(0.D0,1.D0)
      WN2=WN*WN
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J) =(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
      DO 150 IL=LL,1,-1
C
        DO 160 I=1,4
          DO 170 J=1,4
            SA(I,J) =AA(I,J)
  170     CONTINUE
  160   CONTINUE
C
        BM2=VS(IL)*VS(IL)
        WP2=OM2/(VP(IL)*VP(IL))
        WS2=OM2/BM2
C
        CMP=CDSQRT(WP2-WN2)
        CMS=CDSQRT(WS2-WN2)
        IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
        IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
        RAM=CMP/WN
        PM=CMP*THK(IL)
        IF(IL.EQ.LL)THEN
          PM=CMP*DEP
          CP=-CMP*CDSIN(PM)
          SP= CMP*CDCOS(PM)
        ELSE
          CP=CDCOS(PM)
          SP=CDSIN(PM)
        END IF
        STP=RAM*SP
        SDP=SP/RAM
C
        RBM=CMS/WN
        QM=CMS*THK(IL)
        IF(IL.EQ.LL)THEN
          QM=CMS*DEP
          CQ=-CMS*CDSIN(QM)
          SQ= CMS*CDCOS(QM)
        ELSE
          CQ=CDCOS(QM)
          SQ=CDSIN(QM)
        END IF
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
        GM=2.D0*BM2/CV2
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CV2
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= C11*CP-C12*CQ
        SB(1,2)= (C12*SDP+C11*STQ)*CPX
        SB(1,3)=-C15*D11
        SB(1,4)= C15*D12*CPX
        SB(2,1)=-(C11*STP+C12*SDQ)*CPX
        SB(2,2)=-C12*CP+C11*CQ
        SB(2,3)= C15*D13*CPX
        SB(2,4)= SB(1,3)
        SB(3,1)= C16*C11*C12*D11
        SB(3,2)= C16*D14*CPX
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= C16*D15*CPX
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
        DO 180 I=1,4
          DO 190 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 195 K=1,4
              AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
  195       CONTINUE
  190     CONTINUE
  180   CONTINUE
C
  150 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE FNC22(AA,LL,DEP,OM2,CV2,WN,DNS,VP,VS,THK,ML,NL)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),THK(ML)
      COMPLEX*16 VP(ML),VS(ML)
     *          ,AA(4,4),SA(4,4),SB(4,4),CPX
     *          ,BM2,WP2,WS2,CMP,CMS,PM,QM,RBM,RAM
     *          ,CQ,SQ,STQ,SDQ,CP,SP,STP,SDP,GM
     *          ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
C
      CPX=DCMPLX(0.D0,1.D0)
      WN2=WN*WN
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J) =(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
      DO 150 IL=NL-1,LL,-1
C
        DO 160 I=1,4
          DO 170 J=1,4
            SA(I,J) =AA(I,J)
  170     CONTINUE
  160   CONTINUE
C
        BM2=VS(IL)*VS(IL)
        WP2=OM2/(VP(IL)*VP(IL))
        WS2=OM2/BM2
C
        CMP=CDSQRT(WP2-WN2)
        CMS=CDSQRT(WS2-WN2)
        IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
        IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
        RAM=CMP/WN
        PM=CMP*THK(IL)
        IF(IL.EQ.LL)PM=CMP*(THK(IL)-DEP)
        CP=CDCOS(PM)
        SP=CDSIN(PM)
        STP=RAM*SP
        SDP=SP/RAM
C
        RBM=CMS/WN
        QM=CMS*THK(IL)
        IF(IL.EQ.LL)QM=CMS*(THK(IL)-DEP)
        CQ=CDCOS(QM)
        SQ=CDSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
        GM=2.D0*BM2/CV2
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CV2
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= C11*CP-C12*CQ
        SB(1,2)= (C12*SDP+C11*STQ)*CPX
        SB(1,3)=-C15*D11
        SB(1,4)= C15*D12*CPX
        SB(2,1)=-(C11*STP+C12*SDQ)*CPX
        SB(2,2)=-C12*CP+C11*CQ
        SB(2,3)= C15*D13*CPX
        SB(2,4)= SB(1,3)
        SB(3,1)= C16*C11*C12*D11
        SB(3,2)= C16*D14*CPX
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= C16*D15*CPX
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
        DO 180 I=1,4
          DO 190 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 195 K=1,4
              AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
  195       CONTINUE
  190     CONTINUE
  180   CONTINUE
C
  150 CONTINUE
C
      RETURN
      END
