C
C *** APROXIMATE LOVE WAVE IN SEDIMENT BASIN (ASUUMED AT (0,0))
C     BY 2-D SH HASKELL MODEL ***
C
      PROGRAM MHASKF
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 XF    : FAULT LOCATION ON X-AXIS 
C ZF    : FAULT LOCATION ON Z-AXIS 
C XO(I) : LOCATION OF OBSERVATION POINT ON X-AXIS OF Ith OMEGA
C ZO(I) : LOCATION OF OBSERVATION POINT ON Z-AXIS OF Ith OMEGA
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 VDO(I,J): WAVES BY SH HASKELL SOURCE OF Ith OMEGA AT Jth 
C
      PARAMETER(MOM=100,MOB=8)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(MOB),ZO(MOB)
      COMPLEX*16 VFO(MOM,MOB),VF
      CHARACTER FNAME*15
C
C *** INPUT OF FAULT RARAMETERS FROM FILE ***
C
      CALL H2DIN(XF,ZF,W,DA,D,VRR,RT)
C
C *** INPUT OF DOM, NOM AND OBSERVATION POINTS ON BASIN ***
C
      WRITE(*,*)
      WRITE(*,1000)
      WRITE(*,1010)
   10 READ(*,*)DOM,NOM
      IF(NOM.GT.MOM)THEN
        WRITE(*,*)'NOM > MOM !!! : Re-Enter NOM : MOM=',MOM
        GOTO 10
      END IF
      WRITE(*,1020)
      READ(*,*)VS
   20 WRITE(*,1030)
      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(*,1040)I
        READ(*,*)XO(I),ZO(I)
  100 CONTINUE
C
 1000 FORMAT(/' From MHASKF : wave by 2-D SH Haskell source'
     * /' Fundamental solusion in infinite medium is used'/)
 1010 FORMAT( ' Enter delta omega and omega number : ',\)
 1020 FORMAT( ' Enter VS (m/sec)  : ',\)
 1030 FORMAT( ' Enter number of observation point  : ',\)
 1040 FORMAT( ' Enter X and Z of OB of Num. =',I3,' : ',\)
C
C *** CALCULATION OF  ***
C
      DO 200 K=1,NOM
        OM=DOM*DFLOAT(K)
        DO 210 IOB=1,NOB
          CALL HKSHF1(OM,XF,ZF,XO(IOB),ZO(IOB),VS,W,DA,D,VRR,RT,VF)
          VFO(K,IOB)=VF
  210   CONTINUE
  200 CONTINUE
C
C *** OUTPUT OF DATA TO FILE ***
C
      WRITE(*,*)
   30 WRITE(*,1210)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 30
      WRITE(*,1220)
 1210 FORMAT( ' Enter file name of SH Haskell wave (PLOTMS) : ',\)
 1220 FORMAT( '  : Data format for PLOTMS 1:OB => 2: VFO roop')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1300)DOM,NOM
        WRITE(8,1310)
        DO 300 IOB=1,NOB
          DO 310 K=1,NOM,5
            WRITE(8,1320)(VFO(K+J,IOB),J=0,4)
  310     CONTINUE
  300   CONTINUE
      CLOSE(8)
C
 1300 FORMAT(F15.8,I15)
 1310 FORMAT(///)
 1320 FORMAT(10F8.5)
C
      STOP
      END
