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

Annotation of /MITgcm/lsopt/cubic.F

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

Revision 1.1.2.1 - (hide annotations) (download)
Tue Feb 5 20:34:33 2002 UTC (22 years, 3 months ago) by heimbach
Branch: ecco-branch
CVS Tags: ecco_c50_e32, ecco_c50_e33, ecco_c50_e30, ecco_c50_e31, ecco_c51_e34d, ecco_c51_e34e, ecco_c51_e34f, ecco_c51_e34g, ecco_c51_e34a, ecco_c51_e34b, ecco_c51_e34c, icebear5, icebear4, icebear3, icebear2, ecco_c50_e29, ecco_c50_e28, ecco_c44_e19, ecco_c44_e18, ecco_c44_e17, ecco_c44_e16, ecco_c50_e33a, ecco_c51_e34, ecco_ice2, ecco_ice1, ecco_c44_e25, ecco_c44_e22, ecco_c44_e23, ecco_c44_e20, ecco_c44_e21, ecco_c44_e26, ecco_c44_e27, ecco_c44_e24, ecco-branch-mod1, ecco-branch-mod2, ecco-branch-mod3, ecco-branch-mod4, ecco-branch-mod5
Branch point for: c24_e25_ice, icebear
Changes since 1.1: +76 -0 lines 
o Updating adjoint/makefile to ECCO code
o Adding optim and lsopt for line search optimization.
o Adding verif. experiments for ECCO
Code will be tagged ecco-branch-mod1.
1 heimbach 1.1.2.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
  ViewVC Help
Powered by ViewVC 1.1.22