C     ***************************************************
C     *  PROGRAM BMAIN (OUTPUT OF SOURCE INPUT MOTION)  *
C     ***************************************************
C
      PROGRAM BMAIN
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (ND1=100,ND2=20,ND3=60,ND4=5)
      COMPLEX*16 V2(ND1,ND3,ND4),P2(ND1,ND3,ND4),VV(ND2,ND1),PP(ND2,ND1)
     *          ,V1(ND2,ND3,ND4),V3(ND2,ND3,ND4),VF(ND2,ND4)
      DIMENSION SZ2(ND1),ES1(ND2),NES2(ND2),ES3(ND2),NCF(10)
     *         ,FX0(ND4),FZ0(ND4),DA(ND4)
C
C      CALL CLOCK
C
      WRITE(6,6000)ND1,ND2,ND3,ND4
C
      CALL INPUT(SZ2,ES1,NES2,ES3,V2,P2,ND1,ND2,ND3,ND4
     *          ,FX0,FZ0,DA,NCF)
C
      CALL CALVV(SZ2,ES1,ES3,ND1,ND2,ND3,V1,V3,VF,V2,P2,VV,PP
     *          ,FX0,FZ0,DA,NCF)
C
      CALL OUTPUT(V1,V2,P2,V3,ES1,NES2,ES3,ND1,ND2,ND3)
C
C      CALL CLOCK(TIME)
C      WRITE(6,6010)TIME
C
      STOP
C
 6000 FORMAT(//10X,'DIMENSION DATA'//
     *         10X,'ND1(2*Elem.Num) = ',I5,', ND2(2*Out Num) = ',I5/
     *         10X,'ND3( Freq.Num ) = ',I5,', ND4( Flt Num ) = ',I5/)
 6010 FORMAT(//10X,'CPU TIME (sec)=',F15.7/)
      END
C     ***********************
C     *  SUBROUTINE INPUT   *
C     ***********************
C
      SUBROUTINE INPUT(SZ2,ES1,NES2,ES3,V2,P2
     *                ,ND1,ND2,ND3,ND4,FX0,FZ0,DA,NCF)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 V2(ND1,ND3,*),P2(ND1,ND3,*)
      DIMENSION SZ2(*),ES1(*),NES2(*),ES3(*),NCF(10)
     *         ,FX0(*),FZ0(*),DA(*)
      COMMON /FREQ/ FT,NO
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
     *       /INT1/ C11,N11,C12,N12
     *       /INT2/ C21,N21,C22,N22
     *       /FALT/ FW,RT,VR,DV
     *       /NFA / NFAI / NS / NS2
     *       / NE /NE1,NE2,NE3
      CHARACTER FILE1*10,FILE2*10
C
C     *****  READ DATA  *****
C
      WRITE(5,4900)
 4900 FORMAT(/1X,'ENTER FILE NAME OF INPUT DATA 1 : ',$)
      READ(5,4910) FILE1
 4910 FORMAT(A10)
C
      OPEN(4,FILE=FILE1,STATUS='OLD')
C
      READ(4,5000) (NCF(I),I=1,9)
      READ(4,5010) FT,NO
      READ(4,5020) DN1,VP1,VS1,Q1
      READ(4,5020) DN2,VP2,VS2,Q2
      READ(4,5030) RT,VR,DV,FW
      READ(4,5040) NFAI
      DO 100 N=1,NFAI
      READ(4,5050) FX0(N),FZ0(N),DA(N)
      DA(N)=DA(N)*DATAN(1.D0)/45.D0
  100 CONTINUE
      READ(4,5060) NE1,NE2,NE3
      READ(4,5070) (ES1(I),I=1,NE1)
      READ(4,5080) (NES2(I),I=1,NE2)
      READ(4,5070) (ES3(I),I=1,NE3)
      READ(4,5090) THK
      READ(4,5100) C11,C12,N11,N12
      READ(4,5100) C21,C22,N21,N22
      IF(C11.LT.C12)WRITE(*,*)' C11<C12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
      IF(C21.LT.C22)WRITE(*,*)' C21<C22 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
      READ(4,5110) NS2
      READ(4,5120) (SZ2(I),SZ2(I+1),SZ2(I+2),SZ2(I+3)
     *    ,SZ2(I+4),SZ2(I+5),I=1,NS2+1,6)
C
      CLOSE(4)
C
      IF(ND1.LE.NS2*2)WRITE(6,*)' ND1 < NS2*2 !!!!!!!!!!!!!!!!!!!!!!'
      IF(ND2.LE.NE1*2)WRITE(6,*)' ND2 < NE1*2 !!!!!!!!!!!!!!!!!!!!!!'
      IF(ND2.LE.NE2*2)WRITE(6,*)' ND2 < NE2*2 !!!!!!!!!!!!!!!!!!!!!!'
      IF(ND2.LE.NE3*2)WRITE(6,*)' ND2 < NE3*2 !!!!!!!!!!!!!!!!!!!!!!'
      IF(ND3.LE.NO )WRITE(6,*)' ND3 < NO  !!!!!!!!!!!!!!!!!!!!!!!!!!'
      IF(ND4.LE.NFAI)WRITE(6,*)' ND2 < NFAI !!!!!!!!!!!!!!!!!!!!!!!!'
C
C     *****  WRITE DATA  *****
C
      WRITE(6,6000) FT,NO
      WRITE(6,6010) RT,VR,DV,FW,VP1,VP2,VS1,VS2,DN1,DN2,Q1,Q2
      WRITE(6,6020) NFAI
      DO 200 N=1,NFAI
      WRITE(6,6030) FX0(N),FZ0(N),DA(N)
  200 CONTINUE
      WRITE(6,6040) THK,C11,N11,C12,N12,C21,N21,C22,N22
      WRITE(6,6050) NS2
      WRITE(6,6060) (SZ2(I),SZ2(I+1),SZ2(I+2),SZ2(I+3)
     *   ,SZ2(I+4),SZ2(I+5),I=1,NS2+1,6)
      WRITE(6,6070) NE1,NE2,NE3
      WRITE(6,6080) (ES1(I),I=1,NE1)
      WRITE(6,6090) (NES2(I),I=1,NE2)
      WRITE(6,6100) (ES3(I),I=1,NE3)
C
C     *****  READ V2,P2  *****
C
      WRITE(5,6900)
 6900 FORMAT(/1X,'ENTER FILE NAME OF INPUT DATA 2 : ',$)
      READ(5,6910) FILE2
 6910 FORMAT(A10)
C
      OPEN(3,FILE=FILE2,STATUS='OLD')
C
      DO 130 N=1,NFAI
      DO 120 KO=1,NO
      DO 110 I=1,NS2*2,3
      READ(3,7000)V2(I,KO,N),V2(I+1,KO,N),V2(I+2,KO,N)
      READ(3,7000)P2(I,KO,N),P2(I+1,KO,N),P2(I+2,KO,N)
  110 CONTINUE
  120 CONTINUE
  130 CONTINUE
C
      CLOSE(3)
C
      IF(NCF(5).EQ.1)THEN
      WRITE(6,8000)
      WRITE(6,8100)(((V2(I,J,N),I=1,NS2*2),J=1,NO),N=1,NFAI)
      WRITE(6,8200)
      WRITE(6,8100)(((P2(I,J,N),I=1,NS2*2),J=1,NO),N=1,NFAI)
      END IF
 8000 FORMAT(//'**** FROM INPUT (NCF=5) MATRIX V2 ****'/)
 8100 FORMAT(8D12.4)
 8200 FORMAT(/'**** MATRIX P2 ****'/)
C
      WRITE(*,*)
C
      RETURN
C
 5000 FORMAT(9I5)
 5010 FORMAT(F10.0,I5)
 5020 FORMAT(4F10.0)
 5030 FORMAT(4F10.0)
 5040 FORMAT(I10)
 5050 FORMAT(3F10.0)
 5060 FORMAT(3I5)
 5070 FORMAT(6F10.0)
 5080 FORMAT(20I3)
 5090 FORMAT(F10.0)
 5100 FORMAT(2F10.0,2I5)
 5110 FORMAT(I15)
 5120 FORMAT(6F10.0)
C
 6000 FORMAT(//2X,'FT                 =',F15.7
     *         ,2X,'NUMBER OF FREQUENCY=',I8)
 6010 FORMAT(  /2X,'RISE TIME          =',F15.7
     *         ,2X,'RUPTURE VELOCITY   =',F15.7/
     *         ,2X,'DISPLACEMENT OF F. =',F15.7
     *         ,2X,'FAULT WIDTH        =',F15.7//
     *         ,2X,'VP      (1)        =',F15.7
     *         ,2X,'VP      (2)        =',F15.7/
     *         ,2X,'VS      (1)        =',F15.7
     *         ,2X,'VS      (2)        =',F15.7/
     *         ,2X,'DENSITY (1)        =',F15.7
     *         ,2X,'DENSITY (2)        =',F15.7/
     *         ,2X,'Q       (1)        =',F15.7
     *         ,2X,'Q       (2)        =',F15.7)
 6020 FORMAT(  /2X,'NUMBER OF FAULT    =',I8)
 6030 FORMAT(   2X,'FXO,FZO,DA         =',3F15.7)
 6040 FORMAT(  /2X,' THICKNESS OF LAYER =', F10.4/
     *   /2X,' FOR MEDIUM 1:'
     *   /2X,' C1(MED C), N1(  DEV NUM FROM 0 TO OM/C1  ) =',F10.4,I5,
     *   /2X,' C1(MIN C), N1(DEV NUM FROM OM/C1 TO OM/C2) =',F10.4,I5/
     *   /2X,' FOR MEDIUM 2:'
     *   /2X,' C1(MED C), N1(  DEV NUM FROM 0 TO OM/C1  ) =',F10.4,I5,
     *   /2X,' C1(MIN C), N1(DEV NUM FROM OM/C1 TO OM/C2) =',F10.4,I5
     *        //2X,' (SZ2)*6 (X=0!!!) ----'/)
 6050 FORMAT(  /2X,' NUMBER OF BOUNDARY2 ELEM.',I4/)
 6060 FORMAT(2F10.2,2X,2F10.2,2X,2F10.2)
 6070 FORMAT(  /2X,' NE1,NE2,NE3 =',3I5)
 6080 FORMAT(  /2X,' ES1 =',6F10.3)
 6090 FORMAT(   2X,' NES3 =',20I3)
 6100 FORMAT(   2X,' ES3 =',6F10.3)
C
 7000 FORMAT(1X,6D13.6)
C
      END
C     ***********************
C     *  SUBROUTINE CALVV   *
C     ***********************
C
      SUBROUTINE CALVV(SZ2,ES1,ES3,ND1,ND2,ND3,V1,V3,VF
     *                ,V2,P2,VV,PP,FX0,FZ0,DA,NCF)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 V2(ND1,ND3,*),P2(ND1,ND3,*),VV(ND2,*),PP(ND2,*)
     *          ,V1(ND2,ND3,*),V3(ND2,ND3,*),VF(ND2,*),DD
      DIMENSION SZ2(ND1),NCF(10),ES1(*),ES3(*)
     *         ,FX0(*),FZ0(*),DA(*)
      COMMON /NS/NS2/FREQ/FT,NO/NFA/NFAI/NE/NE1,NE2,NE3
      CHARACTER FILEU*10, FILEW*10
C
      WRITE(5,*)'FROM CALVV; OUTPUT OF INCIDENT WAVE'
      WRITE(5,*)'NOTE!; NUMS OF FLT AND OBS MUST BE 1 !! (OK? YES=0)'
      READ(5,*)NYES
      IF(NYES.NE.0)STOP
C
      WRITE(5,4900)
 4900 FORMAT(/1X,'ENTER OUTPUT FILE NAME FOR U: ',$)
      READ(5,4910) FILEU
      WRITE(5,4901)
 4901 FORMAT( 1X,'ENTER OUTPUT FILE NAME FOR W: ',$)
      READ(5,4910) FILEW
 4910 FORMAT(A10)
C
      OPEN(7,FILE=FILEU,STATUS='NEW')
      OPEN(8,FILE=FILEW,STATUS='NEW')
C
      DO 100 KO=1,NO
C
        OM=DFLOAT(KO)*FT
        WRITE(5,*)' NOW, OM =',OM
C
C       *** MATRIX V1 (BASIN) ***
C
        CALL CALPV1(SZ2,ND1,ES1,ND2,VV,PP,OM,NCF)
C
        DO 320 N=1,NFAI
          DO 310 J=1,NE1*2
            DD=(0.D0,0.D0)
            DO 300 I=1,NS2*2
              DD=DD-PP(J,I)*V2(I,KO,N)+VV(J,I)*P2(I,KO,N)
  300       CONTINUE
            V1(J,KO,N)=DD
  310     CONTINUE
  320   CONTINUE
C
C       *** MATRIX V3 (ROCK) ***
C
      CALL CALPV3(SZ2,ND1,ES3,ND2,VV,PP,OM,NCF)
      CALL CALVF(ES3,OM,VF,ND2,FX0,FZ0,DA,NCF)
C
      WRITE(7,1111)VF(1,1)
      WRITE(8,1111)VF(2,1)
 1111 FORMAT(1X,F8.5,1X,F8.5)
C
      DO 220 N=1,NFAI
        DO 210 J=1,NE3*2
          DD=VF(J,N)
          DO 200 I=1,NS2*2
            DD=DD-PP(J,I)*V2(I,KO,N)-VV(J,I)*P2(I,KO,N)
  200     CONTINUE
          V3(J,KO,N)=DD
  210   CONTINUE
  220 CONTINUE
C
  100 CONTINUE
C
      CLOSE(7)
      CLOSE(8)
      WRITE(5,*)'USE MINV2!!!'
      STOP
C
C     *****  CHECKER'S FLAG NO.9  *****
C
      IF(NCF(9).EQ.1) THEN
      WRITE(6,*) '*****  FROM CALVV  (NCF=9) MATRIX V1  *****'
      WRITE(6,*) (((V1(I,J,N),I=1,NO),J=1,NS2*2),N=1,NFAI)
      WRITE(6,*) '*****  VECTOR V3  *****'
      WRITE(6,*) (((V3(I,J,N),I=1,NO),J=1,NS2*2),N=1,NFAI)
      END IF
C
      RETURN
      END
C     ***********************
C     *  SUBROUTINE CALPV1  *
C     ***********************
C
      SUBROUTINE CALPV1(SZ2,ND1,ES1,ND2,V,P,OM,NCF)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SZ2(*),NCF(10),ES1(*)
      COMPLEX*16 V(ND2,ND1),P(ND2,ND1)
     *          ,SQ1,SQ2,SP1,SP2,WP1,WP2,WS1,WS2,SM,RAMD
     *          ,UXX1,UZZ1,UZX1,UXZ1,UXX2,UZZ2,UZX2,UXZ2
     *          ,PXX1,PZZ1,PZX1,PXZ1,PXX2,PZZ2,PZX2,PXZ2
      COMMON /FDAT/ WP1,WP2,WS1,WS2,OM2,XI,ZJ1,ZJ2,SM,RAMD
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
     *       /INT2/ C1,N1,C2,N2
     *       / NS / NS2  / NE / NE1,NE2,NE3
      EXTERNAL FB21,FB22
C
      OM2=OM*OM
C
      SP1=VP1*DCMPLX(1.D0,0.5D0/Q1)
      SP2=VP2*DCMPLX(1.D0,0.5D0/Q2)
      WP1=OM2/(SP1*SP1)
      WP2=OM2/(SP2*SP2)
C
      SQ1=VS1*DCMPLX(1.D0,0.5D0/Q1)
      SQ2=VS2*DCMPLX(1.D0,0.5D0/Q2)
      WS1=OM2/(SQ1*SQ1)
      WS2=OM2/(SQ2*SQ2)
C
      WNB=OM/C1
      WNA=WNB/10000.D0
      WNC=OM/C2
C
      DO 110 I=1,NE1
C
        XI=ES1(I)
C
        DO 100 J=1,NS2
C
          ZJ1=SZ2(J)
          ZJ2=SZ2(J+1)
C
          IF(ZJ2.LE.THK)THEN
            SM=DN1*SQ1*SQ1
            RAMD=DN1*SP1*SP1-2.D0*SM
            CALL NEWTNB(FB21,WNA,WNB,N1,UXX1,UZZ1,UZX1,UXZ1
     *                 ,PXX1,PZZ1,PZX1,PXZ1,IER)
            IF(IER.NE.0)WRITE(*,*)' FB21 1: IER=',IER
            CALL NEWTNB(FB21,WNB,WNC,N2,UXX2,UZZ2,UZX2,UXZ2
     *                 ,PXX2,PZZ2,PZX2,PXZ2,IER)
            IF(IER.NE.0)WRITE(*,*)' FB21 2: IER=',IER
          ELSE
            SM=DN2*SQ2*SQ2
            RAMD=DN2*SP2*SP2-2.D0*SM
            CALL NEWTNB(FB22,WNA,WNB,N1,UXX1,UZZ1,UZX1,UXZ1
     *                 ,PXX1,PZZ1,PZX1,PXZ1,IER)
            IF(IER.NE.0)WRITE(*,*)' FB22 1: IER=',IER
            CALL NEWTNB(FB22,WNB,WNC,N2,UXX2,UZZ2,UZX2,UXZ2
     *                 ,PXX2,PZZ2,PZX2,PXZ2,IER)
            IF(IER.NE.0)WRITE(*,*)' FB22 2: IER=',IER
          END IF
C
          V(2*I-1,2*J-1)=UXX1+UXX2
          V(2*I-1,2*J  )=UZX1+UZX2
          V(2*I  ,2*J-1)=UXZ1+UXZ2
          V(2*I  ,2*J  )=UZZ1+UZZ2
C
          P(2*I-1,2*J-1)=PXX1+PXX2
          P(2*I-1,2*J  )=PZX1+PZX2
          P(2*I  ,2*J-1)=PXZ1+PXZ2
          P(2*I  ,2*J  )=PZZ1+PZZ2
C
  100   CONTINUE
C
  110 CONTINUE
C
C *** CHECK ***
C
      IF(NCF(6).EQ.1)THEN
      WRITE(*,*)' *** FROM CALPV1 (NCF=6): V ***'
      WRITE(*,*)((V(I,J),J=1,2*NS2),I=1,2*NE1)
      WRITE(*,*)' *** FROM CALPV1 (NCF=6): P ***'
      WRITE(*,*)((P(I,J),J=1,2*NS2),I=1,2*NE1)
      END IF
C
      RETURN
      END
      SUBROUTINE NEWTNB(SUB,A,B,N,UXX,UZZ,UZX,UXZ
     *                 ,PXX,PZZ,PZX,PXZ,IER)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 UYXX,UYZZ,UYXZ,UYZX,PYXX,PYZX,PYXZ,PYZZ
     *          ,UFXX,UFZZ,UFXZ,UFZX,PFXX,PFZX,PFXZ,PFZZ
     *          ,UQXX,UQZZ,UQXZ,UQZX,PQXX,PQZX,PQXZ,PQZZ
     *          ,UXX,UZZ,UXZ,UZX,PXX,PZX,PXZ,PZZ
C
      IF(N.LT.1)THEN
        IER=1000
        UXX=(0.D0,0.D0)
        UZX=(0.D0,0.D0)
        UXZ=(0.D0,0.D0)
        UZZ=(0.D0,0.D0)
        PXX=(0.D0,0.D0)
        PZX=(0.D0,0.D0)
        PXZ=(0.D0,0.D0)
        PZZ=(0.D0,0.D0)
        GO TO 9999
      END IF
      IF(A.EQ.B)THEN
        IER=2000
        UXX=(0.D0,0.D0)
        UZX=(0.D0,0.D0)
        UXZ=(0.D0,0.D0)
        UZZ=(0.D0,0.D0)
        PXX=(0.D0,0.D0)
        PZX=(0.D0,0.D0)
        PXZ=(0.D0,0.D0)
        PZZ=(0.D0,0.D0)
        GO TO 9999
      END IF
C
      H=(B-A)/N
      CALL SUB(A,UYXX,UYZX,UYXZ,UYZZ,PYXX,PYZX,PYXZ,PYZZ)
      UQXX=(0.D0,0.D0)
      UQZX=(0.D0,0.D0)
      UQXZ=(0.D0,0.D0)
      UQZZ=(0.D0,0.D0)
      PQXX=(0.D0,0.D0)
      PQZX=(0.D0,0.D0)
      PQXZ=(0.D0,0.D0)
      PQZZ=(0.D0,0.D0)
C
      DO  100 I=1,N
      X=A+H*I
      CALL SUB(X,UFXX,UFZX,UFXZ,UFZZ,PFXX,PFZX,PFXZ,PFZZ)
      UQXX = UQXX +H*(UYXX +UFXX )/2.D0
      UQZX = UQZX +H*(UYZX +UFZX )/2.D0
      UQXZ = UQXZ +H*(UYXZ +UFXZ )/2.D0
      UQZZ = UQZZ +H*(UYZZ +UFZZ )/2.D0
      PQXX = PQXX +H*(PYXX +PFXX )/2.D0
      PQZX = PQZX +H*(PYZX +PFZX )/2.D0
      PQXZ = PQXZ +H*(PYXZ +PFXZ )/2.D0
      PQZZ = PQZZ +H*(PYZZ +PFZZ )/2.D0
      UYXX = UFXX
      UYZX = UFZX
      UYXZ = UFXZ
      UYZZ = UFZZ
      PYXX = PFXX
      PYZX = PFZX
      PYXZ = PFXZ
      PYZZ = PFZZ
  100 CONTINUE
C
      UXX = UQXX
      UZX = UQZX
      UXZ = UQXZ
      UZZ = UQZZ
      PXX = PQXX
      PZX = PQZX
      PXZ = PQXZ
      PZZ = PQZZ
      IER=0
C
 9999 RETURN
C
      END
C
C *** THE CASE FOR DOMAIN 2 OF THAT SP AND OP ARE IN THE SURFACE LAYER,
C     AND SP IS ON THE FREE SURFACE
C
      SUBROUTINE FB21(WN,UXX,UZX,UXZ,UZZ,PXX,PZX,PXZ,PZZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 UXX,UZZ,UZX,UXZ,PXX,PZZ,PZX,PXZ
     *          ,SM,RAMD,WP1,WP2,WS1,WS2
     *          ,CPX,CMP,CMS,CMPN,CMSN,CCC
     *          ,AA(4,2),DAR2(4,4)
     *          ,RAM,PM,CP,SP,STP,SDP,RBM,QM,CQ,SQ,STQ,SDQ
     *          ,GM,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
     *          ,GN,RAN,RBN,B11,B12,B13,B14,B21,B22,B23,B24
     *          ,RL,RK,RN,RM,FR,GR2,RR2,HR2,SR2
     *          ,FWZ0,FWX0,FUZ0,FUX0
     *          ,UXXX,UZXX,UXZX,UZZX,UXXZ,UZXZ,UXZZ,UZZZ
     *          ,DFWZ,DFWX,DFUZ,DFUX,DFWZZ,DFWXZ,DFUZZ,DFUXZ
      COMMON /FDAT/ WP1,WP2,WS1,WS2,OM2,XI,ODP1,ODP2,SM,RAMD
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
C
      WN2=WN*WN
      CV2=OM2/WN2
      CPX=DCMPLX(0.D0,1.D0)
C
      CMP=CDSQRT(WP1-WN2)
      CMS=CDSQRT(WS1-WN2)
      IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
      IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
      RAM=CMP/WN
      PM=CMP*THK
      CP=CDCOS(PM)
      SP=CDSIN(PM)
      STP=RAM*SP
      SDP=SP/RAM
C
      RBM=CMS/WN
      QM=CMS*THK
      CQ=CDCOS(QM)
      SQ=CDSIN(QM)
      STQ=RBM*SQ
      SDQ=SQ/RBM
C
      GM=2.D0*WN2/WS1
      C11=GM
      C12=GM-1.D0
      C13=GM*GM
      C14=C12*C12
      C16=DN1*CV2
      C15=1.D0/C16
C
      D11=CP-CQ
      D14=C14*SDP+C13*STQ
      D15=C13*STP+C14*SDQ
C
      AA(1,1)= C11*CP-C12*CQ
      AA(1,2)= (C12*SDP+C11*STQ)*CPX
      AA(2,1)=-(C11*STP+C12*SDQ)*CPX
      AA(2,2)=-C12*CP+C11*CQ
      AA(3,1)= C16*C11*C12*D11
      AA(3,2)= C16*D14*CPX
      AA(4,1)= C16*D15*CPX
      AA(4,2)= AA(3,1)
C
      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
      RAN=CMPN/WN
      RBN=CMSN/WN
C
      GN=2.D0*WN2/WS2
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DN2*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 *** OP,Z (IN THE SURFACE LAYER) ***
C
      CP=CDCOS(CMP*(THK-ODP2))-CDCOS(CMP*(THK-ODP1))
      SP=CDSIN(CMP*(THK-ODP2))-CDSIN(CMP*(THK-ODP1))
      STP=RAM*SP
      SDP=SP/RAM
C
      CQ=CDCOS(CMS*(THK-ODP2))-CDCOS(CMS*(THK-ODP1))
      SQ=CDSIN(CMS*(THK-ODP2))-CDSIN(CMS*(THK-ODP1))
      STQ=RBM*SQ
      SDQ=SQ/RBM
C
      D11=CP-CQ
      D12=SDP+STQ
      D13=STP+SDQ
C
      DAR2(1,3)=-C15*D11
      DAR2(1,4)= C15*D12*CPX
      DAR2(2,3)= C15*D13*CPX
      DAR2(2,4)= DAR2(1,3)
      DAR2(3,3)=-C12*CP+C11*CQ
      DAR2(3,4)= (C12*SDP+C11*STQ)*CPX
      DAR2(4,3)=-(C11*STP+C12*SDQ)*CPX
      DAR2(4,4)= C11*CP-C12*CQ
C
      GR2=B11*DAR2(1,3)+B12*DAR2(2,3)+B13*DAR2(3,3)+B14*DAR2(4,3)
      RR2=B11*DAR2(1,4)+B12*DAR2(2,4)+B13*DAR2(3,4)+B14*DAR2(4,4)
      HR2=B21*DAR2(1,3)+B22*DAR2(2,3)+B23*DAR2(3,3)+B24*DAR2(4,3)
      SR2=B21*DAR2(1,4)+B22*DAR2(2,4)+B23*DAR2(3,4)+B24*DAR2(4,4)
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      DFWZZ= FWZ0
      DFWXZ= FUZ0
      DFUZZ= FWX0
      DFUXZ= FUX0
C
C *** OP (IN THE SURFACE LAYER) ***
C
      PM=-SP/CMP
      SP= CP/CMP
      CP=PM
      STP=RAM*SP
      SDP=SP/RAM
C
      QM=-SQ/CMS
      SQ= CQ/CMS
      CQ=QM
      STQ=RBM*SQ
      SDQ=SQ/RBM
C
      D11=CP-CQ
      D12=SDP+STQ
      D13=STP+SDQ
C
      DAR2(1,3)=-C15*D11
      DAR2(1,4)= C15*D12*CPX
      DAR2(2,3)= C15*D13*CPX
      DAR2(2,4)= DAR2(1,3)
      DAR2(3,3)=-C12*CP+C11*CQ
      DAR2(3,4)= (C12*SDP+C11*STQ)*CPX
      DAR2(4,3)=-(C11*STP+C12*SDQ)*CPX
      DAR2(4,4)= C11*CP-C12*CQ
C
      GR2=B11*DAR2(1,3)+B12*DAR2(2,3)+B13*DAR2(3,3)+B14*DAR2(4,3)
      RR2=B11*DAR2(1,4)+B12*DAR2(2,4)+B13*DAR2(3,4)+B14*DAR2(4,4)
      HR2=B21*DAR2(1,3)+B22*DAR2(2,3)+B23*DAR2(3,3)+B24*DAR2(4,3)
      SR2=B21*DAR2(1,4)+B22*DAR2(2,4)+B23*DAR2(3,4)+B24*DAR2(4,4)
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      DFWZ= FWZ0
      DFWX= FUZ0
      DFUZ= FWX0
      DFUX= FUX0
C
C *** UXX - PZZ ***
C
      PAI=3.1415927D0
      SF=DSIN(WN*XI)
      CF=DCOS(WN*XI)
C
      CCC=1.D0/(PAI*FR)
      UZXX=-CCC*DFWX*CF
      UXZX= CCC*DFUZ*CF
      UXXX= CPX*CCC*DFUX*SF
      UZZX=-CPX*CCC*DFWZ*SF
C
      CCC=CCC/WN
      UZX = CCC*DFWX*SF
      UXZ =-CCC*DFUZ*SF
      UZXZ= CCC*DFWXZ*SF
      UXZZ=-CCC*DFUZZ*SF
C
      CCC=CCC*CPX
      UXX = CCC*DFUX*CF
      UZZ =-CCC*DFWZ*CF
      UXXZ= CCC*DFUXZ*CF
      UZZZ=-CCC*DFWZZ*CF
C
      PXX=-RAMD*(UXXX+UZXZ)-2.D0*SM*UXXX
      PZX=-SM*(UXXZ+UZXX)
      PXZ=-RAMD*(UXZX+UZZZ)-2.D0*SM*UXZX
      PZZ=-SM*(UXZZ+UZZX)
C
      RETURN
      END
C
C *** THE CASE FOR DOMAIN 2 OF THAT SP IS ON THE FREE SURFACE,
C     AND OP IS IN THE BOTTOM BEDROCK
C
      SUBROUTINE FB22(WN,UXX,UZX,UXZ,UZZ,PXX,PZX,PXZ,PZZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 UXX,UZZ,UZX,UXZ,PXX,PZZ,PZX,PXZ
     *          ,SM,RAMD,WP1,WP2,WS1,WS2
     *          ,CPX,CMP,CMS,CMPN,CMSN,CCC
     *          ,AA(4,2)
     *          ,RAM,PM,CP,SP,STP,SDP,RBM,QM,CQ,SQ,STQ,SDQ
     *          ,GM,C11,C12,C13,C14,C16,D11,D14,D15
     *          ,GN,RAN,RBN,B11,B12,B13,B14,B21,B22,B23,B24
     *          ,RL,RK,RN,RM,FR,GR2,RR2,HR2,SR2
     *          ,FWZ0,FWX0,FUZ0,FUX0
     *          ,UXXX,UZXX,UXZX,UZZX,UXXZ,UZXZ,UXZZ,UZZZ
     *          ,DFWZ,DFWX,DFUZ,DFUX,DFWZZ,DFWXZ,DFUZZ,DFUXZ
      COMMON /FDAT/ WP1,WP2,WS1,WS2,OM2,XI,ODP1,ODP2,SM,RAMD
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
C
      WN2=WN*WN
      CV2=OM2/WN2
      CPX=DCMPLX(0.D0,1.D0)
C
      CMP=CDSQRT(WP1-WN2)
      CMS=CDSQRT(WS1-WN2)
      IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
      IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
      RAM=CMP/WN
      PM=CMP*THK
      CP=CDCOS(PM)
      SP=CDSIN(PM)
      STP=RAM*SP
      SDP=SP/RAM
C
      RBM=CMS/WN
      QM=CMS*THK
      CQ=CDCOS(QM)
      SQ=CDSIN(QM)
      STQ=RBM*SQ
      SDQ=SQ/RBM
C
      GM=2.D0*WN2/WS1
      C11=GM
      C12=GM-1.D0
      C13=GM*GM
      C14=C12*C12
      C16=DN1*CV2
C
      D11=CP-CQ
      D14=C14*SDP+C13*STQ
      D15=C13*STP+C14*SDQ
C
      AA(1,1)= C11*CP-C12*CQ
      AA(1,2)= (C12*SDP+C11*STQ)*CPX
      AA(2,1)=-(C11*STP+C12*SDQ)*CPX
      AA(2,2)=-C12*CP+C11*CQ
      AA(3,1)= C16*C11*C12*D11
      AA(3,2)= C16*D14*CPX
      AA(4,1)= C16*D15*CPX
      AA(4,2)= AA(3,1)
C
      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
      RAN=CMPN/WN
      RBN=CMSN/WN
C
      GN=2.D0*WN2/WS2
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DN2*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 *** OP,Z (IN THE BOTTOM BEDROCK)  ***
C
      CP=CDEXP(-CPX*CMPN*(ODP2-THK))-CDEXP(-CPX*CMPN*(ODP1-THK))
      CQ=CDEXP(-CPX*CMSN*(ODP2-THK))-CDEXP(-CPX*CMSN*(ODP1-THK))
C
      RR2= B14*CP
      GR2=-RAN*RR2
      HR2= B14*CQ
      SR2= RBN*HR2
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      DFWZZ= FWZ0
      DFWXZ= FUZ0
      DFUZZ= FWX0
      DFUXZ= FUX0
C
C *** OP (IN THE BOTTOM BEDROCK)  ***
C
      CP=CPX*CP/CMPN
      CQ=CPX*CQ/CMSN
C
      RR2= B14*CP
      GR2=-RAN*RR2
      HR2= B14*CQ
      SR2= RBN*HR2
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      DFWZ= FWZ0
      DFWX= FUZ0
      DFUZ= FWX0
      DFUX= FUX0
C
C *** UXX - PXZ ***
C
      PAI=3.1415927D0
      SF=DSIN(WN*XI)
      CF=DCOS(WN*XI)
C
      CCC=1.D0/(PAI*FR)
      UZXX=-CCC*DFWX*CF
      UXZX= CCC*DFUZ*CF
      UXXX= CPX*CCC*DFUX*SF
      UZZX=-CPX*CCC*DFWZ*SF
C
      CCC=CCC/WN
      UZX = CCC*DFWX*SF
      UXZ =-CCC*DFUZ*SF
      UZXZ= CCC*DFWXZ*SF
      UXZZ=-CCC*DFUZZ*SF
C
      CCC=CCC*CPX
      UXX = CCC*DFUX*CF
      UZZ =-CCC*DFWZ*CF
      UXXZ= CCC*DFUXZ*CF
      UZZZ=-CCC*DFWZZ*CF
C
      PXX=-RAMD*(UXXX+UZXZ)-2.D0*SM*UXXX
      PZX=-SM*(UXXZ+UZXX)
      PXZ=-RAMD*(UXZX+UZZZ)-2.D0*SM*UXZX
      PZZ=-SM*(UXZZ+UZZX)
C
      RETURN
      END
C     ***********************
C     *  SUBROUTINE CALPV3  *
C     ***********************
C
      SUBROUTINE CALPV3(SZ2,ND1,ES3,ND2,V,P,OM,NCF)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SZ2(*),NCF(10),ES3(*)
      COMPLEX*16 V(ND2,ND1),P(ND2,ND1)
     *          ,SQ2,SP2,WP2,WS2,SM,RAMD
     *          ,UXX1,UZZ1,UZX1,UXZ1,UXX2,UZZ2,UZX2,UXZ2
     *          ,PXX1,PZZ1,PZX1,PXZ1,PXX2,PZZ2,PZX2,PXZ2
      COMMON /FDT1/ WP2,WS2,OM2,XI,ZJ1,ZJ2,SM,RAMD
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
     *       /INT1/ C1,N1,C2,N2
     *       / NS / NS2  / NE / NE1,NE2,NE3
      EXTERNAL FB1
C
      OM2=OM*OM
C
      SP2=VP2*DCMPLX(1.D0,0.5D0/Q2)
      WP2=OM2/(SP2*SP2)
      SQ2=VS2*DCMPLX(1.D0,0.5D0/Q2)
      WS2=OM2/(SQ2*SQ2)
      SM=DN2*SQ2*SQ2
      RAMD=DN2*SP2*SP2-2.D0*SM
C
      WNB=OM/C1
      WNA=WNB/10000.D0
      WNC=OM/C2
C
      DO 110 I=1,NE3
C
        XI=ES3(I)
C
        DO 100 J=1,NS2
C
          ZJ1=SZ2(J)
          ZJ2=SZ2(J+1)
C
          CALL NEWTNB(FB1,WNA,WNB,N1,UXX1,UZZ1,UZX1,UXZ1
     *               ,PXX1,PZZ1,PZX1,PXZ1,IER)
          IF(IER.NE.0)WRITE(*,*)' FB1 1: IER=',IER
          CALL NEWTNB(FB1,WNB,WNC,N2,UXX2,UZZ2,UZX2,UXZ2
     *               ,PXX2,PZZ2,PZX2,PXZ2,IER)
          IF(IER.NE.0)WRITE(*,*)' FB1 2: IER=',IER
C
          V(2*I-1,2*J-1)=UXX1+UXX2
          V(2*I-1,2*J  )=UZX1+UZX2
          V(2*I  ,2*J-1)=UXZ1+UXZ2
          V(2*I  ,2*J  )=UZZ1+UZZ2
C
          P(2*I-1,2*J-1)=PXX1+PXX2
          P(2*I-1,2*J  )=PZX1+PZX2
          P(2*I  ,2*J-1)=PXZ1+PXZ2
          P(2*I  ,2*J  )=PZZ1+PZZ2
C
  100   CONTINUE
C
  110 CONTINUE
C
C *** CHECK ***
C
      IF(NCF(7).EQ.1)THEN
      WRITE(*,*)' *** FROM CALPV3 (NCF=7): V ***'
      WRITE(*,*)((V(I,J),J=1,2*NS2),I=1,2*NE3)
      WRITE(*,*)' *** FROM CALPV3 (NCF=7): P ***'
      WRITE(*,*)((P(I,J),J=1,2*NS2),I=1,2*NE3)
      END IF
C
      RETURN
      END
C
C *** THE CASE FOR DOMAIN 1 OF THAT SP IS ON THE FREE SURFACE
C
      SUBROUTINE FB1(WN,UXX,UZX,UXZ,UZZ,PXX,PZX,PXZ,PZZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 UXX,UZZ,UZX,UXZ,PXX,PZZ,PZX,PXZ
     *          ,SM,RAMD,WP2,WS2
     *          ,CPX,CMPN,CMSN,CCC
     *          ,CP,CQ
     *          ,GN,RAN,RBN,B11,B12,B14,B21,B22
     *          ,RL,RK,RN,RM,FR,GR2,RR2,HR2,SR2
     *          ,FWZ0,FWX0,FUZ0,FUX0
     *          ,UXXX,UZXX,UXZX,UZZX,UXXZ,UZXZ,UXZZ,UZZZ
     *          ,DFWZ,DFWX,DFUZ,DFUX,DFWZZ,DFWXZ,DFUZZ,DFUXZ
      COMMON /FDT1/ WP2,WS2,OM2,XI,ODP1,ODP2,SM,RAMD
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
C
      WN2=WN*WN
      CV2=OM2/WN2
      CPX=DCMPLX(0.D0,1.D0)
C
      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
      RAN=CMPN/WN
      RBN=CMSN/WN
C
      GN=2.D0*WN2/WS2
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DN2*CV2)
      B21=-B12
      B22=GN*RBN
C
      RL=B11
      RK=B12
      RN=B21
      RM=B22
C
      FR=RN*RK-RL*RM
C
C *** OP,Z ***
C
      CP=CDEXP(-CPX*CMPN*ODP2)-CDEXP(-CPX*CMPN*ODP1)
      CQ=CDEXP(-CPX*CMSN*ODP2)-CDEXP(-CPX*CMSN*ODP1)
C
      RR2= B14*CP
      GR2=-RAN*RR2
      HR2= B14*CQ
      SR2= RBN*HR2
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      DFWZZ= FWZ0
      DFWXZ= FUZ0
      DFUZZ= FWX0
      DFUXZ= FUX0
C
C *** OP ***
C
      CP=CPX*CP/CMPN
      CQ=CPX*CQ/CMSN
C
      RR2= B14*CP
      GR2=-RAN*RR2
      HR2= B14*CQ
      SR2= RBN*HR2
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      DFWZ= FWZ0
      DFWX= FUZ0
      DFUZ= FWX0
      DFUX= FUX0
C
C *** UXX - PXZ ***
C
      PAI=3.1415927D0
      SF=DSIN(WN*XI)
      CF=DCOS(WN*XI)
C
      CCC=1.D0/(PAI*FR)
      UZXX=-CCC*DFWX*CF
      UXZX= CCC*DFUZ*CF
      UXXX= CPX*CCC*DFUX*SF
      UZZX=-CPX*CCC*DFWZ*SF
C
      CCC=CCC/WN
      UZX = CCC*DFWX*SF
      UXZ =-CCC*DFUZ*SF
      UZXZ= CCC*DFWXZ*SF
      UXZZ=-CCC*DFUZZ*SF
C
      CCC=CCC*CPX
      UXX = CCC*DFUX*CF
      UZZ =-CCC*DFWZ*CF
      UXXZ= CCC*DFUXZ*CF
      UZZZ=-CCC*DFWZZ*CF
C
      PXX=RAMD*(UXXX+UZXZ)+2.D0*SM*UXXX
      PZX=SM*(UXXZ+UZXX)
      PXZ=RAMD*(UXZX+UZZZ)+2.D0*SM*UXZX
      PZZ=SM*(UXZZ+UZZX)
C
      RETURN
      END
C     **********************
C     *  SUBROUTINE CALVF  *
C     **********************
C
      SUBROUTINE CALVF(ES3,OM,F,ND2,FX0,FZ0,DA,NCF)
C
      IMPLICIT REAL*8(A-H,M,O-Z)
      COMPLEX*16 F(ND2,*),CPX,FSS,FWW,AA
     *          ,SQ2,SP2,WP2,WS2,SM
     *          ,UXXX1,UZXX1,UXZX1,UZZX1,UXXZ1,UZXZ1,UXZZ1,UZZZ1
     *          ,UXXX2,UZXX2,UXZX2,UZZX2,UXXZ2,UZXZ2,UXZZ2,UZZZ2
      DIMENSION ES3(*),NCF(10),FX0(*),FZ0(*),DA(*)
      COMMON /FDF1/ WP2,WS2,OM2,XF,ZF,XI
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
     *       /INT1/ C1,N1,C2,N2
     *       /FALT/ FW,RT,VR,DV
     *       /NFA / NFAI /NE/NE1,NE2,NE3
      EXTERNAL FFB
C
      OM2=OM*OM
      CPX=(0.D0,1.D0)
C
      SP2=VP2*DCMPLX(1.D0,0.5D0/Q2)
      WP2=OM2/(SP2*SP2)
      SQ2=VS2*DCMPLX(1.D0,0.5D0/Q2)
      WS2=OM2/(SQ2*SQ2)
      SM=DN2*SQ2*SQ2
C
      TW=FW/VR
      XW=OM*TW/2.D0
      XT=OM*RT/2.D0
      FWW=(DSIN(XW)/XW)*CDEXP(-CPX*XW)
      FSS=(DSIN(XT)/(XT*CPX*OM))*CDEXP(-CPX*XT)
      AA=SM*DV*FW*FWW*FSS
C
      WNB=OM/C1
      WNA=WNB/10000.D0
      WNC=OM/C2
C
      DO 110 N=1,NFAI
C
        XF=FX0(N)
        ZF=FZ0(N)
C
        SF=DSIN(2.D0*DA(N))
        CF=DCOS(2.D0*DA(N))
C
        DO 100 I=1,NE3
C
          XI=ES3(I)
C
          CALL NEWTNF(FFB,WNA,WNB,N1,IER
     *           ,UXXX1,UZXX1,UXZX1,UXXZ1,UZZX1,UZXZ1,UXZZ1,UZZZ1)
          IF(IER.NE.0)WRITE(*,*)' FFB 1: IER=',IER
          CALL NEWTNF(FFB,WNB,WNC,N2,IER
     *           ,UXXX2,UZXX2,UXZX2,UXXZ2,UZZX2,UZXZ2,UXZZ2,UZZZ2)
          IF(IER.NE.0)WRITE(*,*)' FFB 2: IER=',IER
C
          F(2*I-1,N)=AA*(SF*(UXXX1+UXXX2-UZXZ1-UZXZ2)
     *                  -CF*(UXXZ1+UXXZ2+UZXX1+UZXX2))
          F(2*I  ,N)=AA*(SF*(UXZX1+UXZX2-UZZZ1-UZZZ2)
     *                  -CF*(UXZZ1+UXZZ2+UZZX1+UZZX2))
C
  100   CONTINUE
C
  110 CONTINUE
C
C     *****  CHECKER'S FLAG NO.2  *****
C
      IF(NCF(8).EQ.1) THEN
      WRITE(6,*) ' *****  FROM CALVF (NCF=8) VECTOR VF  *****'
      WRITE(6,*) ((F(I,N),I=1,2*NE3),N=1,NFAI)
      WRITE(6,*) ' ***** VECTOR ABS(VF)  *****'
      WRITE(6,*) ((CDABS(F(I,N)),I=1,2*NE3),N=1,NFAI)
      END IF
C
      RETURN
      END
      SUBROUTINE NEWTNF (SUB,A,B,N,IER
     *             ,UXXX,UZXX,UXZX,UXXZ,UZZX,UZXZ,UXZZ,UZZZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 YXXX,YZXX,YXZX,YXXZ,YZZX,YZXZ,YXZZ,YZZZ
     *          ,FXXX,FZXX,FXZX,FXXZ,FZZX,FZXZ,FXZZ,FZZZ
     *          ,QXXX,QZXX,QXZX,QXXZ,QZZX,QZXZ,QXZZ,QZZZ
     *          ,UXXX,UZXX,UXZX,UXXZ,UZZX,UZXZ,UXZZ,UZZZ
C
      IF(N.LT.1)THEN
        IER=1000
        UXXX=(0.D0,0.D0)
        UZXX=(0.D0,0.D0)
        UXZX=(0.D0,0.D0)
        UXXZ=(0.D0,0.D0)
        UZZX=(0.D0,0.D0)
        UZXZ=(0.D0,0.D0)
        UXZZ=(0.D0,0.D0)
        UZZZ=(0.D0,0.D0)
        GO TO 9999
      END IF
      IF(A.EQ.B)THEN
        IER=2000
        UXXX=(0.D0,0.D0)
        UZXX=(0.D0,0.D0)
        UXZX=(0.D0,0.D0)
        UXXZ=(0.D0,0.D0)
        UZZX=(0.D0,0.D0)
        UZXZ=(0.D0,0.D0)
        UXZZ=(0.D0,0.D0)
        UZZZ=(0.D0,0.D0)
        GO TO 9999
      END IF
C
      H=(B-A)/N
      CALL SUB(A,YXXX,YZXX,YXZX,YXXZ,YZZX,YZXZ,YXZZ,YZZZ)
      QXXX=(0.D0,0.D0)
      QZXX=(0.D0,0.D0)
      QXZX=(0.D0,0.D0)
      QXXZ=(0.D0,0.D0)
      QZZX=(0.D0,0.D0)
      QZXZ=(0.D0,0.D0)
      QXZZ=(0.D0,0.D0)
      QZZZ=(0.D0,0.D0)
C
      DO  100 I=1,N
      X=A+H*I
      CALL SUB(X,FXXX,FZXX,FXZX,FXXZ,FZZX,FZXZ,FXZZ,FZZZ)
      QXXX = QXXX+H*(YXXX+FXXX)/2.D0
      QZXX = QZXX+H*(YZXX+FZXX)/2.D0
      QXZX = QXZX+H*(YXZX+FXZX)/2.D0
      QXXZ = QXXZ+H*(YXXZ+FXXZ)/2.D0
      QZZX = QZZX+H*(YZZX+FZZX)/2.D0
      QZXZ = QZXZ+H*(YZXZ+FZXZ)/2.D0
      QXZZ = QXZZ+H*(YXZZ+FXZZ)/2.D0
      QZZZ = QZZZ+H*(YZZZ+FZZZ)/2.D0
      YXXX = FXXX
      YZXX = FZXX
      YXZX = FXZX
      YXXZ = FXXZ
      YZZX = FZZX
      YZXZ = FZXZ
      YXZZ = FXZZ
      YZZZ = FZZZ
  100 CONTINUE
C
      UXXX = QXXX
      UZXX = QZXX
      UXZX = QXZX
      UXXZ = QXXZ
      UZZX = QZZX
      UZXZ = QZXZ
      UXZZ = QXZZ
      UZZZ = QZZZ
      IER=0
C
 9999 RETURN
C
      END
C
C *** THE CASE FOR DOMAIN 1 OF THAT SP IS ON THE FREE-SURFACE ***
C
      SUBROUTINE FFB(WN,UXXX,UZXX,UXZX,UXXZ,UZZX,UZXZ,UXZZ,UZZZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 WP2,WS2
     *          ,CPX,CMPN,CMSN,CCC
     *          ,CP,CQ
     *          ,GN,RAN,RBN,B11,B12,B14,B21,B22
     *          ,RL,RK,RN,RM,FR,GR2,RR2,HR2,SR2
     *          ,FWZ0,FWX0,FUZ0,FUX0
     *          ,FWZ,FWX,FUZ,FUX,FWZZ,FWXZ,FUZZ,FUXZ
     *          ,UXXX,UZXX,UXZX,UXXZ,UZZX,UZXZ,UXZZ,UZZZ
      COMMON /FDF1/ WP2,WS2,OM2,XF,ZF,XI
     *       /CONS/ THK,DN1,DN2,VP1,VP2,VS1,VS2,Q1,Q2
C
      WN2=WN*WN
      CV2=OM2/WN2
      CPX=DCMPLX(0.D0,1.D0)
C
      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
      RAN=CMPN/WN
      RBN=CMSN/WN
C
      GN=2.D0*WN2/WS2
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DN2*CV2)
      B21=-B12
      B22=GN*RBN
C
      RL=B11
      RK=B12
      RN=B21
      RM=B22
C
      FR=RN*RK-RL*RM
C
C *** OP ***
C
      CP=CDEXP(-CPX*CMPN*ZF)
      CQ=CDEXP(-CPX*CMSN*ZF)
C
      RR2= B14*CP
      GR2=-RAN*RR2
      HR2= B14*CQ
      SR2= RBN*HR2
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      FWZ= FWZ0
      FWX= FUZ0
      FUZ= FWX0
      FUX= FUX0
C
C *** OP,Z ***
C
      CP=-CPX*CP*CMPN
      CQ=-CPX*CQ*CMSN
C
      RR2= B14*CP
      GR2=-RAN*RR2
      HR2= B14*CQ
      SR2= RBN*HR2
C
      FWZ0=GR2*RN-HR2*RL
      FWX0=RR2*RN-SR2*RL
      FUZ0=GR2*RM-HR2*RK
      FUX0=RR2*RM-SR2*RK
C
      FWZZ= FWZ0
      FWXZ= FUZ0
      FUZZ= FWX0
      FUXZ= FUX0
C
C *** UXXX - UZZZ ***
C
      PAI=3.1415927D0
      SF=DSIN(WN*(XF-XI))
      CF=DCOS(WN*(XF-XI))
C
      CCC=1.D0/(PAI*FR)
      UZXX=-CCC*FWX*CF
      UXZX= CCC*FUZ*CF
      UXXX=-CPX*CCC*FUX*SF
      UZZX= CPX*CCC*FWZ*SF
C
      CCC=CCC/WN
      UZXZ=-CCC*FWXZ*SF
      UXZZ= CCC*FUZZ*SF
      IF(ZF.EQ.0.D0)THEN
        UXXZ=( CPX*CCC*FUXZ +WS2/(PAI*DN2*OM2))*CF
        UZZZ=(-CPX*CCC*FWZZ +WP2/(PAI*DN2*OM2))*CF
      ELSE
        UXXZ= CPX*CCC*FUXZ*CF
        UZZZ=-CPX*CCC*FWZZ*CF
      END IF
C
      RETURN
      END
C     ***********************
C     *  SUBROUTINE OUTPUT  *
C     ***********************
C
      SUBROUTINE OUTPUT(V1,V2,P2,V3,ES1,NES2,ES3,ND1,ND2,ND3)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 V1(ND2,ND3,*),V3(ND2,ND3,*),V2(ND1,ND3,*),P2(ND1,ND3,*)
      DIMENSION ES1(*),NES2(*),ES3(*)
      COMMON /FREQ/ FT,NO / NFA / NFAI
     *       / NS / NS2
     *       / NE / NE1,NE2,NE3
      CHARACTER FILE1*10
C
      PAI2=6.283185307D0
C
      WRITE(5,4900)
 4900 FORMAT(/1X,'ENTER FILE NAME OF OUTPUT DATA : ',$)
      READ(5,4910) FILE1
 4910 FORMAT(A10)
C
      OPEN(7,FILE=FILE1,STATUS='NEW')
C
        WRITE(7,7000) FT,NO,NFAI
        WRITE(7,7010) NE1,NE2,NE3
        WRITE(7,7025) (ES1(I),I=1,NE1)
        WRITE(7,7020) (NES2(I),I=1,NE2)
        WRITE(7,7025) (ES3(I),I=1,NE3)
C
      DO 500 N=1,NFAI
C
      WRITE(6,6000)N
C
C     *** BASIN (U) ***
C
      DO 110 MM=1,NE1
        WRITE(6,6015)ES1(MM)
        M=MM*2-1
        DO 100 L=1,NO,5
          WRITE(6,6020) (FT*(L+J),PAI2/(FT*(L+J)),V1(M,L+J,N)
     *      ,CDABS(V1(M,L+J,N)),CDABS(V1(M,L+J,N))*FT*(L+J),J=0,4)
          WRITE(7,7030) (DREAL(V1(M,L+J,N)),
     *                   DIMAG(V1(M,L+J,N)),J=0,4)
  100   CONTINUE
  110 CONTINUE
C
C     *** BASIN (W) ***
C
      DO 160 MM=1,NE1
        WRITE(6,6016)ES1(MM)
        M=MM*2
        DO 150 L=1,NO,5
          WRITE(6,6020) (FT*(L+J),PAI2/(FT*(L+J)),V1(M,L+J,N)
     *      ,CDABS(V1(M,L+J,N)),CDABS(V1(M,L+J,N))*FT*(L+J),J=0,4)
          WRITE(7,7030) (DREAL(V1(M,L+J,N)),
     *                   DIMAG(V1(M,L+J,N)),J=0,4)
  150   CONTINUE
  160 CONTINUE
C
C     *** VERTICAL BOUNDARY (U) ***
C
      WRITE(6,6030)
      DO 210 MM=1,NE2
        M=2*NES2(MM)-1
        WRITE(6,6010)NES2(MM)
        DO 200 L=1,NO,5
          WRITE(6,6020) (FT*(L+J),PAI2/(FT*(L+J)),V2(M,L+J,N)
     *      ,CDABS(V2(M,L+J,N)),CDABS(V2(M,L+J,N))*FT*(L+J),J=0,4)
          WRITE(7,7030) (DREAL(V2(M,L+J,N)),
     *                   DIMAG(V2(M,L+J,N)),J=0,4)
  200   CONTINUE
  210 CONTINUE
C
C     *** VERTICAL BOUNDARY (W) ***
C
      DO 260 MM=1,NE2
        M=2*NES2(MM)
        WRITE(6,6011)NES2(MM)
        DO 250 L=1,NO,5
          WRITE(6,6020) (FT*(L+J),PAI2/(FT*(L+J)),V2(M,L+J,N)
     *      ,CDABS(V2(M,L+J,N)),CDABS(V2(M,L+J,N))*FT*(L+J),J=0,4)
          WRITE(7,7030) (DREAL(V2(M,L+J,N)),
     *                   DIMAG(V2(M,L+J,N)),J=0,4)
  250   CONTINUE
  260 CONTINUE
C
C     *** ROCK (U) ***
C
      WRITE(6,6040)
      DO 310 MM=1,NE3
        WRITE(6,6015)ES3(MM)
        M=2*MM-1
        DO 300 L=1,NO,5
          WRITE(6,6020) (FT*(L+J),PAI2/(FT*(L+J)),V3(M,L+J,N)
     *      ,CDABS(V3(M,L+J,N)),CDABS(V3(M,L+J,N))*FT*(L+J),J=0,4)
          WRITE(7,7030) (DREAL(V3(M,L+J,N)),
     *                   DIMAG(V3(M,L+J,N)),J=0,4)
  300   CONTINUE
  310 CONTINUE
C
C     *** ROCK (W) ***
C
      DO 360 MM=1,NE3
        WRITE(6,6016)ES3(MM)
        M=2*MM
        DO 350 L=1,NO,5
          WRITE(6,6020) (FT*(L+J),PAI2/(FT*(L+J)),V3(M,L+J,N)
     *      ,CDABS(V3(M,L+J,N)),CDABS(V3(M,L+J,N))*FT*(L+J),J=0,4)
          WRITE(7,7030) (DREAL(V3(M,L+J,N)),
     *                   DIMAG(V3(M,L+J,N)),J=0,4)
  350   CONTINUE
  360 CONTINUE
C
C     *** VERTICAL BOUNDARY (PX) ***
C
      WRITE(6,6050)
      DO 410 M=1,NE2
        MM=2*NES2(M)-1
        WRITE(6,6010)NES2(M)
        DO 400 L=1,NO,5
          WRITE(6,6025)(FT*(J+L),PAI2/(FT*(J+L)),P2(MM,L+J,N),
     *       CDABS(P2(MM,L+J,N)),CDABS(P2(MM,L+J,N))*FT*(L+J),J=0,4)
  400   CONTINUE
  410 CONTINUE
C
C     *** VERTICAL BOUNDARY (PZ) ***
C
      DO 460 M=1,NE2
        MM=2*NES2(M)
        WRITE(6,6011)NES2(M)
        DO 450 L=1,NO,5
          WRITE(6,6025)(FT*(J+L),PAI2/(FT*(J+L)),P2(MM,L+J,N),
     *       CDABS(P2(MM,L+J,N)),CDABS(P2(MM,L+J,N))*FT*(L+J),J=0,4)
  450   CONTINUE
  460 CONTINUE
C
  500 CONTINUE
C
      CLOSE(7)
C
      RETURN
C
 6000 FORMAT(/////5X,'NUMBER OF FAULT = ',I5, ///
     *          10X,'DISP OF FREE-SURFACE OF BASIN'/ )
 6010 FORMAT(  /18X,'U: ELEMENT NUMBER  =  ',I4//,5X
     * ,'OMEGA   PERIOD',7X,'REAL',8X,'IMAG',8X,'ABS',7X,'ABS*OMEGA'/)
 6011 FORMAT(  /18X,'W: ELEMENT NUMBER  =  ',I4//,5X
     * ,'OMEGA   PERIOD',7X,'REAL',8X,'IMAG',8X,'ABS',7X,'ABS*OMEGA'/)
 6015 FORMAT(  /18X,'U: ELEMENT POSITION =  ',F10.3//,5X
     * ,'OMEGA   PERIOD',7X,'REAL',8X,'IMAG',8X,'ABS',7X,'ABS*OMEGA'/)
 6016 FORMAT(  /18X,'W: ELEMENT POSITION =  ',F10.3//,5X
     * ,'OMEGA   PERIOD',7X,'REAL',8X,'IMAG',8X,'ABS',7X,'ABS*OMEGA'/)
 6020 FORMAT(    2X,2F9.4,4F12.4)
 6025 FORMAT(    2X,2F9.4,4F12.4)
 6030 FORMAT(///10X,'DISP OF VERICAL BOUNDARY'//)
 6040 FORMAT(///10X,'DISP OF FREE-SURFACE OF ROCK'//)
 6050 FORMAT(///10X,'STRESS OF VERTICAL BOUNDARY'//)
 7000 FORMAT(1X,F15.7,2I15)
 7010 FORMAT(1X,3I15)
 7020 FORMAT(1X,20I3)
 7025 FORMAT(1X,7F10.3)
 7030 FORMAT(1X,10F8.5)
C
      END
