C
C *** APROXIMATE LOVE WAVE IN SEDIMENT BASIN (ASUUMED AT (0,0))
C     BY 2-D SH HASKELL MODEL ***
C
      PROGRAM MLVKIN
C
C ** MARAMETER **
C MOM   : NOM IN MAIN PROGRAM DIMENSION
C       : IF MOM > 2**13(=8192), CAN NOT USE IN SUBROUTINE !!!
C MMD   : NEMD IN MAIN PROGRAM DIMENSION
C ** INPUT DATA **
C DOM   : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM   : TOTAL NUMBER OF CIRCULAR FREQUENCY
C DNS(I): DENSITY OF Ith LAYER ==> SM(I):SHEAR MODULUS => COMMOM SM(35)
C VS (I): VS OF Ith LAYER                              => COMMON VS(35)
C DEP(I): DEPTH ( THICKNESS ) OF Ith LAYER             => COMMON DEP(35)
C NL    : TOTAL NUMBER OF LAYERS                       => COMMON NL
C NDC   : DIVISION NUMBER FROM MIN VS TO MAX VS TO PHASE VELOCITY
C NSMD  : START OUTPUT MODE OF LOVE WAVE
C NEMD  : END OUTPUT MODE OF LOVE WAVE
C NDOM  : DIVISION NUMBER FROM 0 TO MAX OMEGA TO CUT-OFF OMEGA
C XF    : FAULT LOCATION ON X-AXIS 
C ZF    : FAULT LOCATION ON Z-AXIS 
C XO    : LOCATION OF OBSERVATION POINT ON X-AXIS 
C ZO    : LOCATION OF OBSERVATION POINT ON Z-AXIS 
C W     : WIDTH OF FAULT
C DA    : DIP ANGLE OF FAULT (DEGREE--->RADIAN IN THIS SUBROUTINE)
C D     : DISPLACEMENT DIFFELENCE OF FAULT
C VRR   : RUPTURE VELOCITY OF FAULT
C RT    : RISE TIME
C ** CALCULATE DATA **
C CFOM(I) : CUT-OFF OMEGA OF Ith MODE LOVE WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C CC(I,J) : LOVE WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): OUTPUT HIGHEST LOVE MODE NUMBER AT Ith OMEGA
C UU(I,J) : GROUP VELOCITY OF Jth MODE LOVE WAVE AT Ith OMEGA
C ALK(I,J): AMPLITUDE SPEC (i*AL/K) OF Jth MODE LOVE WAVE AT Ith OMEGA
C         : TO OBTAIN AN ACTUAL AL/K, MULTIPLE AL/K BY (-i)
C VF(I)   : Y DISPLACEMENT OF Ith OMEGA AT (XO,ZO) BY FAULT
C VO(I,J) : V DISPLACEMENT OF Ith OMEGA AT Jth POINT ON BAISN 
C
      PARAMETER(MOM=100,MMD=10,MOB=10)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM),XOB(MOB)
     *         ,UU(MOM,0:MMD),ALK(MOM,0:MMD)
      COMPLEX*16 VF(MOM),VO(MOM,MOB),EX,VOA
      CHARACTER FNAME*15
      COMMON /LAYER/NL,VS(35),DEP(35),SM(35) /DNSITY/DNS(35)
C
C *** SOIL DATA INPUT ***
C
   10 CALL LOVEIN(DOM,NOM,NDOM,NSMD,NEMD,NDC)
C
C *** OBSERVATION POINT INPUT ***
C
      WRITE(*,*)
      WRITE(*,1000)
   20 WRITE(*,1010)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 20
      END IF
      IF(NOB.EQ.0)GOTO 20
C
      DO 100 I=1,NOB
        WRITE(*,1020)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
 1000 FORMAT(/' From MLVKIN : Basin edge is located at (0,0)'
     *       /' Approximate Lave Wave on Basin by 2-D Haskell'/)
 1010 FORMAT( ' Enter Number of Observation Point : ',\)
 1020 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',\)
C
C *** SEISMIC FAULT RARAMETER INPUT ***
C
      CALL FALT2D(XF,ZF,W,DA,D,VRR,RT)
C
C *** PARAMETER CHECK ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MLVKIN :'
        WRITE(*,*)'MOM < NOM THUS STOP !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        STOP
      END IF
C
      IF(MMD.LT.NEMD)THEN
        WRITE(*,*)'FROM MAIN MLVKIN :'
        WRITE(*,*)'MMD < NEMD THUS STOP !!!'
        WRITE(*,*)'MMD,NEMD =',MMD,NEMD
        STOP
      END IF
C
C *** INPUT OF DISPLACEMENT BY 2-D SH HASKELL ***
C
      XO=0.D0
      ZO=0.D0
      VSN=VS(NL)
      CALL HASKSH(DOM,NOM,MOM,XF,ZF,XO,ZO,VSN,W,DA,D,VRR,RT,VF)
C
C *** CUT-OFF OMEGA ***
C
       CALL LOVECT(DOM,NOM,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C *** PHASE VELOCITY ***
C
      CALL LOVECC(DOM,NOM,MOM,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
       CALL LOVEUA(DOM,NOM,MOM,CC,NMODE,CFOM,NSMD,MMD,UU,ALK)
C
C *** Calculation of approximate Love waves gererated at basin edge ***
C
      TOKA=2.D0/(1.D0+DNS(1)*VS(1)/(DNS(NL)*VS(NL)))
C
      DO 200 K=1,NOM
C
        OM=DOM*DFLOAT(K)
C
        DO 210 IOB=1,NOB
C
          VO(K,IOB)=(0.D0,0.D0)
C
          DO 220 MODE=NSMD,NCFOM
C
            IF(CC(K,MODE).EQ.0.D0)THEN
              WN=0.D0
              EX=(0.D0,0.D0)
            ELSE
              WN=OM/CC(K,MODE)
              EX=CDEXP(DCMPLX(0.D0,-OM*XOB(IOB)/CC(K,MODE)))
            END IF
C
            VOA=-4.D0*TOKA*SM(1)*DEP(1)*VF(K)*WN*ALK(K,MODE)*EX
            VO(K,IOB)=VO(K,IOB)+VOA
C
  220     CONTINUE
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** WAVE DATA OUTPUT TO FILE ***
C
      WRITE(*,*)
   30 WRITE(*,1200)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 30
      WRITE(*,1210)
 1200 FORMAT(/' Enter name of wave data output file for PLOTMS : ',\)
 1210 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VO roop')
C      
      OPEN(7,FILE=FNAME)
C
        WRITE(7,1300)DOM,NOM
        WRITE(7,1310)
C
        DO 300 IOB=1,NOB
C
          DO 310 K=1,NOM,5
            WRITE(7,1320)(VO(K+J,IOB),J=0,4)
  310     CONTINUE
C
  300   CONTINUE
C
      CLOSE(7)
C
 1300 FORMAT(F15.8,I15)
 1310 FORMAT(///)
 1320 FORMAT(10F8.5)
C
C *** CHECK OF DATA ***
C
      WRITE(*,1400)
      READ(*,*)NYES
      IF(NYES.NE.0)GOTO 40
 1400 FORMAT(/' Output Love wave data on screen ? (yes=0) :',\) 
C
C *** DATA OUTPUT ON SCREEN ***
C
      CALL LOVEOT(DOM,NOM,MOM,CFOM,NCFOM,CC,NSMD,MMD,UU,ALK)
C
   40 STOP
      END
C
C
C *** 2-D SH HASKELL'S FAULT MODEL ***
C
      SUBROUTINE HASKSH(DOM,NOM,MOM,XF,ZF,XO,ZO,VS,W,DA,D,VRR,RT,VF)
C
C ** MARAMETER **
C MOM : NOM IN MAIN PROGRAM DIMENSION
C     : IF MOM > 2**13(=8192), CAN NOT USE IN SUBROUTINE !!!
C ** INPUT DATA **
C DOM : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM : TOTAL NUMBER OF CIRCULAR FREQUENCY
C XF  : FAULT LOCATION ON X-AXIS 
C ZF  : FAULT LOCATION ON Z-AXIS 
C XO  : LOCATION OF OBSERVATION POINT ON X-AXIS 
C ZO  : LOCATION OF OBSERVATION POINT ON Z-AXIS 
C VS  : SHEAR WAVE VELOCITY
C W   : WIDTH OF FAULT
C DA  : DIP ANGLE OF FAULT (DEGREE--->RADIAN IN THIS SUBROUTINE)
C D   : DISPLACEMENT DIFFELENCE OF FAULT
C VRR : RUPTURE VELOCITY OF FAULT
C RT  : RISE TIME
C ** CALCULATE DATA **
C VF(I) : Y DISPLACEMENT OF Ith OMEGA AT (XO,ZO) BY FAULT
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 H1S,VF(MOM)
C
      PAI=3.141592654D0
      DA=DA*PAI/180.0D0
      RW=W/VRR
      CC=D*W/(4.D0*VS)
C
      RFX=XF-XO
      RFZ=ZF-ZO
      XN= DSIN(DA)
      ZN=-DCOS(DA)
      RF=DSQRT(RFX*RFX+RFZ*RFZ)
      RXZ=(RFX*XN+RFZ*ZN)/RF
C
      DO 100 K=1,NOM
C
        OM=DOM*DFLOAT(K)
        AS=OM*RF/VS
        Y1S=BESL21(AS)
        B1S=BESL11(AS)
        H1S=DCMPLX(B1S,-Y1S)
        XT=OM*RT/2.D0
        XW=OM*RW/2.D0
        ST=DSIN(XT)/XT
        SW=DSIN(XW)/XW
        VF(K)=CC*H1S*RXZ*ST*SW*CDEXP(DCMPLX(0.D0,-XT-XW))
C
  100 CONTINUE
C
      RETURN
      END
C
C
      REAL FUNCTION BESL11*8 (X)
C
C   FUNCTION            - BY CHEBYSHEV INTERPOLATION, WE OBTAIN THE
C                           VALUE OF THE BESSEL FUNCTION OF THE FIRST
C                           KIND OF ORDER ONE.
C
C   NOTICE              - WHEN !X!>=1.0D12, ERROR MESSAGE IS PRINTED
C                           OUT AND BESL11 IS SET TO 0.0.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION NAME(2)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9
     *   /  1.2249 4472           D-8, -3.9299 65711          D-7,
     *      8.9514 96872 0        D-6, -1.6124 40490 740      D-4,
     *      2.2574 92590 3842     D-3, -2.3703 70267 63209    D-2,
     *      1.7777 77775 58267 6  D-1, -8.8888 88888 61512 4  D-1,
     *      2.6666 66666 66491 50 D 0, -3.9999 99999 99995 63 D 0 /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14
     *   /  2.2193 1              D-12,  4.1163 4              D-12,
     *     -4.3816 515            D-10, -5.7599 652            D-10,
     *      6.2306 13774          D- 8,  4.8857 13666          D- 8,
     *     -6.1029 92681 89       D- 6, -1.7066 40506 78       D- 6,
     *      3.7656 50660 9458     D- 4, -6.7213 28811 703      D- 5,
     *     -1.2601 05985 12045 4  D- 2,  6.8295 50387 40438    D- 3,
     *      1.6844 60804 23507 80 D- 1, -1.1208 09437 96045 28 D- 1,
     *     -3.2757 91375 91465 22 D- 1                              /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14
     *   / -1.0746 3             D-12,  2.7780 03             D-11,
     *      2.0080 677           D-10, -4.7578 8198           D- 9,
     *     -2.6460 77937         D- 8,  5.7071 42412 9        D- 7,
     *      2.3218 13149 97      D- 6, -4.5285 96911 743      D- 5,
     *     -1.2058 28825 8765    D- 4,  2.1170 94601 37864    D- 3,
     *      2.9657 44869 54234   D- 3, -4.7160 51742 84081 7  D- 2,
     *     -1.9188 38969 35370 8 D- 2,  3.0074 82453 02747 75 D- 1,
     *     -4.6828 23482 34583   D- 3                              /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     *   / -3.0207 732            D-10,  2.0806 2912        D-9,
     *     -7.0003 5699           D- 9,  1.6557 68548       D-8,
     *     -3.5582 32648          D- 8,  8.9550 40706       D-8,
     *     -3.2722 78326 5        D- 7,  2.0593 18935 93    D-6,
     *     -2.8088 72167 333      D- 5,  1.4609 70265 11633 D-3,
     *      7.9788 45608 02865 46 D- 1                          /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     *   /  9.4650 569            D-10, -7.0536 9184        D-9,
     *      2.4646 12368          D- 8, -5.5142 32274       D-8,
     *      9.3432 85420          D- 8, -1.4024 48404 3     D-7,
     *      2.2297 66956 2        D- 7, -4.4752 03358 4     D-7,
     *      1.2957 31867 99       D- 6, -6.0676 75972 38    D-6,
     *      5.4070 79074 732      D- 5, -1.2783 48981 98997 D-3,
     *      2.9920 67103 01074 48 D- 1                          /
      DATA PI4  / 0.78539 81633 97448 31D0 /
      DATA NAME / '$DMJ', '1M  ' /
      AX=ABS(X)
      IF (AX .GE. 1.0D12) GO TO 50
      IF (AX .GE. 8.0D0) GO TO 30
      IF (AX .GE. 6.0D0) GO TO 20
      IF (AX .GE. 4.0D0) GO TO 10
C
C     !X! < 4
C
      Y=X*0.25D0
      Z=Y*Y
      BESL11=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z +2.0D0) *Y
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESL11=(((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z
     *       +D6)*Z +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z +D12)*Z
     *       +D13)*Z +D14
      GO TO 40
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESL11=(((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z
     *       +E13)*Z +E14
      GO TO 40
C
C     8 <= !X! < 1.0D12
C
   30 Z=(8.0D0/AX)**2
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z +A7)*Z
     *  +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *  +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=AX-PI4
      BESL11=(A*SIN(Q) + B*COS(Q)) /SQRT(AX)
   40 IF (X .LT. 0.0D0) BESL11=-BESL11
      GO TO 9999
C
C     PARAMETER ERROR
C
   50 BESL11=0.0D0
      IER=2000
C
 9999 RETURN
      END
C
C
      REAL FUNCTION BESL21*8 (X)
C
C   FUNCTION            - BY CHEBYSHEV INTERPOLATION, WE OBTAIN THE
C                           VALUE OF THE BESSEL FUNCTION OF THE SECOND
C                           KIND OF ORDER ONE.
C
C   NOTICE              - WHEN !X! DOES NOT SATISFY 8.797D-77<=!X!<
C                           1.0D12, ERROR MESSAGE IS PRINTED OUT AND
C                           $DMY0M IS SET TO 0.0 FOR !X!>=1.0D12 OR
C                           -1.0D75 FOR !X!<8.797D-77.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION NAME(2)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9
     *   /  1.2249 4472           D-8, -3.9299 65711          D-7,
     *      8.9514 96872 0        D-6, -1.6124 40490 740      D-4,
     *      2.2574 92590 3842     D-3, -2.3703 70267 63209    D-2,
     *      1.7777 77775 58267 6  D-1, -8.8888 88888 61512 4  D-1,
     *      2.6666 66666 66491 50 D 0, -3.9999 99999 99995 63 D 0 /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11
     *   /  2.3353 113            D-9, -8.7149 8135           D-8,
     *      2.3316 55657 5        D-6, -5.0097 83252 95       D-5,
     *      8.5332 26991 090      D-4, -1.1176 62851 97065    D-2,
     *      1.0801 30262 83770 8  D-1, -7.2710 13435 63856 6  D-1,
     *      3.1262 10900 22178 23 D 0, -7.3980 37853 32564 88 D 0,
     *      6.8529 24964 20573 60 D 0,  3.9325 61521 02602 95 D-1 /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14,E15
     *   /  4.8468 2              D-12, -2.2800 90             D-11,
     *      6.3132 56             D-11, -5.2574 110            D-10,
     *      7.8189 4200           D- 9, -5.9240 4185           D- 9,
     *     -5.9222 98307 9        D- 7, -6.1872 86946 2        D- 7,
     *      5.2495 83490 136      D- 5, -4.1370 37951 76       D- 6,
     *     -2.3535 74442 27453    D- 3,  7.6458 54114 1893     D- 4,
     *      5.3923 88944 88126 4  D- 2, -3.7165 28343 50646 6  D- 2,
     *     -3.3809 02539 27279 20 D- 1,  1.4786 31433 91226 85 D- 1 /
      DATA F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14
     *   /  1.7873 7              D-12,  1.6011 35            D-11,
     *     -3.8677 585            D-10, -2.3735 6163          D- 9,
     *      5.4322 38759          D- 8,  2.6412 21207 5       D- 7,
     *     -5.3968 07752 06       D- 6, -1.7970 21894 410     D- 5,
     *      3.3358 85263 4720     D- 4,  6.7436 67691 6589    D- 4,
     *     -1.1206 70091 19909 2  D- 2, -9.4701 49780 74827   D- 3,
     *      1.4701 02892 64358 67 D- 1,  1.7288 43275 05313 8 D- 2,
     *     -3.0266 72370 24184 87 D- 1                             /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     *   / -3.0207 732            D-10,  2.0806 2912        D-9,
     *     -7.0003 5699           D- 9,  1.6557 68548       D-8,
     *     -3.5582 32648          D- 8,  8.9550 40706       D-8,
     *     -3.2722 78326 5        D- 7,  2.0593 18935 93    D-6,
     *     -2.8088 72167 333      D- 5,  1.4609 70265 11633 D-3,
     *      7.9788 45608 02865 46 D- 1                          /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     *   /  9.4650 569            D-10, -7.0536 9184        D-9,
     *      2.4646 12368          D- 8, -5.5142 32274       D-8,
     *      9.3432 85420          D- 8, -1.4024 48404 3     D-7,
     *      2.2297 66956 2        D- 7, -4.4752 03358 4     D-7,
     *      1.2957 31867 99       D- 6, -6.0676 75972 38    D-6,
     *      5.4070 79074 732      D- 5, -1.2783 48981 98997 D-3,
     *      2.9920 67103 01074 48 D- 1                          /
      DATA PI2,PI4
     *   /  0.63661 97723 67581 34D0, 0.78539 81633 97448 31D0 /
      DATA NAME / '$DMY', '1M  ' /
      AX=ABS(X)
      IF (AX .LT. 8.797D-77) GO TO 60
      IF (AX .GE. 1.0D12) GO TO 50
      IF (AX .GE. 8.0D0) GO TO 30
      IF (AX .GE. 6.0D0) GO TO 20
      IF (AX .GE. 4.0D0) GO TO 10
C
C     0 < !X! < 4
C
      Y=X*0.25D0
      Z=Y*Y
      BJ1=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z +2.0D0) *Y
      ADD=((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z +D6)*Z
     *       +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z -PI2) /X
      BESL21=PI2*BJ1*LOG(AX*0.5D0) +ADD
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESL21=((((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z
     *       +E13)*Z +E14)*Z +E15
      GO TO 40
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESL21=(((((((((((((F0*Z +F1)*Z +F2)*Z +F3)*Z +F4)*Z +F5)*Z
     *       +F6)*Z +F7)*Z +F8)*Z +F9)*Z +F10)*Z +F11)*Z +F12)*Z
     *       +F13)*Z +F14
      GO TO 40
C
C     8 <= !X! < 1.0D12
C
   30 Z=(8.0D0/AX)**2
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z +A7)*Z
     *       +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *       +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=AX-PI4
      BESL21=(B*SIN(Q) - A*COS(Q)) /SQRT(AX)
   40 IF (X .LT. 0.0D0) BESL21=-BESL21
      GO TO 9999
C
C     PARAMETER ERROR
C
   50 BESL21=0.0D0
      GO TO 70
C
   60 BESL21=-1.0D75
   70 IER=2000
C
 9999 RETURN
      END
