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

Contents of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 139 - (show annotations)
Tue May 26 17:46:03 2015 UTC (9 years ago) by guez
Original Path: trunk/Sources/dyn3d/fxhyp.f
File size: 6843 byte(s) 
dynetat0 read rlonu, rlatu, rlonv, rlatv, cu_2d, cv_2d, aire_2d from
"start.nc" and then these variables were overwritten by
inigeom. Corrected this. Now, inigeom does not compute rlonu, rlatu,
rlonv and rlatv. Moreover, cu_2d, cv_2d, aire_2d are not written to
"restart.nc". Since xprimu, xprimv, xprimm025, xprimp025, rlatu1,
rlatu2, yprimu1, yprimu2 are computed at the same time as rlonu,
rlatu, rlonv, rlatv, and since it would not be convenient to separate
those computations, we decide to write xprimu, xprimv, xprimm025,
xprimp025, rlatu1, rlatu2, yprimu1, yprimu2 into "restart.nc", read
them from "start.nc" and not compute them in inigeom. So, in summary,
"start.nc" contains all the coordinates and their derivatives, and
inigeom only computes the 2D-variables.

Technical details:

Moved variables rlatu, rlonv, rlonu, rlatv, xprimu, xprimv from module
comgeom to module dynetat0_m. Upgraded local variables rlatu1,
yprimu1, rlatu2, yprimu2, xprimm025, xprimp025 of procedure inigeom to
variables of module dynetat0_m.

Removed unused local variable yprimu of procedure inigeom and
corresponding argument yyprimu of fyhyp.

Moved variables clat, clon, grossismx, grossismy, dzoomx, dzoomy,
taux, tauy from module serre to module dynetat0_m (since they are read
from "start.nc"). The default values are now defined in read_serre
instead of in the declarations. Changed name of module serre to
read_serre_m, no more module variable here.

The calls to fxhyp and fyhyp are moved from inigeom to etat0.

Side effects in programs other than gcm: etat0 and read_serre write
variables of module dynetat0; the programs test_fxyp and
test_inter_barxy need more source files.

Removed unused arguments len and nd of cv3_tracer. Removed unused
argument PPSOL of LWU.

Bug fix in test_inter_barxy: forgotten call to read_serre.
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 dynetat0_m, only: clon, grossismx, dzoomx, taux
22 use invert_zoom_x_m, only: invert_zoom_x, nmax
23 use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24 use principal_cshift_m, only: principal_cshift
25
26 REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
27 real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
28
29 ! Local:
30 real rlonm025(iim + 1), rlonp025(iim + 1)
31 REAL dzoom, step
32 real d_rlonv(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 fa, fb
37 INTEGER i, is2
38 DOUBLE PRECISION xmoy, fxm
39
40 !----------------------------------------------------------------------
41
42 print *, "Call sequence information: fxhyp"
43
44 test_grossismx: if (grossismx == 1.) then
45 step = twopi / iim
46
47 xprimm025(:iim) = step
48 xprimp025(:iim) = step
49 xprimv(:iim) = step
50 xprimu(:iim) = step
51
52 rlonv(:iim) = arth(- pi + clon, step, iim)
53 rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54 rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55 rlonu(:iim) = rlonv(:iim) + 0.5 * step
56 else test_grossismx
57 dzoom = dzoomx * twopi_d
58 xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
59
60 ! Compute fhyp:
61 DO i = nmax, 2 * nmax
62 fa = taux * (dzoom / 2. - xtild(i))
63 fb = xtild(i) * (pi_d - xtild(i))
64
65 IF (200. * fb < - fa) THEN
66 fhyp(i) = - 1.
67 ELSE IF (200. * fb < fa) THEN
68 fhyp(i) = 1.
69 ELSE
70 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
71 IF (200. * fb + fa < 1e-10) THEN
72 fhyp(i) = - 1.
73 ELSE IF (200. * fb - fa < 1e-10) THEN
74 fhyp(i) = 1.
75 END IF
76 ELSE
77 fhyp(i) = TANH(fa / fb)
78 END IF
79 END IF
80
81 IF (xtild(i) == 0.) fhyp(i) = 1.
82 IF (xtild(i) == pi_d) fhyp(i) = -1.
83 END DO
84
85 ! Calcul de beta
86
87 ffdx = 0.
88
89 DO i = nmax + 1, 2 * nmax
90 xmoy = 0.5 * (xtild(i-1) + xtild(i))
91 fa = taux * (dzoom / 2. - xmoy)
92 fb = xmoy * (pi_d - xmoy)
93
94 IF (200. * fb < - fa) THEN
95 fxm = - 1.
96 ELSE IF (200. * fb < fa) THEN
97 fxm = 1.
98 ELSE
99 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
100 IF (200. * fb + fa < 1e-10) THEN
101 fxm = - 1.
102 ELSE IF (200. * fb - fa < 1e-10) THEN
103 fxm = 1.
104 END IF
105 ELSE
106 fxm = TANH(fa / fb)
107 END IF
108 END IF
109
110 IF (xmoy == 0.) fxm = 1.
111 IF (xmoy == pi_d) fxm = -1.
112
113 ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
114 END DO
115
116 print *, "ffdx = ", ffdx
117 beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
118 print *, "beta = ", beta
119
120 IF (2. * beta - grossismx <= 0.) THEN
121 print *, 'Bad choice of grossismx, taux, dzoomx.'
122 print *, 'Decrease dzoomx or grossismx.'
123 STOP 1
124 END IF
125
126 ! calcul de Xprimt
127 Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
128 xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
129
130 ! Calcul de Xf
131
132 DO i = nmax + 1, 2 * nmax
133 xmoy = 0.5 * (xtild(i-1) + xtild(i))
134 fa = taux * (dzoom / 2. - xmoy)
135 fb = xmoy * (pi_d - xmoy)
136
137 IF (200. * fb < - fa) THEN
138 fxm = - 1.
139 ELSE IF (200. * fb < fa) THEN
140 fxm = 1.
141 ELSE
142 fxm = TANH(fa / fb)
143 END IF
144
145 IF (xmoy == 0.) fxm = 1.
146 IF (xmoy == pi_d) fxm = -1.
147 xxpr(i) = beta + (grossismx - beta) * fxm
148 END DO
149
150 xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
151
152 Xf(0) = - pi_d
153
154 DO i=1, 2 * nmax - 1
155 Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
156 END DO
157
158 Xf(2 * nmax) = pi_d
159
160 call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
161 xuv = - 0.25d0)
162 call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
163 xuv = 0d0)
164 call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
165 xuv = 0.5d0)
166 call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
167 xuv = 0.25d0)
168 end if test_grossismx
169
170 is2 = 0
171
172 IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
173 .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
174 IF (clon <= 0.) THEN
175 is2 = 1
176
177 do while (rlonm025(is2) < - pi .and. is2 < iim)
178 is2 = is2 + 1
179 end do
180
181 if (rlonm025(is2) < - pi) then
182 print *, 'Rlonm025 plus petit que - pi !'
183 STOP 1
184 end if
185 ELSE
186 is2 = iim
187
188 do while (rlonm025(is2) > pi .and. is2 > 1)
189 is2 = is2 - 1
190 end do
191
192 if (rlonm025(is2) > pi) then
193 print *, 'Rlonm025 plus grand que pi !'
194 STOP 1
195 end if
196 END IF
197 END IF
198
199 call principal_cshift(is2, rlonm025, xprimm025)
200 call principal_cshift(is2, rlonv, xprimv)
201 call principal_cshift(is2, rlonu, xprimu)
202 call principal_cshift(is2, rlonp025, xprimp025)
203
204 forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
205 print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
206 print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
207
208 ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
209 DO i = 1, iim + 1
210 IF (rlonp025(i) < rlonv(i)) THEN
211 print *, 'rlonp025(', i, ') = ', rlonp025(i)
212 print *, "< rlonv(", i, ") = ", rlonv(i)
213 STOP 1
214 END IF
215
216 IF (rlonv(i) < rlonm025(i)) THEN
217 print *, 'rlonv(', i, ') = ', rlonv(i)
218 print *, "< rlonm025(", i, ") = ", rlonm025(i)
219 STOP 1
220 END IF
221
222 IF (rlonp025(i) > rlonu(i)) THEN
223 print *, 'rlonp025(', i, ') = ', rlonp025(i)
224 print *, "> rlonu(", i, ") = ", rlonu(i)
225 STOP 1
226 END IF
227 END DO
228
229 END SUBROUTINE fxhyp
230
231 end module fxhyp_m
  ViewVC Help
Powered by ViewVC 1.1.21