MITgcm/lsopt/cubic.F

4cee17c1be Patr*0001
                0002       subroutine cubic( t, f, fp, ta, fa, fpa, tlower, tupper )
                0003
                0004 c-----------------------------------------
                0005 c arguments
                0006 c-----------------------------------------
                0007       double precision    t, f, fp, ta, fa, fpa, tlower, tupper
                0008
                0009 c-----------------------------------------
                0010 c local variables
                0011 c-----------------------------------------
                0012       double precision sign, den, anum
                0013       double precision z1, b, discri
                0014
                0015 c-----------------------------------------
                0016 c Using f and fp at t and ta,
                0017 c computes new t by cubic formula
                0018 c safeguarded inside [tlower,tupper].
                0019 c-----------------------------------------
                0020       z1 = dble(fp) + dble(fpa) - 3.d0*dble(fa-f)/dble(ta-t)
                0021       b  = z1 + dble(fp)
                0022
                0023 c-----------------------------------------
                0024 c first compute the discriminant
                0025 c (without overflow)
                0026 c-----------------------------------------
                0027       if (abs(z1).le.1.) then
                0028          discri = z1*z1-dble(fp)*dble(fpa)
                0029          if (discri .lt. 0.d0) then
                0030             if (fp.lt.0.) t = tupper
                0031             if (fp.ge.0.) t = tlower
                0032             return
                0033          else
                0034             discri = dsqrt(discri)
                0035          end if
                0036       else
                0037          discri = dble(fp)/z1
                0038          discri = discri*dble(fpa)
                0039          discri = z1-discri
                0040          if (z1.ge.0.d0 .and. discri.ge.0.d0) then
                0041             discri = dsqrt(z1)*dsqrt(discri)
                0042          else if (z1.le.0.d0 .and. discri.le.0.d0) then
                0043             discri = dsqrt(-z1)*dsqrt(-discri)
                0044          else
                0045             if (fp.lt.0.) t = tupper
                0046             if (fp.ge.0.) t = tlower
                0047             return
                0048          end if
                0049       end if
                0050
                0051 c-----------------------------------------
                0052 c discriminant nonnegative,
                0053 c compute solution (without overflow)
                0054 c-----------------------------------------
                0055       if (t-ta .lt. 0.0) then
                0056          discri = -discri
                0057       end if
                0058
                0059       sign = (t-ta)/abs(t-ta)
                0060       if (sngl(b)*sign .gt. 0.0) then
                0061          t    = t + fp*(ta-t)/sngl(b+discri)
                0062       else
                0063          den  = sngl(z1+b+dble(fpa))
                0064          anum = sngl(b-discri)
                0065          if (abs((t-ta)*anum).lt.(tupper-tlower)*abs(den)) then
                0066             t = t + anum*(ta-t)/den
                0067          else
                0068             t = tupper
                0069          end if
                0070       end if
                0071
                0072       t = max( t, tlower )
                0073       t = min( t, tupper )
                0074
                0075       return
                0076       end