C
C *** PHASE VELOCITY (CC) OF SURFACE WAVE ***
C
      SUBROUTINE PHSE(FSCC,DOM,NOM,MOM,ML,NL,DNS,VP,VS,THK
     *               ,VMAX,VMIN,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C ** INPUT DATA
C DOM    : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM    : TOTAL NUMBER OF CIRCULAR FREQUENCY
C MOM    : NOM IN MAIN PROGRAM DIMENSION
C        : IF MOM > 2**13(=8192), CAN NOT USE !!!
C NDC    : DIVISION NUMVER FROM MIN CC TO MAX CC TO GET ANY CC
C NSMD   : START OUTPUT MODE OF RAYLEIGH WAVE
C NEMD   : END OUTPUT MODE OF RAYLEIGH WAVE
C CFOM(I): CUT-OFF OMEGA OF Ith MODE
C MMD    : NEMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA
C CC(I,J) : RAYLEIGH WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): TOTAL NUMBER OF PHASE VELOCITY AT Ith OMEGA
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),THK(ML),DNS(ML)
     *         ,CC(MOM,0:MMD),CFOM(0:MMD),NMODE(MOM)
      COMMON / OMEGA / OM 
C
      WRITE(*,*)
C
C *** INCREMENT OF CC (DC) BETWEEN MAX VS AND MIN CC ***
C
      DC=(VMAX-VMIN)/DFLOAT(NDC)
C
C *** INPUT OF CONVERGENCE JUDGMENT CRITERIA FOR ROOTN ***
C
      DCEPS=1000.0D0
      CEPSO=DC/DCEPS
      DFEPS=100000.0D0
      ITMAX=100
C     WRITE(*,*)
C     WRITE(*,*)'From SUBROUTINE PHSE :'
C     WRITE(*,*)
C     WRITE(*,*)'Convergence Criteria for Phase Velocity :'
C     WRITE(*,1000)CEPSO
C     WRITE(*,1010)DFEPS
C     WRITE(*,1020)ITMAX
C     WRITE(*,1030)
C     READ(*,*)IYES
C     IF(IYES.NE.0)THEN
C       WRITE(*,*)
C       WRITE(*,1040)
C       READ(*,*)CEPSO,DFEPS,ITMAX
C     END IF
C     WRITE(*,*)
C     WRITE(*,*)
C1000 FORMAT(/' CEPS = DC/DCEPS            : CEPS  =',F15.10)
C1010 FORMAT( ' FEPS = DABS(FS1-FS2)/DFEPS : DFEPS =',F15.3)
C1020 FORMAT( ' Maximum Iteration Number   : ITMAX =',I15/)
C1030 FORMAT( ' Are These Criteria OK ? ( yes = 0 ) :',\)
C1040 FORMAT( ' Enter New CEPS, DFEPS, ITMAX :',\)
C
C *** PHASE VELOCITY FOR EACH OMEGA ***
C
      DO 100 IOM=1,NOM
C
C  ** OMEGA AND THE LOWEST USED OMEGA **
C
        OM=DFLOAT(IOM)*DOM
        IF(OM.LT.CFOM(NSMD))GOTO 100
C
C  ** LOWER TRIAL CC AND EIGEN VALUE EQUATION **
C
        CC1=VMIN-DC
        FS1=FSCC(CC1,ML,NL,DNS,VP,VS,THK)
C
C  ** CHECK OF LOWER TRIAL CC **
C
        IFC=2
        ICC=0
        IF(FS1.EQ.0.0D0)THEN
          CC(IOM,0)=CC1
          IF(NEMD.EQ.0)GOTO 40
          CC1=VMIN
          FS1=FSCC(CC1,ML,NL,DNS,VP,VS,THK)
          IFC=3
          ICC=1
        ENDIF
C
C  ** SEARCH OF CORRECT CC FROM MINIMUM VS TO MAXIMUM VS(=VS(NL)) **
C
        DO 110 IC=IFC,NDC+2
C
C   * HIGHER TRIAL CC AND EIGEN VALUE EQUARION *
C
          CC2=VMIN+DFLOAT(IC-2)*DC
          FS2=FSCC(CC2,ML,NL,DNS,VP,VS,THK)
C
C   * CHECK OF TRIAL CC ( LOWER AND HIGER CC ) *
C
          IF(FS1*FS2)10,20,30
C
C     * CASE OF LOWER CC < CORRECT CC < HIGHER CC *
C
C       * SEARCH OF CORRECT CC BY ROOTN *
C
   10     CEPS=CEPSO
          FEPS=DABS(FS1-FS2)/DFEPS
          ITR=ITMAX
          CL=CC1
          CR=CC2
          CALL ROOTN(FSCC,ML,NL,DNS,VP,VS,THK
     *              ,CL,CR,CEPS,FEPS,ITR,CCLR,IER)
C
C       * ROOTN ERROR OR NOT CONVERGENCE WORK *
C
          IF(IER.NE.0)THEN
            WRITE(*,*)
            WRITE(*,*)' FROM SUBROUTINE PHSE :'
            WRITE(*,*)' ROOTN ERROR, IER =',IER
            WRITE(*,*)' OMEGA =',OM
            WRITE(*,2000)CEPS,FEPS,ITR
            WRITE(*,*)' CC1,CC2 =',CC1,CC2
            WRITE(*,*)' FS1,FS2 =',FS1,FS2
            IF(IER.EQ.3000)THEN
              WRITE(*,*)' NOT CONVERGENCE (IER=3000) !!!'
              WRITE(*,*)' CL,CR =',CL,CR
              FSL=FSCC(CL,ML,NL,DNS,VP,VS,THK)
              FSR=FSCC(CR,ML,NL,DNS,VP,VS,THK)
              WRITE(*,*)' FSL,FSR =',FSL,FSR
              WRITE(*,*)' CC IS CHANGED BY MEAN CC !!!'
              CCLR=(CR+CL)/2.0D0
              FSLR=FSCC(CCLR,ML,NL,DNS,VP,VS,THK)
              WRITE(*,*)' NEW CC, FS =',CCLR,FSLR
            ELSE
              WRITE(*,*)' IER=2000 INPUT DATA ERROR  !!!'
              WRITE(*,*)' IER=2100 FS1*FS2 > 0       !!!'
              WRITE(*,*)' STOP  !!!'
              STOP
            ENDIF
          ENDIF
 2000     FORMAT(' CEPS,FEPS,ITRATION =',2D15.7,I5)
C
C       * CORRECT CC BY ROOTN *
C
          CC(IOM,ICC)=CCLR
          IF(NEMD.EQ.ICC)GOTO 40
          ICC=ICC+1
          GOTO 30
C
C    * CASE OF CORRECT CC = HIGHER CC *
C
   20     CC(IOM,ICC)=CC2
          IF(NEMD.EQ.ICC)GOTO 40
          ICC=ICC+1
          CC2=CC1+DC/2.0D0
          FS2=FSCC(CC2,ML,NL,DNS,VP,VS,THK)
C
C    * CASE OF LOWER CC < HIGER CC < CORRECT CC *
C
C      * LOWER CC ==> HIGHER CC AND NEXT STEP *
C
   30     CC1=CC2
          FS1=FS2
C
  110   CONTINUE
C
C  ** MODE NUMBER AT EACH OMEGA **
C
        NMODE(IOM)=ICC-1
        GOTO 50
   40   NMODE(IOM)=NEMD
C
C *** IMFORMATION OUTPUT ON SCREEN ***
C
   50   WRITE(*,3000)OM,DOM*DFLOAT(NOM)
 3000   FORMAT('+Now C of R(L) Wave at OMEGA / Final OMEGA =',2F8.3)
C
C *** NEXT OMEGA ***
C
  100 CONTINUE
C
      WRITE(*,*)
C
      RETURN
      END
