C
C *** ROOT FINDING OF FUNC(X)= 0 BY TWO DIVISION AND LINEAR
C     INTERPOLATION COMBINED METHOD
C *** SOURCE PROGRAM IS $DNTBM (MSL2) AND ARRENGED IN 1988
C
      SUBROUTINE LVROOT(FUNC,ML,NL,DNS,VS,DEP
     *                 ,XL,XR,XEPS,FEPS,ITMAX,X,IER)
C
C ** INPUT DATA **
C FUNC  : FUNCTION OF X (NEED EXTERNAL IN MAIN PROGRAM !!!)
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 XL    : LEFT(LOWER LIMIT) X (XL<XS)
C XR    : RIGHT(HIGHER LIMIT) X (XS<XR)
C XEPS  : CRITERION VALUE TO JUGDE CONVERGENCE OF X (XEPS>0)
C FEPS  : CRITERION VALUE TO JUGDE CONVERGENCE OF FUNC(X,*) (FEPS>0)
C ITMAX : MAXIMUM ITERATION NUMBER
C ** OUTPUT DATA **
C ITMAX : FINAL ITERATION NUMBER IN THE SUBROUTINE
C X     : ROOT (SOLUSION) OF FUNC(X,ML,NL,DNS,VS,DEP)=0
C IER   : RETURN CODE (IF ERROR, IER.NE.0.  SEE $DNTBM,MSL2)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),VS(ML),DEP(ML)
C
C *** INPUT DATA CHECK ***
C
      IER = 0
      IF(XEPS.LT.0.0D0 .OR. FEPS.LT.0.0D0 .OR. ITMAX.LT.1) GO TO 260
C
      N   = 1
      X   = XL
      FL  = FUNC(X,ML,NL,DNS,VS,DEP)
      IF( FL ) 100,555,100
  100 X   = XR
      FR  =  FUNC(X,ML,NL,DNS,VS,DEP)
      IF( FR  ) 110,555,110
  110 IF( FR*FL ) 130,555,120
  120 X   = (XL+XR)*.5D0
  125 CONTINUE
      IER=2100
      N   = ITMAX
      GO TO 555
C
C *** ROOT FINDING OF FUNC(X,ML,NL,DNS,VS,DEP)= 0 ***
C
  130 DO  250 N=1,ITMAX
C
C   * TRIAL X AND FUNC(X,ML,NL,DNS,VS,DEP) BY TWO DIVISION METHOD *
C
        X   = (XL+XR)*.5D0
        F   = FUNC(X,ML,NL,DNS,VS,DEP)
C
C   * REPLACE OF XL OR XR BY TRIAL X *
C
        IF( F ) 140,555,140
C
  140     IF( F*FL ) 150,555,160
C
  150       XR  = X
            FR  = F
            GO TO 170
  160       XL  = X
            FL  = F
C
C   * TRIAL X AND FUNC BY LINEAR INTERPOLATION METHOD *
C
  170   DF  = FL-FR
        DX  = (XR-XL)*FL
        DX  = DX/DF
        X   = XL+DX
        F   = FUNC(X,ML,NL,DNS,VS,DEP)
C
C   * REPLACE OF XL OR XR BY TRIAL X *
C
        IF( F ) 175,555,175
C
  175     IF( FL*F ) 180,190,190
C
  180       XR  = X
            FR  = F
            GO TO 200
  190       XL  = X
            FL  = F
C
C   * CONVERGENCE JUDGMENT OF TRIAL X AND FUNC(X,ML,NL,DNS,VS,DEP) *
C
  200   SX  = XR-XL
        IF( ABS(X).LE.1.0D0 ) GO TO 210
          SX  = SX/X
  210   IF( ABS(SX)-XEPS ) 220,220,250
  220     X   = XR
        IF( ABS(FR)-FEPS ) 555,555,230
  230     X   = XL
        IF( ABS(FL)-FEPS ) 555,555,250
C
  250 CONTINUE
C
C *** ITERATION OVER ERROR ***
C
      IER=3000
      N = ITMAX
      GO TO 555
C
C *** INPUT DATA ERROR ***
C
  260 IER=2000
      N = ITMAX
C
C *** FINAL WORK ***
C
  555 CONTINUE
      ITMAX = N
C
      RETURN
      END

