--- MITgcm/lsopt/cubic.F	2002/02/05 20:34:33	1.1
+++ MITgcm/lsopt/cubic.F	2002/11/15 04:03:24	1.2
@@ -0,0 +1,76 @@
+
+      subroutine cubic( t, f, fp, ta, fa, fpa, tlower, tupper )
+
+c-----------------------------------------
+c arguments
+c-----------------------------------------
+      double precision    t, f, fp, ta, fa, fpa, tlower, tupper
+
+c-----------------------------------------
+c local variables
+c-----------------------------------------
+      double precision sign, den, anum
+      double precision z1, b, discri
+
+c-----------------------------------------
+c Using f and fp at t and ta,
+c computes new t by cubic formula
+c safeguarded inside [tlower,tupper].
+c-----------------------------------------
+      z1 = dble(fp) + dble(fpa) - 3.d0*dble(fa-f)/dble(ta-t)
+      b  = z1 + dble(fp)
+
+c-----------------------------------------
+c first compute the discriminant
+c (without overflow)
+c-----------------------------------------
+      if (abs(z1).le.1.) then
+         discri = z1*z1-dble(fp)*dble(fpa)
+         if (discri .lt. 0.d0) then
+            if (fp.lt.0.) t = tupper
+            if (fp.ge.0.) t = tlower
+            return
+         else
+            discri = dsqrt(discri)
+         end if
+      else
+         discri = dble(fp)/z1
+         discri = discri*dble(fpa)
+         discri = z1-discri
+         if (z1.ge.0.d0 .and. discri.ge.0.d0) then
+            discri = dsqrt(z1)*dsqrt(discri)
+         else if (z1.le.0.d0 .and. discri.le.0.d0) then
+            discri = dsqrt(-z1)*dsqrt(-discri)
+         else
+            if (fp.lt.0.) t = tupper
+            if (fp.ge.0.) t = tlower
+            return
+         end if
+      end if
+
+c-----------------------------------------
+c discriminant nonnegative,
+c compute solution (without overflow)
+c-----------------------------------------
+      if (t-ta .lt. 0.0) then
+         discri = -discri
+      end if
+
+      sign = (t-ta)/abs(t-ta)
+      if (sngl(b)*sign .gt. 0.0) then
+         t    = t + fp*(ta-t)/sngl(b+discri)
+      else
+         den  = sngl(z1+b+dble(fpa))
+         anum = sngl(b-discri)
+         if (abs((t-ta)*anum).lt.(tupper-tlower)*abs(den)) then
+            t = t + anum*(ta-t)/den
+         else
+            t = tupper
+         end if
+      end if
+
+      t = max( t, tlower )
+      t = min( t, tupper )
+      
+      return
+      end