MITgcm/lsopt/cubic.F


      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
1
2	subroutine cubic( t, f, fp, ta, fa, fpa, tlower, tupper )
3
4	c-----------------------------------------
5	c arguments
6	c-----------------------------------------
7	double precision t, f, fp, ta, fa, fpa, tlower, tupper
8
9	c-----------------------------------------
10	c local variables
11	c-----------------------------------------
12	double precision sign, den, anum
13	double precision z1, b, discri
14
15	c-----------------------------------------
16	c Using f and fp at t and ta,
17	c computes new t by cubic formula
18	c safeguarded inside [tlower,tupper].
19	c-----------------------------------------
20	z1 = dble(fp) + dble(fpa) - 3.d0*dble(fa-f)/dble(ta-t)
21	b = z1 + dble(fp)
22
23	c-----------------------------------------
24	c first compute the discriminant
25	c (without overflow)
26	c-----------------------------------------
27	if (abs(z1).le.1.) then
28	discri = z1z1-dble(fp)dble(fpa)
29	if (discri .lt. 0.d0) then
30	if (fp.lt.0.) t = tupper
31	if (fp.ge.0.) t = tlower
32	return
33	else
34	discri = dsqrt(discri)
35	end if
36	else
37	discri = dble(fp)/z1
38	discri = discri*dble(fpa)
39	discri = z1-discri
40	if (z1.ge.0.d0 .and. discri.ge.0.d0) then
41	discri = dsqrt(z1)*dsqrt(discri)
42	else if (z1.le.0.d0 .and. discri.le.0.d0) then
43	discri = dsqrt(-z1)*dsqrt(-discri)
44	else
45	if (fp.lt.0.) t = tupper
46	if (fp.ge.0.) t = tlower
47	return
48	end if
49	end if
50
51	c-----------------------------------------
52	c discriminant nonnegative,
53	c compute solution (without overflow)
54	c-----------------------------------------
55	if (t-ta .lt. 0.0) then
56	discri = -discri
57	end if
58
59	sign = (t-ta)/abs(t-ta)
60	if (sngl(b)*sign .gt. 0.0) then
61	t = t + fp*(ta-t)/sngl(b+discri)
62	else
63	den = sngl(z1+b+dble(fpa))
64	anum = sngl(b-discri)
65	if (abs((t-ta)anum).lt.(tupper-tlower)abs(den)) then
66	t = t + anum*(ta-t)/den
67	else
68	t = tupper
69	end if
70	end if
71
72	t = max( t, tlower )
73	t = min( t, tupper )
74
75	return
76	end