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

Contents of /MITgcm/lsopt/cubic.F

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


Revision 1.1.2.1 - (show annotations) (download)
Tue Feb 5 20:34:33 2002 UTC (22 years, 10 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
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