C
C *** CUT-OFF OMEGA OF LOVE WAVE (C=VS(NL))
C
      SUBROUTINE LVCT(DOM,NOM,ML,NL,DNS,VS,DEP,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C ** INPUT DATA
C DOM  : INCREMENT OF CIRCULAR FREQUENCY (DELTA OMEGA)
C NOM  : TOTAL NUMBER OF CIRCULAR FREQUENCY
C ML    : NL IN MAIN 
C NL    : TOTAL NUMBER OF LAYERS           
C DNS(I): DENSITY OF Ith LAYER 
C VS (I): VS OF Ith LAYER      
C DEP(I): DEPTH (THICKNESS) OF Ith LAYER  
C NDOM : DIVISION NUMBER from 0 to MAX OMEGA TO GET CUT-OFF OMEGA
C NEMD : HIGHEST OUTPUT MODE OF LOVE WAVE
C MMD  : NMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA
C CFOM(I): CUT-OFF OMEGA OF Ith MODE LOVE WAVE
C NCFOM  : FINAL MODE
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),DNS(ML),VS(ML),DEP(ML)
      EXTERNAL FA21
C
C *** 0 th MODE CUT-OFF OMEGA ***
C
      ICO=0
      CFOM(ICO)=0.0D0
      IF(NEMD.EQ.0)GOTO 40
C
C *** INCREMENT OF OMEGA ***
C
      DDOM=DFLOAT(NOM)*DOM/DFLOAT(NDOM)
C
C *** INPUT OF CONVERGENCE JUDGMENT CRITERIA FOR LVROOT ***
C
      DXEPS=100.0D0
      XEPSO=DDOM/DXEPS
      DFEPS=1000.0D0
      ITMAX=100
      WRITE(*,*)
C     WRITE(*,*)'From SUBROUTINE LVCT :'
C     WRITE(*,*)
C     WRITE(*,*)'Convergence Criteria for Cut-Off Omega :'
C     WRITE(*,1000)XEPSO
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(*,*)XEPSO,DFEPS,ITMAX
C     END IF
C     WRITE(*,*)
C1000 FORMAT(/' XEPS = DDOM/DXEPS          : XEPS =',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 DDOM, DXEPS, DFEPS, ITMAX :',\)
C
C *** LOWER TRIAL OMEGA AND EIGEN VALUE EQUATION ***
C
      OM1=DDOM
      FL1=FA21(OM1,ML,NL,DNS,VS,DEP)
C
C *** LOWER TRIAL OMEGA CHECK ***
C
      ICO=1
      IFOM=2
      IF(FL1.EQ.0.0D0)THEN
        CFOM(1)=OM1
        IF(NEMD.EQ.1)GOTO 40
        OM1=2.0D0*DDOM
        FL1=FA21(OM1,ML,NL,DNS,VS,DEP)
        ICO=2
        IFOM=3
      ENDIF
C
C *** SEARCH OF CUT-OFF OMEGA OF HIGHER LOVE MODE ***
C
      DO 100 IOM=IFOM,NDOM
C
C  ** HIGHER TRIAL OMEGA AND EIGEN VALUE EQUATION **
C
        OM2=DFLOAT(IOM)*DDOM
        FL2=FA21(OM2,ML,NL,DNS,VS,DEP)
C
C  ** LOWER AND HIGER OMEGA CHECK **
C
        IF(FL1*FL2)10,20,30
C
C   * CASE OF LOWER OM < CUT-OFF OM < HIGHER OM *
C
C     * SEARCH OF CUT-OFF OM BY D2ROOT *
C
   10     XEPS=XEPSO
          FEPS=DABS(FL1-FL2)/DFEPS
          ITR=ITMAX
          OML=OM1
          OMR=OM2
          CALL LVROOT(FA21,ML,NL,DNS,VS,DEP
     *               ,OML,OMR,XEPS,FEPS,ITR,OMLR,IER)
C
C     * D2ROOT ERROR OR NOT CONVERGENCE WORK *
C
          IF(IER.NE.0)THEN
            WRITE(*,*)
            WRITE(*,*)' FROM SUBROUTINE LVCT :'
            WRITE(*,*)' D2ROOT ERROR, IER =',IER
            WRITE(*,*)' XEPS,FEPS,ITMAX =',XEPS,FEPS,ITR
            WRITE(*,*)' OM1,OM2 =',OM1,OM2
            WRITE(*,*)' FL1,FL2 =',FL1,FL2
            IF(IER.EQ.3000)THEN
              WRITE(*,*)' NOT CONVERGENCE (IER=3000) !!!'
              WRITE(*,*)' OML,OMR =',OML,OMR
              FLL=FA21(OML,ML,NL,DNS,VS,DEP)
              FLR=FA21(OMR,ML,NL,DNS,VS,DEP)
              WRITE(*,*)' FLL,FLR =',FLL,FLR
              WRITE(*,*)' OMEGA IS CHANGED BY MEAN OMEGA !!!'
              OMLR=(OMR+OML)/2.0D0
              FLLR=FA21(OMLR,ML,NL,DNS,VS,DEP)
              WRITE(*,*)' NEW OM, FL =',OMLR,FLLR
            ELSE
              WRITE(*,*)' IER=2000 INPUT DATA ERROR  !!!'
              WRITE(*,*)' IER=2100 FL1*FL2 > 0       !!!'
              WRITE(*,*)' STOP  !!!'
              STOP
            ENDIF
          ENDIF
C
C     * ENTER OF CUT-OFF OMEGA *
C
          CFOM(ICO)=OMLR
          IF(NEMD.EQ.ICO)GOTO 40
          ICO=ICO+1
          GOTO 30
C
C   * CASE OF HIGER TRIAL OM = CUT-OFF OM *
C
   20     CFOM(ICO)=OM2
          IF(NEMD.EQ.ICO)GOTO 40
          ICO=ICO+1
          OM2=OM1+DDOM/2.0D0
          FL2=FA21(OM2,ML,NL,DNS,VS,DEP)
C
C   * CASE OF LOWER OM < HIGHER OM < CUT-OFF OM *
C
C     * HIGHER OM ==> LOWER OM AND NEXT STEP *
C
   30     OM1=OM2
          FL1=FL2
C
  100 CONTINUE
C
C *** NUMBER OF CUT-OFF OMEGA ***
C
      NCFOM=ICO-1
      RETURN
C
   40 NCFOM=ICO
C
      RETURN
      END
