C
C *** PHASE VELOCITY (CC) OF LOVE WAVE
C
      SUBROUTINE PHAS2(OM,NDC,MMD,CC,NMODE)
C
C ** INPUT DATA **
C OM    : CIRCULAR FREQUENCY (OMEGA)
C NDC   : DIVISION NUMVER FROM MIN VS TO MAX VS TO GET CC
C MMD   : NEMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA
C CC(J) : LOVE WAVE PHASE VELOCITY OF Jth MODE
C NMODE : TOTAL NUMBER OF PHASE VELOCITY
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/OMEGA/OM1
      COMMON/MATP/VS1,VS2,DN1,DN2,Q1,Q2 /DEPTH/DEP
     *      /NS/NS1,NS2,NS3,NS4,NS5
      DIMENSION CC(0:MMD)
      EXTERNAL FNFL2
C
C *** COMMON ***
C
      OM1=OM
C
C *** MINIMUM VS ***
C
      VSMIN=VS1
C
C *** INCREMENT OF CC (DC) BETWEEN VS(NL)(= MAX VS) AND MINIMUM VS ***
C
      BN=VS2
      DC=(BN-VSMIN)/DFLOAT(NDC)
C
C *** INPUT OF CONVERGENCE JUDGMENT CRITERIA FOR LVROOT ***
C
      DCEPS=100000.0D0
      CEPSO=DC/DCEPS
      DFEPS=100000.0D0
      ITMAX=100
C     WRITE(*,*)
C     WRITE(*,*)'From SUBROUTINE LVCC :'
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(FL1-FL2)/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 ***
C
C  ** LOWER TRIAL CC AND EIGEN VALUE EQUATION **
C
      CC1=VSMIN-DC
      FL1=FNFL2(CC1)
C
C  ** CHECK OF LOWER TRIAL CC **
C
      IFC=2
      ICC=0
      IF(FL1.EQ.0.0D0)THEN
        CC(0)=CC1
        CC1=VSMIN
        FL1=FNFL2(CC1)
        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=VSMIN+DFLOAT(IC-2)*DC
        FL2=FNFL2(CC2)
C
C   * CHECK OF TRIAL CC (LOWER AND HIGER CC) *
C
        IF(FL1*FL2)10,20,30
C
C     * CASE OF LOWER CC < CORRECT CC < HIGHER CC *
C
C       * SEARCH OF CORRECT CC BY LVROOT *
C
   10     CEPS=CEPSO
          FEPS=DABS(FL1-FL2)/DFEPS
          ITR=ITMAX
          CL=CC1
          CR=CC2
          CALL LVRT2(FNFL2,CL,CR,CEPS,FEPS,ITR,CCLR,IER)
C
C       * LVROOT ERROR OR NOT CONVERGENCE WORK *
C
          IF(IER.NE.0)THEN
            WRITE(*,*)
            WRITE(*,*)' FROM SUBROUTINE PHAS2 :'
            WRITE(*,*)' LVROOT ERROR, IER =',IER
            WRITE(*,*)' OMEGA =',OM
            WRITE(*,2000)CEPS,FEPS,ITR
            WRITE(*,*)' CC1,CC2 =',CC1,CC2
            WRITE(*,*)' FL1,FL2 =',FL1,FL2
            IF(IER.EQ.3000)THEN
              WRITE(*,*)' NOT CONVERGENCE (IER=3000) !!!'
              WRITE(*,*)' CL,CR =',CL,CR
              FLL=FNFL2(CL)
              FLR=FNFL2(CR)
              WRITE(*,*)' FLL,FLR =',FLL,FLR
              WRITE(*,*)' CC IS CHANGED BY MEAN CC !!!'
              CCLR=(CR+CL)/2.0D0
              FLLR=FNFL2(CCLR)
              WRITE(*,*)' NEW CC, FL =',CCLR,FLLR
            ELSE
              WRITE(*,*)' IER=2000 INPUT DATA ERROR  !!!'
              WRITE(*,*)' IER=2100 FL1*FL2 > 0       !!!'
              WRITE(*,*)' STOP  !!!'
              STOP
            ENDIF
          ENDIF
 2000     FORMAT(' CEPS,FEPS,ITRATION =',2D15.7,I5)
C
C       * CORRECT CC BY LVROOT *
C
          CC(ICC)=CCLR
          ICC=ICC+1
          GOTO 30
C
C    * CASE OF CORRECT CC = HIGHER CC *
C
   20     CC(ICC)=CC2
          ICC=ICC+1
          CC2=CC1+DC/2.0D0
          FL2=FNFL2(CC2)
C
C    * CASE OF LOWER CC < HIGER CC < CORRECT CC *
C
C      * LOWER CC ==> HIGHER CC AND NEXT STEP *
C
   30     CC1=CC2
          FL1=FL2
C
  110 CONTINUE
C
C  ** NUMBER OF MODE **
C
      NMODE=ICC-1
C
      RETURN
      END
