/[lmdze]/trunk/dyn3d/fxhyp.f
ViewVC logotype

Contents of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (show annotations)
Thu Feb 5 12:41:08 2015 UTC (9 years, 3 months ago) by guez
File size: 5327 byte(s)
Added some test programs.

In fxhyp_loop_ik, changed precision from 1e-3 to 1e-6. Reset initial
value of xo1 to first guess of xi instead of final value of xi for
previous i (better logic).

1 module fxhyp_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8
9 ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10 ! Author: P. Le Van, from formulas by R. Sadourny
11
12 ! Calcule les longitudes et dérivées dans la grille du GCM pour
13 ! une fonction f(x) à dérivée tangente hyperbolique.
14
15 ! Il vaut mieux avoir : grossismx \times dzoom < pi
16
17 ! Le premier point scalaire pour une grille regulière (grossismx =
18 ! 1., taux=0., clon=0.) est à - 180 degrés.
19
20 USE dimens_m, ONLY: iim
21 use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22 use nr_util, only: pi_d, twopi_d, arth
23 use serre, only: clon, grossismx, dzoomx, taux
24
25 REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
26 real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
27
28 ! Local:
29
30 real rlonm025(iim + 1), rlonp025(iim + 1)
31 REAL dzoom
32 DOUBLE PRECISION xlon(iim)
33 DOUBLE PRECISION xtild(0:2 * nmax)
34 DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35 DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36 DOUBLE PRECISION xzoom, fa, fb
37 INTEGER i
38 DOUBLE PRECISION xmoy, fxm
39 DOUBLE PRECISION decalx
40
41 !----------------------------------------------------------------------
42
43 print *, "Call sequence information: fxhyp"
44
45 dzoom = dzoomx * twopi_d
46 xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
47
48 ! Compute fhyp:
49 DO i = nmax, 2 * nmax
50 fa = taux * (dzoom / 2. - xtild(i))
51 fb = xtild(i) * (pi_d - xtild(i))
52
53 IF (200. * fb < - fa) THEN
54 fhyp(i) = - 1.
55 ELSE IF (200. * fb < fa) THEN
56 fhyp(i) = 1.
57 ELSE
58 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
59 IF (200. * fb + fa < 1e-10) THEN
60 fhyp(i) = - 1.
61 ELSE IF (200. * fb - fa < 1e-10) THEN
62 fhyp(i) = 1.
63 END IF
64 ELSE
65 fhyp(i) = TANH(fa / fb)
66 END IF
67 END IF
68
69 IF (xtild(i) == 0.) fhyp(i) = 1.
70 IF (xtild(i) == pi_d) fhyp(i) = -1.
71 END DO
72
73 ! Calcul de beta
74
75 ffdx = 0.
76
77 DO i = nmax + 1, 2 * nmax
78 xmoy = 0.5 * (xtild(i-1) + xtild(i))
79 fa = taux * (dzoom / 2. - xmoy)
80 fb = xmoy * (pi_d - xmoy)
81
82 IF (200. * fb < - fa) THEN
83 fxm = - 1.
84 ELSE IF (200. * fb < fa) THEN
85 fxm = 1.
86 ELSE
87 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
88 IF (200. * fb + fa < 1e-10) THEN
89 fxm = - 1.
90 ELSE IF (200. * fb - fa < 1e-10) THEN
91 fxm = 1.
92 END IF
93 ELSE
94 fxm = TANH(fa / fb)
95 END IF
96 END IF
97
98 IF (xmoy == 0.) fxm = 1.
99 IF (xmoy == pi_d) fxm = -1.
100
101 ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
102 END DO
103
104 print *, "ffdx = ", ffdx
105 beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
106 print *, "beta = ", beta
107
108 IF (2. * beta - grossismx <= 0.) THEN
109 print *, 'Bad choice of grossismx, taux, dzoomx.'
110 print *, 'Decrease dzoomx or grossismx.'
111 STOP 1
112 END IF
113
114 ! calcul de Xprimt
115 Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
116 xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
117
118 ! Calcul de Xf
119
120 DO i = nmax + 1, 2 * nmax
121 xmoy = 0.5 * (xtild(i-1) + xtild(i))
122 fa = taux * (dzoom / 2. - xmoy)
123 fb = xmoy * (pi_d - xmoy)
124
125 IF (200. * fb < - fa) THEN
126 fxm = - 1.
127 ELSE IF (200. * fb < fa) THEN
128 fxm = 1.
129 ELSE
130 fxm = TANH(fa / fb)
131 END IF
132
133 IF (xmoy == 0.) fxm = 1.
134 IF (xmoy == pi_d) fxm = -1.
135 xxpr(i) = beta + (grossismx - beta) * fxm
136 END DO
137
138 xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
139
140 Xf(0) = - pi_d
141
142 DO i=1, 2 * nmax - 1
143 Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
144 END DO
145
146 Xf(2 * nmax) = pi_d
147
148 IF (grossismx == 1.) THEN
149 decalx = 1d0
150 else
151 decalx = 0.75d0
152 END IF
153
154 xzoom = clon * pi_d / 180d0
155 call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
156 xprimm025, xuv = - 0.25d0)
157 call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
158 xuv = 0d0)
159 call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
160 xuv = 0.5d0)
161 call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
162 xprimp025, xuv = 0.25d0)
163
164 print *
165
166 forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)
167 print *, "Minimum longitude step:", MINval(xlon) * 180. / pi_d, "degrees"
168 print *, "Maximum longitude step:", MAXval(xlon) * 180. / pi_d, "degrees"
169
170 DO i = 1, iim + 1
171 IF (rlonp025(i) < rlonv(i)) THEN
172 print *, 'rlonp025(', i, ') = ', rlonp025(i)
173 print *, "< rlonv(", i, ") = ", rlonv(i)
174 STOP 1
175 END IF
176
177 IF (rlonv(i) < rlonm025(i)) THEN
178 print *, 'rlonv(', i, ') = ', rlonv(i)
179 print *, "< rlonm025(", i, ") = ", rlonm025(i)
180 STOP 1
181 END IF
182
183 IF (rlonp025(i) > rlonu(i)) THEN
184 print *, 'rlonp025(', i, ') = ', rlonp025(i)
185 print *, "> rlonu(", i, ") = ", rlonu(i)
186 STOP 1
187 END IF
188 END DO
189
190 END SUBROUTINE fxhyp
191
192 end module fxhyp_m

  ViewVC Help
Powered by ViewVC 1.1.21