/[MITgcm]/MITgcm/lsopt/cubic.F
ViewVC logotype

Diff of /MITgcm/lsopt/cubic.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by heimbach, Tue Feb 5 20:34:33 2002 UTC revision 1.1.2.1 by heimbach, Tue Feb 5 20:34:33 2002 UTC
# Line 0  Line 1 
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 = z1*z1-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

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.2.1

  ViewVC Help
Powered by ViewVC 1.1.22