C
      PROGRAM M
C
      PARAMETER(MOM=200,MOB=21)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VS1,VS2,SM1,SM2,WB1,WB2,SS,FNC,VV(MOM,MOB)
      DIMENSION XOB(MOB)
      COMMON /DATA/ OM,X,WB1,WB2,SM1,SM2,THK
      CHARACTER FNAME*15
      EXTERNAL FNC
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)
 1000 FORMAT(/' From M : SOURCE POINT is located at (0,0)'
     * /' EXACT GREEN FUNTION IN 2-LAYER SH MEDIUM')
C
      WRITE(*,*)'Q1 ?'
      READ(*,*)Q1
      WRITE(*,*)'Q2 ?'
      READ(*,*)Q2
      WRITE(*,*)'DN1 ?'
      READ(*,*)DN1
      WRITE(*,*)'DN2 ?'
      READ(*,*)DN2
      WRITE(*,*)'V1  ?'
      READ(*,*)V1 
      WRITE(*,*)'V2  ?'
      READ(*,*)V2 
      WRITE(*,*)'THK ?'
      READ(*,*)THK
      WRITE(*,1010)Q1,DN1,V1,THK
      WRITE(*,1020)Q2,DN2,V2
 1010 FORMAT(/' Q1,DN1,V1,THK=',4F8.1)
 1020 FORMAT( ' Q2,DN2,V2    =',4F8.1)
C
C *** DATA INPUT ***
C
   10 WRITE(*,*)'OMEGA NUMBER?'
      READ(*,*)NOM
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN M'
        WRITE(*,*)'MOM < NOM RE-ENTRER !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        GOTO 10
      END IF
C
      WRITE(*,*)'FIRST OM?'
      READ(*,*)FOM
      WRITE(*,*)'DELTA OM?'
      READ(*,*)DOM
C
   20 WRITE(*,1110)
      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(*,1120)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
 1110 FORMAT(/' Enter Number of Observation Point : ',\)
 1120 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',\)
C
      WRITE(*,1210)
      READ(*,*)SCC
      WRITE(*,1220)
      READ(*,*)NS
C
 1210 FORMAT(/' Enter the lowest phase velocity : ',\)
 1220 FORMAT( ' Enter div. num. for intergral   : ',\)
C
C *** DATA CAL ***
C
      VS1=V1*DCMPLX(1.D0,0.5D0/Q1)
      VS2=V2*DCMPLX(1.D0,0.5D0/Q2)
      SM1=DN1*VS1*VS1
      SM2=DN2*VS2*VS2
C
C *** GREEN FUNC ***
C
      DO 200 IOB=1,NOB
C
        X=XOB(IOB)
C
        DO 210 IOM=1,NOM
C
          OM=FOM+DOM*DFLOAT(IOM-1)
          WB1=OM*OM/(VS1*VS1)
          WB2=OM*OM/(VS2*VS2)
C
          WNA=0.D0
          WNB=OM/SCC
          IF(OM.EQ.0.D0)WRITE(*,*)'OM=0 !!!!!!!!!!!!!!!!!!!!'
C
          CALL NEWTON(FNC,WNA,WNB,NS,SS,IER)
          IF(IER.NE.0)WRITE(*,*)' IOB,IOM,IER',IOB,IOM,IER,'!!!!!!!!!!'
C
          VV(IOM,IOB)=SS/3.1415927D0
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** OUTPUT OF DATA ***
C
        WRITE(*,1400)DOM,NOM
        DO 300 IOB=1,NOB
        WRITE(*,1410)XOB(IOB)
          DO 310 K=1,NOM,4
            WRITE(*,1420)(VV(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
C
 1400 FORMAT(/' DOM,NOM=',F6.2,I5)
 1410 FORMAT(/' X=',F10.2)
 1420 FORMAT(8D10.3)
C
C *** OUTPUT OF DATA FILE ***
C
      WRITE(*,*)
   30 WRITE(*,1310)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 30
      WRITE(*,1320)
 1310 FORMAT( ' Enter file name : ',\)
 1320 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VV roop')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1500)DOM,NOM
        WRITE(8,1510)
        DO 500 IOB=1,NOB
          DO 510 K=1,NOM,5
            WRITE(8,1520)(VV(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 1500 FORMAT(F15.8,I15)
 1510 FORMAT(///)
 1520 FORMAT(10F8.5)
C
      STOP
      END
