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

Contents of /trunk/dyn3d/comgeom.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 265 - (show annotations)
Tue Mar 20 09:35:59 2018 UTC (6 years, 1 month ago) by guez
File size: 14170 byte(s)
Rename module dimens_m to dimensions.
1 module comgeom
2
3 use dimensions, only: iim, jjm
4
5 implicit none
6
7 private iim, jjm
8
9 real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m
10 real cu((iim + 1) * (jjm + 1)), cv((iim + 1) * jjm) ! in m
11 equivalence (cu, cu_2d), (cv, cv_2d)
12
13 real unscu2_2d(iim + 1, jjm + 1) ! in m-2
14 real unscu2((iim + 1) * (jjm + 1)) ! in m-2
15 equivalence (unscu2, unscu2_2d)
16
17 real unscv2_2d(iim + 1, jjm) ! in m-2
18 real unscv2((iim + 1) * jjm) ! in m-2
19 equivalence (unscv2, unscv2_2d)
20
21 real aire((iim + 1) * (jjm + 1)), aire_2d(iim + 1, jjm + 1) ! in m2
22 real airesurg_2d(iim + 1, jjm + 1), airesurg((iim + 1) * (jjm + 1))
23 equivalence (aire, aire_2d), (airesurg, airesurg_2d)
24
25 real aireu_2d(iim + 1, jjm + 1) ! in m2
26 real aireu((iim + 1) * (jjm + 1)) ! in m2
27 equivalence (aireu, aireu_2d)
28
29 real airev((iim + 1) * jjm), airev_2d(iim + 1, jjm) ! in m2
30 real unsaire((iim + 1) * (jjm + 1)), unsaire_2d(iim + 1, jjm + 1) ! in m-2
31 equivalence (airev, airev_2d), (unsaire, unsaire_2d)
32
33 real apoln, apols ! in m2
34
35 real unsairez_2d(iim + 1, jjm)
36 real unsairez((iim + 1) * jjm)
37 equivalence (unsairez, unsairez_2d)
38
39 real alpha1_2d(iim + 1, jjm + 1)
40 real alpha1((iim + 1) * (jjm + 1))
41 equivalence (alpha1, alpha1_2d)
42
43 real alpha2_2d(iim + 1, jjm + 1)
44 real alpha2((iim + 1) * (jjm + 1))
45 equivalence (alpha2, alpha2_2d)
46
47 real alpha3_2d(iim + 1, jjm + 1), alpha4_2d(iim + 1, jjm + 1)
48 real alpha3((iim + 1) * (jjm + 1)), alpha4((iim + 1) * (jjm + 1))
49 equivalence (alpha3, alpha3_2d), (alpha4, alpha4_2d)
50
51 real alpha1p2_2d(iim + 1, jjm + 1)
52 real alpha1p2((iim + 1) * (jjm + 1))
53 equivalence (alpha1p2, alpha1p2_2d)
54
55 real alpha1p4_2d(iim + 1, jjm + 1), alpha2p3_2d(iim + 1, jjm + 1)
56 real alpha1p4((iim + 1) * (jjm + 1)), alpha2p3((iim + 1) * (jjm + 1))
57 equivalence (alpha1p4, alpha1p4_2d), (alpha2p3, alpha2p3_2d)
58
59 real alpha3p4((iim + 1) * (jjm + 1))
60 real alpha3p4_2d(iim + 1, jjm + 1)
61 equivalence (alpha3p4, alpha3p4_2d)
62
63 real fext_2d(iim + 1, jjm), constang_2d(iim + 1, jjm + 1)
64 real fext((iim + 1) * jjm), constang((iim + 1) * (jjm + 1))
65 equivalence (fext, fext_2d), (constang, constang_2d)
66
67 real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
68 real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
69 equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
70
71 real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)
72 ! no dimension
73 real cvusurcu((iim + 1) * (jjm + 1)), cusurcvu((iim + 1) * (jjm + 1))
74 ! no dimension
75 equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)
76
77 real cuvscvgam1_2d(iim + 1, jjm)
78 real cuvscvgam1((iim + 1) * jjm)
79 equivalence (cuvscvgam1, cuvscvgam1_2d)
80
81 real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
82 real cuvscvgam2((iim + 1) * jjm), cvuscugam1((iim + 1) * (jjm + 1))
83 equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)
84
85 real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)
86 real cvuscugam2((iim + 1) * (jjm + 1)), cvscuvgam((iim + 1) * jjm)
87 equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)
88
89 real cuscvugam((iim + 1) * (jjm + 1))
90 real cuscvugam_2d(iim + 1, jjm + 1)
91 equivalence (cuscvugam, cuscvugam_2d)
92
93 real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2
94
95 real unsair_gam1_2d(iim + 1, jjm + 1), unsair_gam2_2d(iim + 1, jjm + 1)
96 real unsair_gam1((iim + 1) * (jjm + 1)), unsair_gam2((iim + 1) * (jjm + 1))
97 equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)
98
99 real unsairz_gam_2d(iim + 1, jjm)
100 real unsairz_gam((iim + 1) * jjm)
101 equivalence (unsairz_gam, unsairz_gam_2d)
102
103 save
104
105 contains
106
107 SUBROUTINE inigeom
108
109 ! Auteur : P. Le Van
110
111 ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
112 ! endroits que les aires aireij1_2d, ..., aireij4_2d.
113
114 ! Calcul des coefficients cu_2d, cv_2d, 1. / cu_2d**2, 1. /
115 ! cv_2d**2. Les coefficients cu_2d et cv_2d permettent de passer
116 ! des vitesses naturelles aux vitesses covariantes et
117 ! contravariantes, ou vice-versa.
118
119 ! On a :
120 ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
121 ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
122
123 ! On en tire :
124 ! u(covariant) = cu_2d * cu_2d * u(contravariant)
125 ! v(covariant) = cv_2d * cv_2d * v(contravariant)
126
127 ! x est la longitude du point en radians.
128 ! y est la latitude du point en radians.
129 !
130 ! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
131 ! cv(j) = rad * dy / dY
132 ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
133 !
134 ! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
135 ! dépendant de j uniquement, sera ici indicé aussi en i pour un
136 ! adressage plus facile en ij.
137
138 ! cv_2d est aux points v. cu_2d est aux points u. Cf. "inigeom.txt".
139
140 USE comconst, ONLY : g, omeg, rad
141 USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
142 use dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
143 yprimu1, yprimu2
144 USE paramet_m, ONLY : iip1, jjp1
145
146 ! Local:
147 INTEGER i, j
148 REAL ai14, ai23, airez, un4rad2
149 REAL coslatm, coslatp, radclatm, radclatp
150 REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
151 REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
152 REAL gamdi_gdiv, gamdi_grot, gamdi_h
153 real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
154 aireij4_2d ! in m2
155
156 !------------------------------------------------------------------
157
158 PRINT *, 'Call sequence information: inigeom'
159
160 IF (nitergdiv /= 2) THEN
161 gamdi_gdiv = coefdis / (nitergdiv - 2)
162 ELSE
163 gamdi_gdiv = 0.
164 END IF
165
166 IF (nitergrot /= 2) THEN
167 gamdi_grot = coefdis / (nitergrot - 2)
168 ELSE
169 gamdi_grot = 0.
170 END IF
171
172 IF (niterh /= 2) THEN
173 gamdi_h = coefdis / (niterh - 2)
174 ELSE
175 gamdi_h = 0.
176 END IF
177
178 print *, 'gamdi_gdiv = ', gamdi_gdiv
179 print *, "gamdi_grot = ", gamdi_grot
180 print *, "gamdi_h = ", gamdi_h
181
182 un4rad2 = 0.25 * rad * rad
183
184 ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
185 ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
186 ! chaque aire_2d(i, j), ainsi que les quatre élongations
187 ! élémentaires cuij et les quatre élongations cvij qui sont
188 ! calculées aux mêmes endroits que les aireij.
189
190 coslatm = cos(rlatu1(1))
191 radclatm = 0.5 * rad * coslatm
192
193 aireij1_2d(:iim, 1) = 0.
194 aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
195 aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
196 aireij4_2d(:iim, 1) = 0.
197
198 cuij1(:iim, 1) = 0.
199 cuij2(:iim, 1) = radclatm * xprimp025(:iim)
200 cuij3(:iim, 1) = radclatm * xprimm025(:iim)
201 cuij4(:iim, 1) = 0.
202
203 cvij1(:iim, 1) = 0.
204 cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
205 cvij3(:iim, 1) = cvij2(:iim, 1)
206 cvij4(:iim, 1) = 0.
207
208 do j = 2, jjm
209 coslatm = cos(rlatu1(j))
210 coslatp = cos(rlatu2(j-1))
211 radclatp = 0.5 * rad * coslatp
212 radclatm = 0.5 * rad * coslatm
213 ai14 = un4rad2 * coslatp * yprimu2(j-1)
214 ai23 = un4rad2 * coslatm * yprimu1(j)
215
216 aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
217 aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
218 aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
219 aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
220 cuij1(:iim, j) = radclatp * xprimp025(:iim)
221 cuij2(:iim, j) = radclatm * xprimp025(:iim)
222 cuij3(:iim, j) = radclatm * xprimm025(:iim)
223 cuij4(:iim, j) = radclatp * xprimm025(:iim)
224 cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
225 cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
226 cvij3(:iim, j) = cvij2(:iim, j)
227 cvij4(:iim, j) = cvij1(:iim, j)
228 end do
229
230 coslatp = cos(rlatu2(jjm))
231 radclatp = 0.5 * rad * coslatp
232
233 aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
234 aireij2_2d(:iim, jjp1) = 0.
235 aireij3_2d(:iim, jjp1) = 0.
236 aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
237
238 cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
239 cuij2(:iim, jjp1) = 0.
240 cuij3(:iim, jjp1) = 0.
241 cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
242
243 cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
244 cvij2(:iim, jjp1) = 0.
245 cvij3(:iim, jjp1) = 0.
246 cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
247
248 ! Périodicité :
249
250 cvij1(iip1, :) = cvij1(1, :)
251 cvij2(iip1, :) = cvij2(1, :)
252 cvij3(iip1, :) = cvij3(1, :)
253 cvij4(iip1, :) = cvij4(1, :)
254
255 cuij1(iip1, :) = cuij1(1, :)
256 cuij2(iip1, :) = cuij2(1, :)
257 cuij3(iip1, :) = cuij3(1, :)
258 cuij4(iip1, :) = cuij4(1, :)
259
260 aireij1_2d(iip1, :) = aireij1_2d(1, :)
261 aireij2_2d(iip1, :) = aireij2_2d(1, :)
262 aireij3_2d(iip1, :) = aireij3_2d(1, :)
263 aireij4_2d(iip1, :) = aireij4_2d(1, :)
264
265 DO j = 1, jjp1
266 DO i = 1, iim
267 aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
268 + aireij3_2d(i, j) + aireij4_2d(i, j)
269 alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
270 alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
271 alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
272 alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
273 alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
274 alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
275 alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
276 alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
277 END DO
278
279 aire_2d(iip1, j) = aire_2d(1, j)
280 alpha1_2d(iip1, j) = alpha1_2d(1, j)
281 alpha2_2d(iip1, j) = alpha2_2d(1, j)
282 alpha3_2d(iip1, j) = alpha3_2d(1, j)
283 alpha4_2d(iip1, j) = alpha4_2d(1, j)
284 alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
285 alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
286 alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
287 alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
288 END DO
289
290 DO j = 1, jjp1
291 DO i = 1, iim
292 aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
293 aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
294 unsaire_2d(i, j) = 1. / aire_2d(i, j)
295 unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
296 unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
297 airesurg_2d(i, j) = aire_2d(i, j) / g
298 END DO
299 aireu_2d(iip1, j) = aireu_2d(1, j)
300 unsaire_2d(iip1, j) = unsaire_2d(1, j)
301 unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
302 unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
303 airesurg_2d(iip1, j) = airesurg_2d(1, j)
304 END DO
305
306 DO j = 1, jjm
307 DO i = 1, iim
308 airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
309 aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
310 END DO
311 DO i = 1, iim
312 airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
313 + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
314 unsairez_2d(i, j) = 1. / airez
315 unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
316 fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
317 END DO
318 airev_2d(iip1, j) = airev_2d(1, j)
319 unsairez_2d(iip1, j) = unsairez_2d(1, j)
320 fext_2d(iip1, j) = fext_2d(1, j)
321 unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
322 END DO
323
324 ! Calcul des élongations cu_2d, cv_2d
325
326 DO j = 1, jjm
327 DO i = 1, iim
328 cv_2d(i, j) = 0.5 * &
329 (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
330 unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
331 END DO
332 DO i = 1, iim
333 cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
334 cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
335 cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
336 cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
337 cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
338 END DO
339 cv_2d(iip1, j) = cv_2d(1, j)
340 unscv2_2d(iip1, j) = unscv2_2d(1, j)
341 cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
342 cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
343 cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
344 cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
345 cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
346 END DO
347
348 DO j = 2, jjm
349 DO i = 1, iim
350 cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
351 + cuij3(i + 1, j))
352 unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
353 cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
354 cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
355 cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
356 cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
357 cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
358 END DO
359 cu_2d(iip1, j) = cu_2d(1, j)
360 unscu2_2d(iip1, j) = unscu2_2d(1, j)
361 cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
362 cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
363 cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
364 cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
365 cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
366 END DO
367
368 ! Calcul aux pôles
369
370 cu_2d(:, 1) = 0.
371 unscu2_2d(:, 1) = 0.
372
373 cu_2d(:, jjp1) = 0.
374 unscu2_2d(:, jjp1) = 0.
375
376 ! Calcul des aires aux pôles :
377
378 apoln = sum(aire_2d(:iim, 1))
379 apols = sum(aire_2d(:iim, jjp1))
380 unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
381 unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
382 unsapolnga2 = 1. / (apoln**(-gamdi_h))
383 unsapolsga2 = 1. / (apols**(-gamdi_h))
384
385 ! Changement F. Hourdin calcul conservatif pour fext_2d
386 ! constang_2d contient le produit a * cos (latitude) * omega
387
388 DO i = 1, iim
389 constang_2d(i, 1) = 0.
390 END DO
391 DO j = 1, jjm - 1
392 DO i = 1, iim
393 constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
394 * cos(rlatu(j + 1))
395 END DO
396 END DO
397 DO i = 1, iim
398 constang_2d(i, jjp1) = 0.
399 END DO
400
401 ! Périodicité en longitude
402 DO j = 1, jjp1
403 constang_2d(iip1, j) = constang_2d(1, j)
404 END DO
405
406 END SUBROUTINE inigeom
407
408 end module comgeom

  ViewVC Help
Powered by ViewVC 1.1.21