c    @(#)erate.f	1.2
      subroutine erate(us,usbs,h,hmax,be,ier)
c
c
c
c  5/22/84 leonard walstad
c
c   routine to calculate the entrainment rate based on garwood's
c   mixed layer model.  
c
c
c
c  us      friction velocity in m/s
c
c  usbs    integrated buoyancy flux in the mixed layer due to 
c          surface processes. 
c
c  h       mixed layer depth in meters.
c
c  be      entrainment rate in m**2/s**3.  if less than zero, then
c          detrainment is occuring and the rate is not valid.
c
c====================================================================
c
c
      real m1, m2, m3, m4, m5, m3us3
c
      common /garwood/  f0, m1, m2, m3, m4, m5
c
      data icall/0/, usl/1000./, usbsl/1000./, hlast/1000./
c
c
      if(icall.ne.0)go to 100
        icall = 1
        a1 = m2 + .666666 * m1
        ta2 = .666666 * m5 * f0
        a3 =  - 3.0 * m2
        b1 =  - .5 * m4
        b2 = m2 - .333333 * m1
        tb3 =  - .333333 * m5 * f0
        c1 = 3.0 * m2 + 2.0 * m1
        tc2 = 1.33333 * m5 * f0
        c3 =  - 6.0 * m2
        c4 = 3.0 * m2 - m1
        c5 =  - m4
        tc6 =  - .666666 * m5 * f0
        d1=b2
        d2=a1
        d3=3.0*d1
        d4=2.0*tb3
        d5=3.0*d2
        d6=-2.0*tc6
        d7=sqrt(1.5/m5/f0)
c
        qzfluxb=(m3/3.0/m2)**.3333333
        hzfluxb=-b2*qzfluxb/tb3
        d8=b2/a1
        d9=(a1+2.0*b2)**.6666666
        d10=(1.0+d8+d8)
        call cubert(-1.333333,3.333333*d8,d8,0.0,s1,s2,s3,ier)
        if(ier.lt.0)stop'erate1'
c
        z1=1000.
        if(s1.le.-.5 .or. s1.ge.0.0)go to 51
        z1=s1*d10*(1.0+2.0*s1)**.6666666/d9/(d8-s1)
   51 continue
        z2=1000.
        if(s2.le.-.5 .or. s2.ge.0.0)go to 52
        z2=s2*d10*(1.0+2.0*s2)**.6666666/d9/(d8-s2)
   52 continue
        z3=1000.
        if(s3.le.-.5 .or. s3.ge.0.0)go to 53
        z3=s3*d10*(1.0+2.0*s3)**.6666666/d9/(d8-s3)
   53 continue
c
        if(z1.gt.0.0 .and. z2.gt.0.0 .and. z3.gt.0.0)stop 'erate2'
c
        z1=min(z1,z2,z3)
        fmin=.6666666*f0*(m3**.66666666)*m5*z1
c
        rnd=1e-10
        rnd2=1e-15
c
c
c
c
  100 continue
      ier=0
c
c --- determine wether or not there is a solution to the problem.
c
      be=0.0
      ier=-2
      if(us.lt.0.0)return
      usbs2=usbs/2.0
      m3us3 = m3 * us * us * us
      emax=.00001*(m3us3+abs(h*.5*usbs))
      if(usbs.ne.0.0)go to 110
        hm=hzfluxb*us
        qm=qzfluxb*us
        go to 190
c
c
  110 continue
      if(usbs.lt.0.0)go to 120
c       if(abs(usbs-usbsl)/(abs(usbsl)+rnd2) .lt. .25   .and.
c    .     abs(us-usl)/(abs(usl)+rnd) .lt. .25)go to 150
        hm=hzfluxb*us
        qm=qzfluxb*us
        if(us.eq.0.0)go to 190
        go to 150
c
c
  120 continue
        bsmin=us*us*fmin
        hm=10000.
        if(usbs.lt.bsmin)go to 190
c       if(abs(usbs-usbsl)/(abs(usbsl)+rnd2) .lt. .25   .and.
c    .     abs(us-usl)/(abs(usl)+rnd) .lt. .25)go to 150
        qm=qzfluxb*us
        hm=hzfluxb*us
        go to 150
c
c
  150 continue
      ncnt=0
c
  151 continue
      q2=qm*qm
      q3=qm*q2
      hmq=hm*qm
      hmq2=hm*q2
c
      e1=d2*q3-tc6*hmq2-m3us3
      e2=d1*q3+tb3*hmq2-usbs2*hm
c
      if(abs(e1).le.emax .and. abs(e2).le.emax)go to 190
c
c
      e1q=d5*q2+d6*hmq
      e1h=-tc6*q2
      e2q=d3*q2+d4*hmq
      e2h=tb3*q2-usbs2
c
      det=e1q*e2h-e1h*e2q
      dq=(-e1*e2h+e1h*e2)/det
      dh=(-e1q*e2+e1*e2q)/det
c
      qm=qm+dq
      hm=hm+dh
      ncnt=ncnt+1
c
      if(ncnt.lt.20)go to 151
      ier=-1
      return
c
c
c
  190 continue
      ier=-1
      usl=-10000.
      if(qm.lt.0.0 .or. hm.lt.0.0)return
      ier=0
      be=0.0
      hmax=hm
      if(hmax.lt.0.0)hmax=0.0
      if(h.lt.hm)go to 200
      usbsl=usbs
      usl=us
        be=0.0
        return
c
c
c
c
  200 continue
      a2 = ta2 * h
      b3 = tb3 * h
      c2 = tc2 * h
      c6 = tc6 * h
      r=0.0
      husbs = .5 * h * usbs
      ncnt = 0
c
c
c     if(abs(usbs-usbsl)/(abs(usbsl)+rnd2) .lt. .25   .and.
c    .   abs(us-usl)/(abs(usl)+rnd) .lt. .25  .and.
c    .   abs(h-hlast)/(abs(hlast)+rnd) .lt. .25)go to 201
      q=(m3us3/a1)**.3333333
      r=.3333333*q + (abs(usbs))**.3333333
      if(q.lt.r)q=r
c
  201 continue
      r2 = r*r
c
      e1 = ((a1 * q + a2) * q + a3 * r2) * q -m3us3
      e2 = ((b1 * r + b3 + b2 * q) * q + a3 *r2) * q - husbs
c
      if(abs(e1).le.emax .and. abs(e2).le.emax) go to 300
c
      e1q = (c1 * q + c2)*q + a3 * r2
      e1r = c3 * q *  r
      e2q = (c5 * r + c6 + c4 * q) * q + a3 * r2
      e2r = (b1 * q + c3 * r) * q
c
      det = e1q * e2r - e1r * e2q
c
      dq = ( - e1 * e2r + e1r * e2) / det
      dr = ( - e1q * e2 + e1 * e2q) / det
c
      q = q + dq
      r = r + dr
      ncnt = ncnt+1
c
c
      if(ncnt.lt.20)go to 201
      be = 0.0
      usl=10000.
      ier = -1
      return
c
c
c
  300 continue
      be = m4 * r * q * q / h 
      usl=us
      usbsl=usbs
      hlast=h
      return
      end
