GCC Code Coverage Report


Directory: ./
File: phys/conlmd.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 828 0.0%
Branches: 0 770 0.0%

Line Branch Exec Source
1
2 ! $Header$
3
4 ✗ SUBROUTINE conlmd(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, &
5 ✗ ibas, itop)
6 USE dimphy
7 IMPLICIT NONE
8 ! ======================================================================
9 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
10 ! Objet: Schema de convection utilis'e dans le modele du LMD
11 ! Ajustement humide (Manabe) + Ajustement convectif (Kuo)
12 ! ======================================================================
13 include "YOMCST.h"
14 include "YOETHF.h"
15
16 ! Arguments:
17
18 REAL dtime ! pas d'integration (s)
19 REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
20 REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
21 REAL t(klon, klev) ! temperature (K)
22 REAL q(klon, klev) ! humidite specifique (kg/kg)
23 REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
24
25 REAL d_t(klon, klev) ! incrementation temperature
26 REAL d_q(klon, klev) ! incrementation humidite
27 REAL rain(klon) ! pluies (mm/s)
28 REAL snow(klon) ! neige (mm/s)
29 INTEGER ibas(klon) ! niveau du bas
30 INTEGER itop(klon) ! niveau du haut
31
32 LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
33 PARAMETER (usekuo=.TRUE.)
34
35 ✗ REAL d_t_bis(klon, klev)
36 ✗ REAL d_q_bis(klon, klev)
37 ✗ REAL rain_bis(klon)
38 ✗ REAL snow_bis(klon)
39 ✗ INTEGER ibas_bis(klon)
40 ✗ INTEGER itop_bis(klon)
41 ✗ REAL d_ql(klon, klev), d_ql_bis(klon, klev)
42 ✗ REAL rneb(klon, klev), rneb_bis(klon, klev)
43
44 INTEGER i, k
45 REAL zlvdcp, zlsdcp, zdelta, zz, za, zb
46
47 ! cc CALL fiajh ! ancienne version de Convection Manabe
48 CALL conman & ! nouvelle version de Convection
49 ! Manabe
50 ✗ (dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, ibas, itop)
51
52 IF (usekuo) THEN
53 ! cc CALL fiajc ! ancienne version de Convection Kuo
54 CALL conkuo & ! nouvelle version de Convection
55 ! Kuo
56 (dtime, paprs, pplay, t, q, conv_q, d_t_bis, d_q_bis, d_ql_bis, &
57 ✗ rneb_bis, rain_bis, snow_bis, ibas_bis, itop_bis)
58 ✗ DO k = 1, klev
59 ✗ DO i = 1, klon
60 ✗ d_t(i, k) = d_t(i, k) + d_t_bis(i, k)
61 ✗ d_q(i, k) = d_q(i, k) + d_q_bis(i, k)
62 ✗ d_ql(i, k) = d_ql(i, k) + d_ql_bis(i, k)
63 END DO
64 END DO
65 ✗ DO i = 1, klon
66 ✗ rain(i) = rain(i) + rain_bis(i)
67 ✗ snow(i) = snow(i) + snow_bis(i)
68 ✗ ibas(i) = min(ibas(i), ibas_bis(i))
69 ✗ itop(i) = max(itop(i), itop_bis(i))
70 END DO
71 END IF
72
73 ! L'eau liquide convective est dispersee dans l'air:
74
75 ✗ DO k = 1, klev
76 ✗ DO i = 1, klon
77 ✗ zlvdcp = rlvtt/rcpd/(1.0+rvtmp2*q(i,k))
78 ✗ zlsdcp = rlstt/rcpd/(1.0+rvtmp2*q(i,k))
79 ✗ zdelta = max(0., sign(1.,rtt-t(i,k)))
80 ✗ zz = d_ql(i, k) ! re-evap. de l'eau liquide
81 ✗ zb = max(0.0, zz)
82 ✗ za = -max(0.0, zz)*(zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
83 ✗ d_t(i, k) = d_t(i, k) + za
84 ✗ d_q(i, k) = d_q(i, k) + zb
85 END DO
86 END DO
87
88 ✗ RETURN
89 END SUBROUTINE conlmd
90 ✗ SUBROUTINE conman(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, &
91 ✗ snow, ibas, itop)
92 ✗ USE dimphy
93 IMPLICIT NONE
94 ! ======================================================================
95 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324
96 ! Objet: ajustement humide convectif avec la possibilite de faire
97 ! l'ajustement sur une fraction de la maille.
98 ! Methode: On impose une distribution uniforme pour la vapeur d'eau
99 ! au sein d'une maille. On applique la procedure d'ajustement
100 ! successivement a la totalite, 75%, 50%, 25% et 5% de la maille
101 ! jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente
102 ! les activites convectives et corrige le biais "trop froid et sec"
103 ! du modele.
104 ! ======================================================================
105 include "YOMCST.h"
106
107 REAL dtime ! pas d'integration (s)
108 REAL t(klon, klev) ! temperature (K)
109 REAL q(klon, klev) ! humidite specifique (kg/kg)
110 REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
111 REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
112
113 REAL d_t(klon, klev) ! incrementation temperature
114 REAL d_q(klon, klev) ! incrementation humidite
115 REAL d_ql(klon, klev) ! incrementation eau liquide
116 REAL rneb(klon, klev) ! nebulosite
117 REAL rain(klon) ! pluies (mm/s)
118 REAL snow(klon) ! neige (mm/s)
119 INTEGER ibas(klon) ! niveau du bas
120 INTEGER itop(klon) ! niveau du haut
121
122 ✗ LOGICAL afaire(klon) ! .TRUE. implique l'ajustement
123 ✗ LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif
124
125 INTEGER nb ! nombre de sous-fractions a considere
126 PARAMETER (nb=1)
127 ! cc PARAMETER (nb=3)
128
129 REAL ratqs ! largeur de la distribution pour vapeur d'eau
130 PARAMETER (ratqs=0.05)
131
132 ✗ REAL w_q(klon, klev)
133 ✗ REAL w_d_t(klon, klev), w_d_q(klon, klev), w_d_ql(klon, klev)
134 ✗ REAL w_rneb(klon, klev)
135 ✗ REAL w_rain(klon), w_snow(klon)
136 ✗ INTEGER w_ibas(klon), w_itop(klon)
137 REAL zq1, zq2
138 INTEGER i, k, n
139
140 REAL t_coup
141 PARAMETER (t_coup=234.0)
142 REAL zdp1, zdp2
143 REAL zqs1, zqs2, zdqs1, zdqs2
144 REAL zgamdz
145 REAL zflo ! flotabilite
146 REAL zsat ! sur-saturation
147 REAL zdelta, zcor, zcvm5
148 LOGICAL imprim
149
150 INTEGER ncpt
151 SAVE ncpt
152 !$OMP THREADPRIVATE(ncpt)
153 REAL frac(nb) ! valeur de la maille fractionnelle
154 SAVE frac
155 !$OMP THREADPRIVATE(frac)
156 INTEGER opt_cld(nb) ! option pour le modele nuageux
157 SAVE opt_cld
158 !$OMP THREADPRIVATE(opt_cld)
159 LOGICAL appel1er
160 SAVE appel1er
161 !$OMP THREADPRIVATE(appel1er)
162
163 ! Fonctions thermodynamiques:
164
165 include "YOETHF.h"
166 include "FCTTRE.h"
167
168 DATA frac/1.0/
169 DATA opt_cld/4/
170 ! cc DATA frac / 1.0, 0.50, 0.25/
171 ! cc DATA opt_cld / 4, 4, 4/
172
173 DATA appel1er/.TRUE./
174 DATA ncpt/0/
175
176 ✗ IF (appel1er) THEN
177 ✗ PRINT *, 'conman, nb:', nb
178 ✗ PRINT *, 'conman, frac:', frac
179 ✗ PRINT *, 'conman, opt_cld:', opt_cld
180 ✗ appel1er = .FALSE.
181 END IF
182
183 ! Initialiser les sorties a zero:
184
185 ✗ DO k = 1, klev
186 ✗ DO i = 1, klon
187 ✗ d_t(i, k) = 0.0
188 ✗ d_q(i, k) = 0.0
189 ✗ d_ql(i, k) = 0.0
190 ✗ rneb(i, k) = 0.0
191 END DO
192 END DO
193 ✗ DO i = 1, klon
194 ✗ ibas(i) = klev
195 ✗ itop(i) = 1
196 ✗ rain(i) = 0.0
197 ✗ snow(i) = 0.0
198 END DO
199
200 ! S'il n'y a pas d'instabilite conditionnelle,
201 ! pas la penne de se fatiguer:
202
203 ✗ DO i = 1, klon
204 ✗ afaire(i) = .FALSE.
205 END DO
206 ✗ DO k = 1, klev - 1
207 ✗ DO i = 1, klon
208 IF (thermcep) THEN
209 ✗ zdelta = max(0., sign(1.,rtt-t(i,k)))
210 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
211 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
212 ✗ zqs1 = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
213 ✗ zqs1 = min(0.5, zqs1)
214 ✗ zcor = 1./(1.-retv*zqs1)
215 ✗ zqs1 = zqs1*zcor
216 ✗ zdqs1 = foede(t(i,k), zdelta, zcvm5, zqs1, zcor)
217
218 ✗ zdelta = max(0., sign(1.,rtt-t(i,k+1)))
219 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
220 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k+1))
221 ✗ zqs2 = r2es*foeew(t(i,k+1), zdelta)/pplay(i, k+1)
222 ✗ zqs2 = min(0.5, zqs2)
223 ✗ zcor = 1./(1.-retv*zqs2)
224 ✗ zqs2 = zqs2*zcor
225 ✗ zdqs2 = foede(t(i,k+1), zdelta, zcvm5, zqs2, zcor)
226 ELSE
227 IF (t(i,k)<t_coup) THEN
228 zqs1 = qsats(t(i,k))/pplay(i, k)
229 zdqs1 = dqsats(t(i,k), zqs1)
230
231 zqs2 = qsats(t(i,k+1))/pplay(i, k+1)
232 zdqs2 = dqsats(t(i,k+1), zqs2)
233 ELSE
234 zqs1 = qsatl(t(i,k))/pplay(i, k)
235 zdqs1 = dqsatl(t(i,k), zqs1)
236
237 zqs2 = qsatl(t(i,k+1))/pplay(i, k+1)
238 zdqs2 = dqsatl(t(i,k+1), zqs2)
239 END IF
240 END IF
241 ✗ zdp1 = paprs(i, k) - paprs(i, k+1)
242 ✗ zdp2 = paprs(i, k+1) - paprs(i, k+2)
243 zgamdz = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd*(rd*(t(i, &
244 k)*zdp1+t(i,k+1)*zdp2)/(zdp1+zdp2)+rlvtt*(zqs1*zdp1+zqs2*zdp2)/(zdp1+ &
245 ✗ zdp2))/(1.0+(zdqs1*zdp1+zdqs2*zdp2)/(zdp1+zdp2))
246 ✗ zflo = t(i, k) + zgamdz - t(i, k+1)
247 zsat = (q(i,k)-zqs1)*zdp1 + (q(i,k+1)-zqs2)*zdp2
248 ✗ IF (zflo>0.0) afaire(i) = .TRUE.
249 ! erreur IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE.
250 END DO
251 END DO
252
253 ✗ imprim = mod(ncpt, 48) == 0
254 ✗ DO n = 1, nb
255
256 ✗ DO k = 1, klev
257 ✗ DO i = 1, klon
258 ✗ IF (afaire(i)) THEN
259 ✗ zq1 = q(i, k)*(1.0-ratqs)
260 ✗ zq2 = q(i, k)*(1.0+ratqs)
261 ✗ w_q(i, k) = zq2 - frac(n)/2.0*(zq2-zq1)
262 END IF
263 END DO
264 END DO
265
266 CALL conmanv(dtime, paprs, pplay, t, w_q, afaire, opt_cld(n), w_d_t, &
267 w_d_q, w_d_ql, w_rneb, w_rain, w_snow, w_ibas, w_itop, accompli, &
268 ✗ imprim)
269 ✗ DO k = 1, klev
270 ✗ DO i = 1, klon
271 ✗ IF (afaire(i) .AND. accompli(i)) THEN
272 ✗ d_t(i, k) = w_d_t(i, k)*frac(n)
273 ✗ d_q(i, k) = w_d_q(i, k)*frac(n)
274 ✗ d_ql(i, k) = w_d_ql(i, k)*frac(n)
275 ✗ IF (nint(w_rneb(i,k))==1) rneb(i, k) = frac(n)
276 END IF
277 END DO
278 END DO
279 ✗ DO i = 1, klon
280 ✗ IF (afaire(i) .AND. accompli(i)) THEN
281 ✗ rain(i) = w_rain(i)*frac(n)
282 ✗ snow(i) = w_snow(i)*frac(n)
283 ✗ ibas(i) = min(ibas(i), w_ibas(i))
284 ✗ itop(i) = max(itop(i), w_itop(i))
285 END IF
286 END DO
287 ✗ DO i = 1, klon
288 ✗ IF (afaire(i) .AND. accompli(i)) afaire(i) = .FALSE.
289 END DO
290
291 END DO
292
293 ✗ ncpt = ncpt + 1
294
295 ✗ RETURN
296 END SUBROUTINE conman
297 ✗ SUBROUTINE conmanv(dtime, paprs, pplay, t, q, afaire, opt_cld, d_t, d_q, &
298 ✗ d_ql, rneb, rain, snow, ibas, itop, accompli, imprim)
299 USE dimphy
300 IMPLICIT NONE
301 ! ======================================================================
302 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
303 ! Objet: ajustement humide (convection proposee par Manabe).
304 ! Pour une colonne verticale, il peut avoir plusieurs blocs
305 ! necessitant l'ajustement. ibas est le bas du plus bas bloc
306 ! et itop est le haut du plus haut bloc
307 ! ======================================================================
308 include "YOMCST.h"
309
310 ! Arguments:
311
312 REAL dtime ! pas d'integration (s)
313 REAL t(klon, klev) ! temperature (K)
314 REAL q(klon, klev) ! humidite specifique (kg/kg)
315 REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
316 REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
317 INTEGER opt_cld ! comment traiter l'eau liquide
318 LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input)
319 LOGICAL imprim ! .T. pour imprimer quelques diagnostiques
320
321 REAL d_t(klon, klev) ! incrementation temperature
322 REAL d_q(klon, klev) ! incrementation humidite
323 REAL d_ql(klon, klev) ! incrementation eau liquide
324 REAL rneb(klon, klev) ! nebulosite
325 REAL rain(klon) ! pluies (mm/s)
326 REAL snow(klon) ! neige (mm/s)
327 INTEGER ibas(klon) ! niveau du bas
328 INTEGER itop(klon) ! niveau du haut
329 LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output)
330
331 ! Quelques options:
332
333 LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait
334 PARAMETER (new_top=.FALSE.)
335 LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
336 PARAMETER (evap_prec=.TRUE.)
337 REAL coef_eva
338 PARAMETER (coef_eva=1.0E-05)
339 REAL t_coup
340 PARAMETER (t_coup=234.0)
341 REAL seuil_vap
342 PARAMETER (seuil_vap=1.0E-10)
343 LOGICAL old_tau ! implique precip nulle, si vrai.
344 PARAMETER (old_tau=.FALSE.)
345 ✗ REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
346 REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
347 PARAMETER (dpmin=0.15, tomax=0.97)
348 REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
349 PARAMETER (dpmax=0.30, tomin=0.05)
350 REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
351 PARAMETER (deep_sig=0.50, deep_to=0.05)
352 LOGICAL exigent ! implique un calcul supplementaire pour Qs
353 PARAMETER (exigent=.FALSE.)
354
355 INTEGER kbase
356 PARAMETER (kbase=0)
357
358 ! Variables locales:
359
360 INTEGER nexpo
361 INTEGER i, k, k1min, k1max, k2min, k2max, is
362 ✗ REAL zgamdz(klon, klev-1)
363 ✗ REAL zt(klon, klev), zq(klon, klev)
364 ✗ REAL zqs(klon, klev), zdqs(klon, klev)
365 ✗ REAL zqmqsdp(klon, klev)
366 ✗ REAL ztnew(klon, klev), zqnew(klon, klev)
367 ✗ REAL zcond(klon), zvapo(klon), zrapp(klon)
368 ✗ REAL zrfl(klon), zrfln, zqev, zqevt
369 ✗ REAL zsat(klon) ! sur-saturation
370 ✗ REAL zflo(klon) ! flotabilite
371 ✗ REAL za(klon), zb(klon), zc(klon)
372 ✗ INTEGER k1(klon), k2(klon)
373 REAL zdelta, zcor, zcvm5
374 ✗ REAL delp(klon, klev)
375 ✗ LOGICAL possible(klon), todo(klon), etendre(klon)
376 ✗ LOGICAL aller(klon), todobis(klon)
377 REAL zalfa
378 INTEGER nbtodo, nbdone
379
380 ! Fonctions thermodynamiques:
381
382 include "YOETHF.h"
383 include "FCTTRE.h"
384
385 ✗ DO k = 1, klev
386 ✗ DO i = 1, klon
387 ✗ delp(i, k) = paprs(i, k) - paprs(i, k+1)
388 END DO
389 END DO
390
391 ! Initialiser les sorties a zero
392
393 ✗ DO k = 1, klev
394 ✗ DO i = 1, klon
395 ✗ d_t(i, k) = 0.0
396 ✗ d_q(i, k) = 0.0
397 ✗ d_ql(i, k) = 0.0
398 ✗ rneb(i, k) = 0.0
399 END DO
400 END DO
401 ✗ DO i = 1, klon
402 ✗ ibas(i) = klev
403 ✗ itop(i) = 1
404 ✗ rain(i) = 0.0
405 ✗ snow(i) = 0.0
406 ✗ accompli(i) = .FALSE.
407 END DO
408
409 ! Preparations
410
411 ✗ DO k = 1, klev
412 ✗ DO i = 1, klon
413 ✗ IF (afaire(i)) THEN
414 ✗ zt(i, k) = t(i, k)
415 ✗ zq(i, k) = q(i, k)
416
417 ! Calculer Qs et L/Cp*dQs/dT
418
419 IF (thermcep) THEN
420 ✗ zdelta = max(0., sign(1.,rtt-zt(i,k)))
421 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
422 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i,k))
423 ✗ zqs(i, k) = r2es*foeew(zt(i,k), zdelta)/pplay(i, k)
424 ✗ zqs(i, k) = min(0.5, zqs(i,k))
425 ✗ zcor = 1./(1.-retv*zqs(i,k))
426 ✗ zqs(i, k) = zqs(i, k)*zcor
427 ✗ zdqs(i, k) = foede(zt(i,k), zdelta, zcvm5, zqs(i,k), zcor)
428 ELSE
429 IF (zt(i,k)<t_coup) THEN
430 zqs(i, k) = qsats(zt(i,k))/pplay(i, k)
431 zdqs(i, k) = dqsats(zt(i,k), zqs(i,k))
432 ELSE
433 zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
434 zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
435 END IF
436 END IF
437
438 ! Calculer (q-qs)*dp
439 ✗ zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
440 END IF
441 END DO
442 END DO
443
444 ! -----zgama is the moist convective lapse rate (-dT/dz).
445 ! -----zgamdz(*,k) est la difference minimale autorisee des temperatures
446 ! -----entre deux couches (k et k+1), c.a.d. si T(k+1)-T(k) est inferieur
447 ! -----a zgamdz(*,k), alors ces 2 couches sont instables conditionnellement
448
449 ✗ DO k = 1, klev - 1
450 ✗ DO i = 1, klon
451 ✗ IF (afaire(i)) THEN
452 zgamdz(i, k) = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd*(rd*(zt( &
453 i,k)*delp(i,k)+zt(i,k+1)*delp(i,k+1))/(delp(i,k)+delp(i, &
454 k+1))+rlvtt*(zqs(i,k)*delp(i,k)+zqs(i,k+1)*delp(i,k+1))/(delp(i, &
455 k)+delp(i,k+1)))/(1.0+(zdqs(i,k)*delp(i,k)+zdqs(i,k+1)*delp(i, &
456 ✗ k+1))/(delp(i,k)+delp(i,k+1)))
457 END IF
458 END DO
459 END DO
460
461 ! On cherche la presence simultanee d'instabilite conditionnelle
462 ! et de sur-saturation. Sinon, pas la penne de se fatiguer:
463
464 ✗ DO i = 1, klon
465 ✗ possible(i) = .FALSE.
466 END DO
467 ✗ DO k = 2, klev
468 ✗ DO i = 1, klon
469 ✗ IF (afaire(i)) THEN
470 ✗ zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
471 ✗ zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
472 ✗ IF (zflo(i)>0.0 .AND. zsat(i)>0.0) possible(i) = .TRUE.
473 END IF
474 END DO
475 END DO
476
477 ✗ DO i = 1, klon
478 ✗ IF (possible(i)) THEN
479 ✗ k1(i) = kbase
480 ✗ k2(i) = k1(i) + 1
481 END IF
482 END DO
483
484 810 CONTINUE ! chercher le bas de la colonne a ajuster
485
486 k2min = klev
487 ✗ DO i = 1, klon
488 ✗ todo(i) = .FALSE.
489 ✗ aller(i) = .TRUE.
490 ✗ IF (possible(i)) k2min = min(k2min, k2(i))
491 END DO
492 ✗ IF (k2min==klev) GO TO 860
493 ✗ DO k = k2min, klev - 1
494 ✗ DO i = 1, klon
495 ✗ IF (possible(i) .AND. k>=k2(i) .AND. aller(i)) THEN
496 ✗ zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
497 ✗ zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k+1)
498 ✗ IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
499 ✗ k1(i) = k
500 ✗ k2(i) = k + 1
501 ✗ todo(i) = .TRUE.
502 ✗ aller(i) = .FALSE.
503 END IF
504 END IF
505 END DO
506 END DO
507 ✗ DO i = 1, klon
508 ✗ IF (possible(i) .AND. aller(i)) THEN
509 ✗ todo(i) = .FALSE.
510 ✗ k1(i) = klev
511 ✗ k2(i) = klev
512 END IF
513 END DO
514
515 ! CC DO i = 1, klon
516 ! CC IF (possible(i)) THEN
517 ! CC 811 k2(i) = k2(i) + 1
518 ! CC IF (k2(i) .GT. klev) THEN
519 ! CC todo(i) = .FALSE.
520 ! CC GOTO 812
521 ! CC ENDIF
522 ! CC k = k2(i)
523 ! CC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
524 ! CC zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1)
525 ! CC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 811
526 ! CC k1(i) = k2(i) - 1
527 ! CC todo(i) = .TRUE.
528 ! CC ENDIF
529 ! CC 812 CONTINUE
530 ! CC ENDDO
531
532 820 CONTINUE ! chercher le haut de la colonne
533
534 k2min = klev
535 ✗ DO i = 1, klon
536 ✗ aller(i) = .TRUE.
537 ✗ IF (todo(i)) k2min = min(k2min, k2(i))
538 END DO
539 ✗ IF (k2min<klev) THEN
540 ✗ DO k = k2min, klev
541 ✗ DO i = 1, klon
542 ✗ IF (todo(i) .AND. k>k2(i) .AND. aller(i)) THEN
543 ✗ zsat(i) = zsat(i) + zqmqsdp(i, k)
544 ✗ zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
545 ✗ IF (zflo(i)<=0.0 .OR. zsat(i)<=0.0) THEN
546 ✗ aller(i) = .FALSE.
547 ELSE
548 ✗ k2(i) = k
549 END IF
550 END IF
551 END DO
552 END DO
553 ! error is = 0
554 ! error DO i = 1, klon
555 ! error IF(todo(i).AND.aller(i)) THEN
556 ! error is = is + 1
557 ! error todo(i) = .FALSE.
558 ! error k2(i) = klev
559 ! error ENDIF
560 ! error ENDDO
561 ! error IF (is.GT.0) THEN
562 ! error PRINT*, "Bizard. je pourrais continuer mais j arrete"
563 ! error CALL abort
564 ! error ENDIF
565 END IF
566
567 ! CC DO i = 1, klon
568 ! CC IF (todo(i)) THEN
569 ! CC 821 CONTINUE
570 ! CC IF (k2(i) .EQ. klev) GOTO 822
571 ! CC k = k2(i) + 1
572 ! CC zsat(i) = zsat(i) + zqmqsdp(i,k)
573 ! CC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
574 ! CC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 822
575 ! CC k2(i) = k
576 ! CC GOTO 821
577 ! CC ENDIF
578 ! CC 822 CONTINUE
579 ! CC ENDDO
580
581 830 CONTINUE ! faire l'ajustement en sachant k1 et k2
582
583 ✗ is = 0
584 ✗ DO i = 1, klon
585 ✗ IF (todo(i)) THEN
586 ✗ IF (k2(i)<=k1(i)) is = is + 1
587 END IF
588 END DO
589 ✗ IF (is>0) THEN
590 ✗ PRINT *, 'Impossible: k1 trop grand ou k2 trop petit'
591 ✗ PRINT *, 'is=', is
592 ✗ CALL abort
593 END IF
594
595 k1min = klev
596 k1max = 1
597 k2max = 1
598 ✗ DO i = 1, klon
599 ✗ IF (todo(i)) THEN
600 ✗ k1min = min(k1min, k1(i))
601 ✗ k1max = max(k1max, k1(i))
602 ✗ k2max = max(k2max, k2(i))
603 END IF
604 END DO
605
606 ✗ DO i = 1, klon
607 ✗ IF (todo(i)) THEN
608 ✗ k = k1(i)
609 ✗ za(i) = 0.
610 zb(i) = (rcpd*(1.+zdqs(i,k))*(zt(i,k)-za(i))-rlvtt*(zqs(i,k)-zq(i, &
611 ✗ k)))*delp(i, k)
612 ✗ zc(i) = delp(i, k)*rcpd*(1.+zdqs(i,k))
613 END IF
614 END DO
615
616 ✗ DO k = k1min, k2max
617 ✗ DO i = 1, klon
618 ✗ IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) THEN
619 ✗ za(i) = za(i) + zgamdz(i, k-1)
620 zb(i) = zb(i) + (rcpd*(1.+zdqs(i,k))*(zt(i,k)-za(i))-rlvtt*(zqs(i, &
621 ✗ k)-zq(i,k)))*delp(i, k)
622 ✗ zc(i) = zc(i) + delp(i, k)*rcpd*(1.+zdqs(i,k))
623 END IF
624 END DO
625 END DO
626
627 ✗ DO i = 1, klon
628 ✗ IF (todo(i)) THEN
629 ✗ k = k1(i)
630 ✗ ztnew(i, k) = zb(i)/zc(i)
631 ✗ zqnew(i, k) = zqs(i, k) + (ztnew(i,k)-zt(i,k))*rcpd/rlvtt*zdqs(i, k)
632 END IF
633 END DO
634
635 ✗ DO k = k1min, k2max
636 ✗ DO i = 1, klon
637 ✗ IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) THEN
638 ✗ ztnew(i, k) = ztnew(i, k-1) + zgamdz(i, k-1)
639 ✗ zqnew(i, k) = zqs(i, k) + (ztnew(i,k)-zt(i,k))*rcpd/rlvtt*zdqs(i, k)
640 END IF
641 END DO
642 END DO
643
644 ! Quantite de condensation produite pendant l'ajustement:
645
646 ✗ DO i = 1, klon
647 ✗ zcond(i) = 0.0
648 END DO
649 ✗ DO k = k1min, k2max
650 ✗ DO i = 1, klon
651 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
652 ✗ rneb(i, k) = 1.0
653 ✗ zcond(i) = zcond(i) + (zq(i,k)-zqnew(i,k))*delp(i, k)/rg
654 END IF
655 END DO
656 END DO
657
658 ! Si condensation negative, effort completement perdu:
659
660 ✗ DO i = 1, klon
661 ✗ IF (todo(i) .AND. zcond(i)<=0.) todo(i) = .FALSE.
662 END DO
663
664 ! L'ajustement a ete accompli, meme les calculs accessoires
665 ! ne sont pas encore faits:
666
667 ✗ DO i = 1, klon
668 ✗ IF (todo(i)) accompli(i) = .TRUE.
669 END DO
670
671 ! =====
672 ! Une fois que la condensation a lieu, on doit construire un
673 ! "modele nuageux" pour partager la condensation entre l'eau
674 ! liquide nuageuse et la precipitation (leur rapport toliq
675 ! est calcule selon l'epaisseur nuageuse). Je suppose que
676 ! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
677 ! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
678 ! lineaire entre dpmin et dpmax).
679 ! =====
680 ✗ DO i = 1, klon
681 ✗ IF (todo(i)) THEN
682 toliq(i) = tomax - ((paprs(i,k1(i))-paprs(i,k2(i)+1))/paprs(i,1)-dpmin) &
683 ✗ *(tomax-tomin)/(dpmax-dpmin)
684 ✗ toliq(i) = max(tomin, min(tomax,toliq(i)))
685 ✗ IF (pplay(i,k2(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
686 IF (old_tau) toliq(i) = 1.0
687 END IF
688 END DO
689 ! =====
690 ! On doit aussi determiner la distribution verticale de
691 ! l'eau nuageuse. Plusieurs options sont proposees:
692
693 ! (0) La condensation precipite integralement (toliq ne sera
694 ! pas utilise).
695 ! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
696 ! a la vapeur d'eau locale.
697 ! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
698 ! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
699 ! est effectivement diminuee pendant le processus d'ajustement.
700 ! (4) Elle est en fonction (lineaire ou exponentielle) de la
701 ! distance (epaisseur en pression) avec le niveau k1 (la couche
702 ! k1 n'aura donc pas d'eau liquide).
703 ! =====
704
705 ✗ IF (opt_cld==0) THEN
706
707 ✗ DO i = 1, klon
708 ✗ IF (todo(i)) zrfl(i) = zcond(i)/dtime
709 END DO
710
711 ✗ ELSE IF (opt_cld==1) THEN
712
713 ✗ DO i = 1, klon
714 ✗ IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
715 END DO
716 ✗ DO k = k1min, k2max
717 ✗ DO i = 1, klon
718 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
719 ✗ zqnew(i, k)*delp(i, k)/rg
720 END DO
721 END DO
722 ✗ DO i = 1, klon
723 ✗ IF (todo(i)) THEN
724 ✗ zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
725 ✗ zrapp(i) = max(0., min(1.,zrapp(i)))
726 ✗ zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
727 END IF
728 END DO
729 ✗ DO k = k1min, k2max
730 ✗ DO i = 1, klon
731 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
732 ✗ d_ql(i, k) = d_ql(i, k) + zrapp(i)*zqnew(i, k)
733 END IF
734 END DO
735 END DO
736
737 ✗ ELSE IF (opt_cld==2) THEN
738
739 ✗ DO i = 1, klon
740 ✗ IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
741 END DO
742 ✗ DO k = k1min, k2max
743 ✗ DO i = 1, klon
744 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
745 ✗ delp(i, k)/rg
746 END DO
747 END DO
748 ✗ DO k = k1min, k2max
749 ✗ DO i = 1, klon
750 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
751 ✗ d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)
752 END IF
753 END DO
754 END DO
755 ✗ DO i = 1, klon
756 ✗ IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
757 END DO
758
759 ✗ ELSE IF (opt_cld==3) THEN
760
761 ✗ DO i = 1, klon
762 ✗ IF (todo(i)) zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
763 END DO
764 ✗ DO k = k1min, k2max
765 ✗ DO i = 1, klon
766 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
767 ✗ max(0.0, zq(i,k)-zqnew(i,k))*delp(i, k)/rg
768 END DO
769 END DO
770 ✗ DO k = k1min, k2max
771 ✗ DO i = 1, klon
772 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i) .AND. zvapo(i)>0.0) d_ql(i, &
773 k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*max(0.0, zq(i,k)-zqnew &
774 ✗ (i,k))
775 END DO
776 END DO
777 ✗ DO i = 1, klon
778 ✗ IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
779 END DO
780
781 ✗ ELSE IF (opt_cld==4) THEN
782
783 nexpo = 3
784 ! cc nexpo = 1 ! distribution lineaire
785
786 ✗ DO i = 1, klon
787 ✗ IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
788 END DO ! (avec ponderation)
789 ✗ DO k = k1min, k2max
790 ✗ DO i = 1, klon
791 ✗ IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
792 ✗ delp(i, k)/rg*(pplay(i,k1(i))-pplay(i,k))**nexpo
793 END DO
794 END DO
795 ✗ DO k = k1min, k2max
796 ✗ DO i = 1, klon
797 ✗ IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) d_ql(i, k) = d_ql(i, &
798 ✗ k) + toliq(i)*zcond(i)/zvapo(i)*(pplay(i,k1(i))-pplay(i,k))**nexpo
799 END DO
800 END DO
801 ✗ DO i = 1, klon
802 ✗ IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
803 END DO
804
805 ELSE ! valeur non-prevue pour opt_cld
806
807 ✗ PRINT *, 'opt_cld est faux:', opt_cld
808 ✗ CALL abort
809
810 END IF ! fin de opt_cld
811
812 ! L'eau precipitante peut etre evaporee:
813
814 zalfa = 0.05
815 ✗ IF (evap_prec .AND. (k1max>=2)) THEN
816 ✗ DO k = k1max - 1, 1, -1
817 ✗ DO i = 1, klon
818 ✗ IF (todo(i) .AND. k<k1(i) .AND. zrfl(i)>0.0) THEN
819 ✗ zqev = max(0.0, (zqs(i,k)-zq(i,k))*zalfa)
820 zqevt = coef_eva*(1.0-zq(i,k)/zqs(i,k))*sqrt(zrfl(i))*delp(i, k)/ &
821 ✗ pplay(i, k)*zt(i, k)*rd/rg
822 ✗ zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/delp(i, k)
823 ✗ zqev = min(zqev, zqevt)
824 ✗ zrfln = zrfl(i) - zqev*(delp(i,k))/rg/dtime
825 ✗ zq(i, k) = zq(i, k) - (zrfln-zrfl(i))*(rg/(delp(i,k)))*dtime
826 zt(i, k) = zt(i, k) + (zrfln-zrfl(i))*(rg/(delp(i, &
827 ✗ k)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i,k))
828 ✗ zrfl(i) = zrfln
829 END IF
830 END DO
831 END DO
832 END IF
833
834 ! La temperature de la premiere couche determine la pluie ou la neige:
835
836 ✗ DO i = 1, klon
837 ✗ IF (todo(i)) THEN
838 ✗ IF (zt(i,1)>rtt) THEN
839 ✗ rain(i) = rain(i) + zrfl(i)
840 ELSE
841 ✗ snow(i) = snow(i) + zrfl(i)
842 END IF
843 END IF
844 END DO
845
846 ! Mise a jour de la temperature et de l'humidite
847
848 ✗ DO k = k1min, k2max
849 ✗ DO i = 1, klon
850 ✗ IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
851 ✗ zt(i, k) = ztnew(i, k)
852 ✗ zq(i, k) = zqnew(i, k)
853 END IF
854 END DO
855 END DO
856
857 ! Re-calculer certaines variables pour etendre et re-ajuster la colonne
858
859 IF (exigent) THEN
860 DO k = 1, klev
861 DO i = 1, klon
862 IF (todo(i)) THEN
863 IF (thermcep) THEN
864 zdelta = max(0., sign(1.,rtt-zt(i,k)))
865 zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
866 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i,k))
867 zqs(i, k) = r2es*foeew(zt(i,k), zdelta)/pplay(i, k)
868 zqs(i, k) = min(0.5, zqs(i,k))
869 zcor = 1./(1.-retv*zqs(i,k))
870 zqs(i, k) = zqs(i, k)*zcor
871 zdqs(i, k) = foede(zt(i,k), zdelta, zcvm5, zqs(i,k), zcor)
872 ELSE
873 IF (zt(i,k)<t_coup) THEN
874 zqs(i, k) = qsats(zt(i,k))/pplay(i, k)
875 zdqs(i, k) = dqsats(zt(i,k), zqs(i,k))
876 ELSE
877 zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
878 zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
879 END IF
880 END IF
881 END IF
882 END DO
883 END DO
884 END IF
885
886 IF (exigent) THEN
887 DO k = 1, klev - 1
888 DO i = 1, klon
889 IF (todo(i)) THEN
890 zgamdz(i, k) = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd*(rd*( &
891 zt(i,k)*delp(i,k)+zt(i,k+1)*delp(i,k+1))/(delp(i,k)+delp(i, &
892 k+1))+rlvtt*(zqs(i,k)*delp(i,k)+zqs(i,k+1)*delp(i,k+1))/(delp(i, &
893 k)+delp(i,k+1)))/(1.0+(zdqs(i,k)*delp(i,k)+zdqs(i,k+1)*delp(i, &
894 k+1))/(delp(i,k)+delp(i,k+1)))
895 END IF
896 END DO
897 END DO
898 END IF
899
900 ! Puisque l'humidite a ete modifiee, on re-fait (q-qs)*dp
901
902 ✗ DO k = 1, klev
903 ✗ DO i = 1, klon
904 ✗ IF (todo(i)) THEN
905 ✗ zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
906 END IF
907 END DO
908 END DO
909
910 ! Verifier si l'on peut etendre le bas de la colonne
911
912 ✗ DO i = 1, klon
913 ✗ etendre(i) = .FALSE.
914 END DO
915
916 k1max = 1
917 ✗ DO i = 1, klon
918 ✗ IF (todo(i) .AND. k1(i)>(kbase+1)) THEN
919 k = k1(i)
920 ✗ zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
921 ✗ zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
922 ! sc voici l'ancienne ligne:
923 ! sc IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN
924 ! sc sylvain: il faut RESPECTER les 2 criteres:
925 ✗ IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
926 ✗ etendre(i) = .TRUE.
927 ✗ k1(i) = k1(i) - 1
928 ✗ k1max = max(k1max, k1(i))
929 ✗ aller(i) = .TRUE.
930 END IF
931 END IF
932 END DO
933
934 ✗ IF (k1max>(kbase+1)) THEN
935 ✗ DO k = k1max, kbase + 1, -1
936 ✗ DO i = 1, klon
937 ✗ IF (etendre(i) .AND. k<k1(i) .AND. aller(i)) THEN
938 ✗ zsat(i) = zsat(i) + zqmqsdp(i, k)
939 ✗ zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
940 ✗ IF (zsat(i)<=0.0 .OR. zflo(i)<=0.0) THEN
941 ✗ aller(i) = .FALSE.
942 ELSE
943 ✗ k1(i) = k
944 END IF
945 END IF
946 END DO
947 END DO
948 ✗ DO i = 1, klon
949 ✗ IF (etendre(i) .AND. aller(i)) THEN
950 ✗ k1(i) = 1
951 END IF
952 END DO
953 END IF
954
955 ! CC DO i = 1, klon
956 ! CC IF (etendre(i)) THEN
957 ! CC 840 k = k1(i)
958 ! CC IF (k.GT.1) THEN
959 ! CC zsat(i) = zsat(i) + zqmqsdp(i,k-1)
960 ! CC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
961 ! CC IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN
962 ! CC k1(i) = k - 1
963 ! CC GOTO 840
964 ! CC ENDIF
965 ! CC ENDIF
966 ! CC ENDIF
967 ! CC ENDDO
968
969 ✗ DO i = 1, klon
970 ✗ todobis(i) = todo(i)
971 ✗ todo(i) = .FALSE.
972 END DO
973 ✗ is = 0
974 ✗ DO i = 1, klon
975 ✗ IF (etendre(i)) THEN
976 ✗ todo(i) = .TRUE.
977 ✗ is = is + 1
978 END IF
979 END DO
980 ✗ IF (is>0) THEN
981 IF (new_top) THEN
982 GO TO 820 ! chercher de nouveau le sommet k2
983 ELSE
984 GO TO 830 ! supposer que le sommet est celui deja trouve
985 END IF
986 END IF
987
988 ✗ DO i = 1, klon
989 ✗ possible(i) = .FALSE.
990 END DO
991 ✗ is = 0
992 ✗ DO i = 1, klon
993 ✗ IF (todobis(i) .AND. k2(i)<klev) THEN
994 ✗ is = is + 1
995 ✗ possible(i) = .TRUE.
996 END IF
997 END DO
998 ✗ IF (is>0) GO TO 810 !on cherche en haut d'autres blocks
999 ! a ajuster a partir du sommet de la colonne precedente
1000
1001 860 CONTINUE ! Calculer les tendances et diagnostiques
1002 ! cc print*, "Apres 860"
1003
1004 ✗ DO k = 1, klev
1005 ✗ DO i = 1, klon
1006 ✗ IF (accompli(i)) THEN
1007 ✗ d_t(i, k) = zt(i, k) - t(i, k)
1008 ✗ zq(i, k) = max(zq(i,k), seuil_vap)
1009 ✗ d_q(i, k) = zq(i, k) - q(i, k)
1010 END IF
1011 END DO
1012 END DO
1013
1014 ✗ DO i = 1, klon
1015 ✗ IF (accompli(i)) THEN
1016 ✗ DO k = 1, klev
1017 ✗ IF (rneb(i,k)>0.0) THEN
1018 ✗ ibas(i) = k
1019 ✗ GO TO 807
1020 END IF
1021 END DO
1022 807 CONTINUE
1023 ✗ DO k = klev, 1, -1
1024 ✗ IF (rneb(i,k)>0.0) THEN
1025 ✗ itop(i) = k
1026 ✗ GO TO 808
1027 END IF
1028 END DO
1029 808 CONTINUE
1030 END IF
1031 END DO
1032
1033 ✗ IF (imprim) THEN
1034 ✗ nbtodo = 0
1035 ✗ nbdone = 0
1036 ✗ DO i = 1, klon
1037 ✗ IF (afaire(i)) nbtodo = nbtodo + 1
1038 ✗ IF (accompli(i)) nbdone = nbdone + 1
1039 END DO
1040 ✗ PRINT *, 'nbTodo, nbDone=', nbtodo, nbdone
1041 END IF
1042
1043 ✗ RETURN
1044 END SUBROUTINE conmanv
1045 ✗ SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
1046 ✗ rain, snow, ibas, itop)
1047 ✗ USE dimphy
1048 IMPLICIT NONE
1049 ! ======================================================================
1050 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1051 ! Objet: Schema de convection de type Kuo (1965).
1052 ! Cette version du code peut calculer le niveau de depart
1053 ! N.B. version vectorielle (le 6 oct. 1997)
1054 ! ======================================================================
1055 include "YOMCST.h"
1056
1057 ! Arguments:
1058
1059 REAL dtime ! intervalle du temps (s)
1060 REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
1061 REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
1062 REAL t(klon, klev) ! temperature (K)
1063 REAL q(klon, klev) ! humidite specifique
1064 REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
1065
1066 REAL d_t(klon, klev) ! incrementation temperature
1067 REAL d_q(klon, klev) ! incrementation humidite
1068 REAL d_ql(klon, klev) ! incrementation eau liquide
1069 REAL rneb(klon, klev) ! nebulosite
1070 REAL rain(klon) ! pluies (mm/s)
1071 REAL snow(klon) ! neige (mm/s)
1072 INTEGER itop(klon) ! niveau du sommet
1073 INTEGER ibas(klon) ! niveau du bas
1074
1075 ✗ LOGICAL ldcum(klon) ! convection existe
1076 ✗ LOGICAL todo(klon)
1077
1078 ! Quelsques options:
1079
1080 LOGICAL calcfcl ! calculer le niveau de convection libre
1081 PARAMETER (calcfcl=.TRUE.)
1082 INTEGER ldepar ! niveau fixe de convection libre
1083 PARAMETER (ldepar=4)
1084 INTEGER opt_cld ! comment traiter l'eau liquide
1085 PARAMETER (opt_cld=4) ! valeur possible: 0, 1, 2, 3 ou 4
1086 LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
1087 PARAMETER (evap_prec=.TRUE.)
1088 REAL coef_eva
1089 PARAMETER (coef_eva=1.0E-05)
1090 LOGICAL new_deh ! nouvelle facon de calculer dH
1091 PARAMETER (new_deh=.FALSE.)
1092 REAL t_coup
1093 PARAMETER (t_coup=234.0)
1094 LOGICAL old_tau ! implique precipitation nulle
1095 PARAMETER (old_tau=.FALSE.)
1096 ✗ REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
1097 REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
1098 PARAMETER (dpmin=0.15, tomax=0.97)
1099 REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
1100 PARAMETER (dpmax=0.30, tomin=0.05)
1101 REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
1102 PARAMETER (deep_sig=0.50, deep_to=0.05)
1103
1104 ! Variables locales:
1105
1106 INTEGER nexpo
1107 ✗ LOGICAL nuage(klon)
1108 INTEGER i, k, kbmin, kbmax, khmax
1109 ✗ REAL ztotal(klon, klev), zdeh(klon, klev)
1110 ✗ REAL zgz(klon, klev)
1111 ✗ REAL zqs(klon, klev)
1112 ✗ REAL zdqs(klon, klev)
1113 ✗ REAL ztemp(klon, klev)
1114 ✗ REAL zpres(klon, klev)
1115 ✗ REAL zconv(klon) ! convergence d'humidite
1116 ✗ REAL zvirt(klon) ! convergence virtuelle d'humidite
1117 ✗ REAL zfrac(klon) ! fraction convective
1118 ✗ INTEGER kb(klon), kh(klon)
1119
1120 ✗ REAL zcond(klon), zvapo(klon), zrapp(klon)
1121 ✗ REAL zrfl(klon), zrfln, zqev, zqevt
1122 REAL zdelta, zcvm5, zcor
1123 REAL zvar
1124
1125 LOGICAL appel1er
1126 SAVE appel1er
1127 !$OMP THREADPRIVATE(appel1er)
1128
1129 ! Fonctions thermodynamiques
1130
1131 include "YOETHF.h"
1132 include "FCTTRE.h"
1133
1134 DATA appel1er/.TRUE./
1135
1136 ✗ IF (appel1er) THEN
1137 ✗ PRINT *, 'conkuo, calcfcl:', calcfcl
1138 IF (.NOT. calcfcl) PRINT *, 'conkuo, ldepar:', ldepar
1139 ✗ PRINT *, 'conkuo, opt_cld:', opt_cld
1140 ✗ PRINT *, 'conkuo, evap_prec:', evap_prec
1141 ✗ PRINT *, 'conkuo, new_deh:', new_deh
1142 ✗ appel1er = .FALSE.
1143 END IF
1144
1145 ! Initialiser les sorties a zero
1146
1147 ✗ DO k = 1, klev
1148 ✗ DO i = 1, klon
1149 ✗ d_q(i, k) = 0.0
1150 ✗ d_t(i, k) = 0.0
1151 ✗ d_ql(i, k) = 0.0
1152 ✗ rneb(i, k) = 0.0
1153 END DO
1154 END DO
1155 ✗ DO i = 1, klon
1156 ✗ rain(i) = 0.0
1157 ✗ snow(i) = 0.0
1158 ✗ ibas(i) = 0
1159 ✗ itop(i) = 0
1160 END DO
1161
1162 ! Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT
1163
1164 ✗ DO k = 1, klev
1165 ✗ DO i = 1, klon
1166 ✗ IF (thermcep) THEN
1167 ✗ zdelta = max(0., sign(1.,rtt-t(i,k)))
1168 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1169 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
1170 ✗ zqs(i, k) = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
1171 ✗ zqs(i, k) = min(0.5, zqs(i,k))
1172 ✗ zcor = 1./(1.-retv*zqs(i,k))
1173 ✗ zqs(i, k) = zqs(i, k)*zcor
1174 ✗ zdqs(i, k) = foede(t(i,k), zdelta, zcvm5, zqs(i,k), zcor)
1175 ELSE
1176 IF (t(i,k)<t_coup) THEN
1177 zqs(i, k) = qsats(t(i,k))/pplay(i, k)
1178 zdqs(i, k) = dqsats(t(i,k), zqs(i,k))
1179 ELSE
1180 zqs(i, k) = qsatl(t(i,k))/pplay(i, k)
1181 zdqs(i, k) = dqsatl(t(i,k), zqs(i,k))
1182 END IF
1183 END IF
1184 END DO
1185 END DO
1186
1187 ! Calculer gz (energie potentielle)
1188
1189 ✗ DO i = 1, klon
1190 zgz(i, 1) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i, &
1191 ✗ 1)))*(paprs(i,1)-pplay(i,1))
1192 END DO
1193 ✗ DO k = 2, klev
1194 ✗ DO i = 1, klon
1195 zgz(i, k) = zgz(i, k-1) + rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i &
1196 ✗ ,k-1)-pplay(i,k))
1197 END DO
1198 END DO
1199
1200 ! Calculer l'energie statique humide saturee (Cp*T + gz + L*Qs)
1201
1202 ✗ DO k = 1, klev
1203 ✗ DO i = 1, klon
1204 ✗ ztotal(i, k) = rcpd*t(i, k) + rlvtt*zqs(i, k) + zgz(i, k)
1205 END DO
1206 END DO
1207
1208 ! Determiner le niveau de depart et calculer la difference de
1209 ! l'energie statique humide saturee (ztotal) entre la couche
1210 ! de depart et chaque couche au-dessus.
1211
1212 IF (calcfcl) THEN
1213 ✗ DO k = 1, klev
1214 ✗ DO i = 1, klon
1215 ✗ zpres(i, k) = pplay(i, k)
1216 ✗ ztemp(i, k) = t(i, k)
1217 END DO
1218 END DO
1219 ✗ CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb)
1220 ✗ DO i = 1, klon
1221 ✗ IF (ldcum(i)) THEN
1222 ✗ k = kb(i)
1223 IF (new_deh) THEN
1224 zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
1225 ELSE
1226 zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/ &
1227 paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
1228 ✗ rlvtt*(zqs(i,k-1)-zqs(i,k))
1229 END IF
1230 ✗ zdeh(i, k) = zdeh(i, k)*0.5
1231 END IF
1232 END DO
1233 ✗ DO k = 1, klev
1234 ✗ DO i = 1, klon
1235 ✗ IF (ldcum(i) .AND. k>=(kb(i)+1)) THEN
1236 ✗ IF (new_deh) THEN
1237 zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
1238 ELSE
1239 zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
1240 rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)* &
1241 ✗ (pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
1242 END IF
1243 END IF
1244 END DO
1245 END DO
1246 ELSE
1247 DO i = 1, klon
1248 k = ldepar
1249 kb(i) = ldepar
1250 ldcum(i) = .TRUE.
1251 IF (new_deh) THEN
1252 zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
1253 ELSE
1254 zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/paprs( &
1255 i, k)*(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
1256 END IF
1257 zdeh(i, k) = zdeh(i, k)*0.5
1258 END DO
1259 DO k = ldepar + 1, klev
1260 DO i = 1, klon
1261 IF (new_deh) THEN
1262 zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
1263 ELSE
1264 zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
1265 rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
1266 rlvtt*(zqs(i,k-1)-zqs(i,k))
1267 END IF
1268 END DO
1269 END DO
1270 END IF
1271
1272 ! -----Chercher le sommet du nuage
1273 ! -----Calculer la convergence de l'humidite (en kg/m**2 a un facteur
1274 ! -----psolpa/RG pres) du bas jusqu'au sommet du nuage.
1275 ! -----Calculer la convergence virtuelle pour que toute la maille
1276 ! -----deviennt nuageuse (du bas jusqu'au sommet du nuage)
1277
1278 ✗ DO i = 1, klon
1279 ✗ nuage(i) = .TRUE.
1280 ✗ zconv(i) = 0.0
1281 ✗ zvirt(i) = 0.0
1282 ✗ kh(i) = -999
1283 END DO
1284 ✗ DO k = 1, klev
1285 ✗ DO i = 1, klon
1286 ✗ IF (k>=kb(i) .AND. ldcum(i)) THEN
1287 ✗ nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
1288 ✗ IF (nuage(i)) THEN
1289 ✗ kh(i) = k
1290 ✗ zconv(i) = zconv(i) + conv_q(i, k)*dtime*(paprs(i,k)-paprs(i,k+1))
1291 zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k) &
1292 ✗ -paprs(i,k+1))
1293 END IF
1294 END IF
1295 END DO
1296 END DO
1297
1298 ✗ DO i = 1, klon
1299 ✗ todo(i) = ldcum(i) .AND. kh(i) > kb(i) .AND. zconv(i) > 0.0
1300 END DO
1301
1302 kbmin = klev
1303 kbmax = 0
1304 khmax = 0
1305 ✗ DO i = 1, klon
1306 ✗ IF (todo(i)) THEN
1307 ✗ kbmin = min(kbmin, kb(i))
1308 ✗ kbmax = max(kbmax, kb(i))
1309 ✗ khmax = max(khmax, kh(i))
1310 END IF
1311 END DO
1312
1313 ! -----Calculer la surface couverte par le nuage
1314
1315 ✗ DO i = 1, klon
1316 ✗ IF (todo(i)) THEN
1317 ✗ zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
1318 END IF
1319 END DO
1320
1321 ! -----Calculs essentiels:
1322
1323 ✗ DO i = 1, klon
1324 ✗ IF (todo(i)) THEN
1325 ✗ zcond(i) = 0.0
1326 END IF
1327 END DO
1328 ✗ DO k = kbmin, khmax
1329 ✗ DO i = 1, klon
1330 ✗ IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1331 ✗ zvar = zdeh(i, k)/(1.+zdqs(i,k))
1332 ✗ d_t(i, k) = zvar*zfrac(i)/rcpd
1333 d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
1334 ✗ conv_q(i, k)*dtime
1335 ✗ zcond(i) = zcond(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
1336 ✗ rneb(i, k) = zfrac(i)
1337 END IF
1338 END DO
1339 END DO
1340
1341 ✗ DO i = 1, klon
1342 ✗ IF (todo(i) .AND. zcond(i)<0.0) THEN
1343 ✗ PRINT *, 'WARNING: cond. negative (Kuo) ', i, kb(i), kh(i), zcond(i)
1344 ✗ zcond(i) = 0.0
1345 ✗ DO k = kb(i), kh(i)
1346 ✗ d_t(i, k) = 0.0
1347 ✗ d_q(i, k) = 0.0
1348 END DO
1349 ✗ todo(i) = .FALSE. ! effort totalement perdu
1350 END IF
1351 END DO
1352
1353 ! =====
1354 ! Une fois que la condensation a lieu, on doit construire un
1355 ! "modele nuageux" pour partager la condensation entre l'eau
1356 ! liquide nuageuse et la precipitation (leur rapport toliq
1357 ! est calcule selon l'epaisseur nuageuse). Je suppose que
1358 ! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
1359 ! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
1360 ! lineaire entre dpmin et dpmax).
1361 ! =====
1362 ✗ DO i = 1, klon
1363 ✗ IF (todo(i)) THEN
1364 toliq(i) = tomax - ((paprs(i,kb(i))-paprs(i,kh(i)+1))/paprs(i,1)-dpmin) &
1365 ✗ *(tomax-tomin)/(dpmax-dpmin)
1366 ✗ toliq(i) = max(tomin, min(tomax,toliq(i)))
1367 ✗ IF (pplay(i,kh(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
1368 IF (old_tau) toliq(i) = 1.0
1369 END IF
1370 END DO
1371 ! =====
1372 ! On doit aussi determiner la distribution verticale de
1373 ! l'eau nuageuse. Plusieurs options sont proposees:
1374
1375 ! (0) La condensation precipite integralement (toliq ne sera
1376 ! pas utilise).
1377 ! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
1378 ! a la vapeur d'eau locale.
1379 ! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
1380 ! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
1381 ! est effectivement diminuee pendant le processus d'ajustement.
1382 ! (4) Elle est en fonction (lineaire ou exponentielle) de la
1383 ! distance (epaisseur en pression) avec le niveau k1 (la couche
1384 ! k1 n'aura donc pas d'eau liquide).
1385 ! =====
1386
1387 IF (opt_cld==0) THEN
1388
1389 DO i = 1, klon
1390 IF (todo(i)) zrfl(i) = zcond(i)/dtime
1391 END DO
1392
1393 ELSE IF (opt_cld==1) THEN
1394
1395 DO i = 1, klon
1396 IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
1397 END DO
1398 DO k = kbmin, khmax
1399 DO i = 1, klon
1400 IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1401 zvapo(i) = zvapo(i) + (q(i,k)+d_q(i,k))*(paprs(i,k)-paprs(i,k+1))/ &
1402 rg
1403 END IF
1404 END DO
1405 END DO
1406 DO i = 1, klon
1407 IF (todo(i)) THEN
1408 zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
1409 zrapp(i) = max(0., min(1.,zrapp(i)))
1410 END IF
1411 END DO
1412 DO k = kbmin, khmax
1413 DO i = 1, klon
1414 IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1415 d_ql(i, k) = zrapp(i)*(q(i,k)+d_q(i,k))
1416 END IF
1417 END DO
1418 END DO
1419 DO i = 1, klon
1420 IF (todo(i)) THEN
1421 zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1422 END IF
1423 END DO
1424
1425 ELSE IF (opt_cld==2) THEN
1426
1427 DO i = 1, klon
1428 IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
1429 END DO
1430 DO k = kbmin, khmax
1431 DO i = 1, klon
1432 IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1433 zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg
1434 END IF
1435 END DO
1436 END DO
1437 DO k = kbmin, khmax
1438 DO i = 1, klon
1439 IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1440 d_ql(i, k) = toliq(i)*zcond(i)/zvapo(i)
1441 END IF
1442 END DO
1443 END DO
1444 DO i = 1, klon
1445 IF (todo(i)) THEN
1446 zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1447 END IF
1448 END DO
1449
1450 ELSE IF (opt_cld==3) THEN
1451
1452 DO i = 1, klon
1453 IF (todo(i)) THEN
1454 zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
1455 END IF
1456 END DO
1457 DO k = kbmin, khmax
1458 DO i = 1, klon
1459 IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
1460 zvapo(i) = zvapo(i) + max(0.0, -d_q(i,k))*(paprs(i,k)-paprs(i,k+1)) &
1461 /rg
1462 END IF
1463 END DO
1464 END DO
1465 DO k = kbmin, khmax
1466 DO i = 1, klon
1467 IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i) .AND. zvapo(i)>0.0) THEN
1468 d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*max(0.0, -d_q( &
1469 i,k))
1470 END IF
1471 END DO
1472 END DO
1473 DO i = 1, klon
1474 IF (todo(i)) THEN
1475 zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1476 END IF
1477 END DO
1478
1479 ELSE IF (opt_cld==4) THEN
1480
1481 nexpo = 3
1482 ! cc nexpo = 1 ! distribution lineaire
1483
1484 ✗ DO i = 1, klon
1485 ✗ IF (todo(i)) THEN
1486 ✗ zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation)
1487 END IF
1488 END DO
1489 ✗ DO k = kbmin, khmax
1490 ✗ DO i = 1, klon
1491 ✗ IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
1492 zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg*(pplay(i,kb(i))- &
1493 ✗ pplay(i,k))**nexpo
1494 END IF
1495 END DO
1496 END DO
1497 ✗ DO k = kbmin, khmax
1498 ✗ DO i = 1, klon
1499 ✗ IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
1500 d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*(pplay(i,kb(i) &
1501 ✗ )-pplay(i,k))**nexpo
1502 END IF
1503 END DO
1504 END DO
1505 ✗ DO i = 1, klon
1506 ✗ IF (todo(i)) THEN
1507 ✗ zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
1508 END IF
1509 END DO
1510
1511 ELSE ! valeur non-prevue pour opt_cld
1512
1513 PRINT *, 'opt_cld est faux:', opt_cld
1514 CALL abort
1515
1516 END IF ! fin de opt_cld
1517
1518 ! L'eau precipitante peut etre re-evaporee:
1519
1520 ✗ IF (evap_prec .AND. kbmax>=2) THEN
1521 ✗ DO k = kbmax, 1, -1
1522 ✗ DO i = 1, klon
1523 ✗ IF (todo(i) .AND. k<=(kb(i)-1) .AND. zrfl(i)>0.0) THEN
1524 ✗ zqev = max(0.0, (zqs(i,k)-q(i,k))*zfrac(i))
1525 zqevt = coef_eva*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i))* &
1526 ✗ (paprs(i,k)-paprs(i,k+1))/pplay(i, k)*t(i, k)*rd/rg
1527 zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/ &
1528 ✗ (paprs(i,k)-paprs(i,k+1))
1529 ✗ zqev = min(zqev, zqevt)
1530 ✗ zrfln = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
1531 ✗ d_q(i, k) = -(zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i,k+1)))*dtime
1532 d_t(i, k) = (zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
1533 ✗ k+1)))*dtime*rlvtt/rcpd
1534 ✗ zrfl(i) = zrfln
1535 END IF
1536 END DO
1537 END DO
1538 END IF
1539
1540 ! La temperature de la premiere couche determine la pluie ou la neige:
1541
1542 ✗ DO i = 1, klon
1543 ✗ IF (todo(i)) THEN
1544 ✗ IF (t(i,1)>rtt) THEN
1545 ✗ rain(i) = rain(i) + zrfl(i)
1546 ELSE
1547 ✗ snow(i) = snow(i) + zrfl(i)
1548 END IF
1549 END IF
1550 END DO
1551
1552 ✗ RETURN
1553 END SUBROUTINE conkuo
1554 ✗ SUBROUTINE kuofcl(pt, pq, pg, pp, ldcum, kcbot)
1555 ✗ USE dimphy
1556 IMPLICIT NONE
1557 ! ======================================================================
1558 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1559 ! adaptation du code de Tiedtke du ECMWF
1560 ! Objet: calculer le niveau de convection libre
1561 ! (FCL: Free Convection Level)
1562 ! ======================================================================
1563 ! Arguments:
1564 ! pt---input-R- temperature (K)
1565 ! pq---input-R- vapeur d'eau (kg/kg)
1566 ! pg---input-R- geopotentiel (g*z ou z est en metre)
1567 ! pp---input-R- pression (Pa)
1568
1569 ! LDCUM---output-L- Y-t-il la convection
1570 ! kcbot---output-I- Niveau du bas de la convection
1571 ! ======================================================================
1572 include "YOMCST.h"
1573 include "YOETHF.h"
1574
1575 REAL pt(klon, klev), pq(klon, klev), pg(klon, klev), pp(klon, klev)
1576 INTEGER kcbot(klon)
1577 LOGICAL ldcum(klon)
1578
1579 ✗ REAL ztu(klon, klev), zqu(klon, klev), zlu(klon, klev)
1580 ✗ REAL zqold(klon), zbuo
1581 INTEGER is, i, k
1582
1583 ! klab=1: on est sous le nuage convectif
1584 ! klab=2: le bas du nuage convectif
1585 ! klab=0: autres couches
1586 ✗ INTEGER klab(klon, klev)
1587
1588 ! quand lflag=.true., on est sous le nuage, il faut donc appliquer
1589 ! le processus d'elevation.
1590 ✗ LOGICAL lflag(klon)
1591
1592 ✗ DO k = 1, klev
1593 ✗ DO i = 1, klon
1594 ✗ ztu(i, k) = pt(i, k)
1595 ✗ zqu(i, k) = pq(i, k)
1596 ✗ zlu(i, k) = 0.0
1597 ✗ klab(i, k) = 0
1598 END DO
1599 END DO
1600 ! ----------------------------------------------------------------------
1601 ✗ DO i = 1, klon
1602 ✗ klab(i, 1) = 1
1603 ✗ kcbot(i) = 2
1604 ✗ ldcum(i) = .FALSE.
1605 END DO
1606
1607 ✗ DO k = 2, klev - 1
1608
1609 is = 0
1610 ✗ DO i = 1, klon
1611 ✗ IF (klab(i,k-1)==1) is = is + 1
1612 ✗ lflag(i) = .FALSE.
1613 ✗ IF (klab(i,k-1)==1) lflag(i) = .TRUE.
1614 END DO
1615 ✗ IF (is==0) GO TO 290
1616
1617 ! on eleve le parcel d'air selon l'adiabatique sec
1618
1619 ✗ DO i = 1, klon
1620 ✗ IF (lflag(i)) THEN
1621 ✗ zqu(i, k) = zqu(i, k-1)
1622 ✗ ztu(i, k) = ztu(i, k-1) + (pg(i,k-1)-pg(i,k))/rcpd
1623 zbuo = ztu(i, k)*(1.+retv*zqu(i,k)) - pt(i, k)*(1.+retv*pq(i,k)) + &
1624 ✗ 0.5
1625 ✗ IF (zbuo>0.) klab(i, k) = 1
1626 ✗ zqold(i) = zqu(i, k)
1627 END IF
1628 END DO
1629
1630 ! on calcule la condensation eventuelle
1631
1632 ✗ CALL adjtq(pp(1,k), ztu(1,k), zqu(1,k), lflag, 1)
1633
1634 ! s'il y a la condensation et la "buoyancy" force est positive
1635 ! c'est bien le bas de la tour de convection
1636
1637 ✗ DO i = 1, klon
1638 ✗ IF (lflag(i) .AND. zqu(i,k)/=zqold(i)) THEN
1639 ✗ klab(i, k) = 2
1640 ✗ zlu(i, k) = zlu(i, k) + zqold(i) - zqu(i, k)
1641 zbuo = ztu(i, k)*(1.+retv*zqu(i,k)) - pt(i, k)*(1.+retv*pq(i,k)) + &
1642 ✗ 0.5
1643 ✗ IF (zbuo>0.) THEN
1644 ✗ kcbot(i) = k
1645 ✗ ldcum(i) = .TRUE.
1646 END IF
1647 END IF
1648 END DO
1649
1650 ✗ 290 END DO
1651
1652 ✗ RETURN
1653 END SUBROUTINE kuofcl
1654 ✗ SUBROUTINE adjtq(pp, pt, pq, ldflag, kcall)
1655 USE dimphy
1656 IMPLICIT NONE
1657 ! ======================================================================
1658 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1659 ! adaptation du code de Tiedtke du ECMWF
1660 ! Objet: ajustement entre T et Q
1661 ! ======================================================================
1662 ! Arguments:
1663 ! pp---input-R- pression (Pa)
1664 ! pt---input/output-R- temperature (K)
1665 ! pq---input/output-R- vapeur d'eau (kg/kg)
1666 ! ======================================================================
1667 ! TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT
1668
1669 ! NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS
1670 ! KCALL=0 ENV. T AND QS IN*CUINI*
1671 ! KCALL=1 CONDENSATION IN UPDRAFTS (E.G. CUBASE, CUASC)
1672 ! KCALL=2 EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1673
1674 include "YOMCST.h"
1675
1676 REAL pt(klon), pq(klon), pp(klon)
1677 LOGICAL ldflag(klon)
1678 INTEGER kcall
1679
1680 REAL t_coup
1681 PARAMETER (t_coup=234.0)
1682
1683 ✗ REAL zcond(klon), zcond1
1684 REAL zdelta, zcvm5, zldcp, zqsat, zcor, zdqsat
1685 INTEGER is, i
1686 include "YOETHF.h"
1687 include "FCTTRE.h"
1688
1689 ✗ DO i = 1, klon
1690 ✗ zcond(i) = 0.0
1691 END DO
1692
1693 ✗ DO i = 1, klon
1694 ✗ IF (ldflag(i)) THEN
1695 ✗ zdelta = max(0., sign(1.,rtt-pt(i)))
1696 ✗ zldcp = rlvtt*(1.-zdelta) + zdelta*rlstt
1697 ✗ zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
1698 IF (thermcep) THEN
1699 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1700 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
1701 ✗ zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
1702 ✗ zqsat = min(0.5, zqsat)
1703 ✗ zcor = 1./(1.-retv*zqsat)
1704 ✗ zqsat = zqsat*zcor
1705 ✗ zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
1706 ELSE
1707 IF (pt(i)<t_coup) THEN
1708 zqsat = qsats(pt(i))/pp(i)
1709 zdqsat = dqsats(pt(i), zqsat)
1710 ELSE
1711 zqsat = qsatl(pt(i))/pp(i)
1712 zdqsat = dqsatl(pt(i), zqsat)
1713 END IF
1714 END IF
1715 ✗ zcond(i) = (pq(i)-zqsat)/(1.+zdqsat)
1716 ✗ IF (kcall==1) zcond(i) = max(zcond(i), 0.)
1717 ✗ IF (kcall==2) zcond(i) = min(zcond(i), 0.)
1718 ✗ pt(i) = pt(i) + zldcp*zcond(i)
1719 ✗ pq(i) = pq(i) - zcond(i)
1720 END IF
1721 END DO
1722
1723 is = 0
1724 ✗ DO i = 1, klon
1725 ✗ IF (zcond(i)/=0.) is = is + 1
1726 END DO
1727 ✗ IF (is==0) GO TO 230
1728
1729 ✗ DO i = 1, klon
1730 ✗ IF (ldflag(i) .AND. zcond(i)/=0.) THEN
1731 ✗ zdelta = max(0., sign(1.,rtt-pt(i)))
1732 ✗ zldcp = rlvtt*(1.-zdelta) + zdelta*rlstt
1733 ✗ zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
1734 IF (thermcep) THEN
1735 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1736 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
1737 ✗ zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
1738 ✗ zqsat = min(0.5, zqsat)
1739 ✗ zcor = 1./(1.-retv*zqsat)
1740 ✗ zqsat = zqsat*zcor
1741 ✗ zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
1742 ELSE
1743 IF (pt(i)<t_coup) THEN
1744 zqsat = qsats(pt(i))/pp(i)
1745 zdqsat = dqsats(pt(i), zqsat)
1746 ELSE
1747 zqsat = qsatl(pt(i))/pp(i)
1748 zdqsat = dqsatl(pt(i), zqsat)
1749 END IF
1750 END IF
1751 ✗ zcond1 = (pq(i)-zqsat)/(1.+zdqsat)
1752 ✗ pt(i) = pt(i) + zldcp*zcond1
1753 ✗ pq(i) = pq(i) - zcond1
1754 END IF
1755 END DO
1756
1757 230 CONTINUE
1758 ✗ RETURN
1759 END SUBROUTINE adjtq
1760 ✗ SUBROUTINE fiajh(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, &
1761 ✗ ibas, itop)
1762 USE dimphy
1763 IMPLICIT NONE
1764
1765 ! Ajustement humide (Schema de convection de Manabe)
1766 ! .
1767 include "YOMCST.h"
1768
1769 ! Arguments:
1770
1771 REAL dtime ! intervalle du temps (s)
1772 REAL t(klon, klev) ! temperature (K)
1773 REAL q(klon, klev) ! humidite specifique (kg/kg)
1774 REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
1775 REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
1776
1777 REAL d_t(klon, klev) ! incrementation pour la temperature
1778 REAL d_q(klon, klev) ! incrementation pour vapeur d'eau
1779 REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
1780 REAL rneb(klon, klev) ! fraction nuageuse
1781
1782 REAL rain(klon) ! variable non utilisee
1783 REAL snow(klon) ! variable non utilisee
1784 INTEGER ibas(klon) ! variable non utilisee
1785 INTEGER itop(klon) ! variable non utilisee
1786
1787 REAL t_coup
1788 PARAMETER (t_coup=234.0)
1789 REAL seuil_vap
1790 PARAMETER (seuil_vap=1.0E-10)
1791
1792 ! Variables locales:
1793
1794 INTEGER i, k
1795 INTEGER k1, k1p, k2, k2p
1796 ✗ LOGICAL itest(klon)
1797 ✗ REAL delta_q(klon, klev)
1798 ✗ REAL cp_new_t(klev)
1799 ✗ REAL cp_delta_t(klev)
1800 ✗ REAL new_qb(klev)
1801 ✗ REAL v_cptj(klev), v_cptjk1, v_ssig
1802 ✗ REAL v_cptt(klon, klev), v_p, v_t
1803 ✗ REAL v_qs(klon, klev), v_qsd(klon, klev)
1804 ✗ REAL zq1(klon), zq2(klon)
1805 ✗ REAL gamcpdz(klon, 2:klev)
1806 REAL zdp, zdpm
1807
1808 REAL zsat ! sur-saturation
1809 REAL zflo ! flotabilite
1810
1811 ✗ REAL local_q(klon, klev), local_t(klon, klev)
1812
1813 REAL zdelta, zcor, zcvm5
1814
1815 include "YOETHF.h"
1816 include "FCTTRE.h"
1817
1818 ✗ DO k = 1, klev
1819 ✗ DO i = 1, klon
1820 ✗ local_q(i, k) = q(i, k)
1821 ✗ local_t(i, k) = t(i, k)
1822 ✗ rneb(i, k) = 0.0
1823 ✗ d_ql(i, k) = 0.0
1824 ✗ d_t(i, k) = 0.0
1825 ✗ d_q(i, k) = 0.0
1826 END DO
1827 END DO
1828 ✗ DO i = 1, klon
1829 ✗ rain(i) = 0.0
1830 ✗ snow(i) = 0.0
1831 ✗ ibas(i) = 0
1832 ✗ itop(i) = 0
1833 END DO
1834
1835 ! Calculer v_qs et v_qsd:
1836
1837 ✗ DO k = 1, klev
1838 ✗ DO i = 1, klon
1839 ✗ v_cptt(i, k) = rcpd*local_t(i, k)
1840 v_t = local_t(i, k)
1841 ✗ v_p = pplay(i, k)
1842
1843 ✗ IF (thermcep) THEN
1844 ✗ zdelta = max(0., sign(1.,rtt-v_t))
1845 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1846 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
1847 ✗ v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
1848 ✗ v_qs(i, k) = min(0.5, v_qs(i,k))
1849 ✗ zcor = 1./(1.-retv*v_qs(i,k))
1850 ✗ v_qs(i, k) = v_qs(i, k)*zcor
1851 ✗ v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
1852 ELSE
1853 IF (v_t<t_coup) THEN
1854 v_qs(i, k) = qsats(v_t)/v_p
1855 v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
1856 ELSE
1857 v_qs(i, k) = qsatl(v_t)/v_p
1858 v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
1859 END IF
1860 END IF
1861 END DO
1862 END DO
1863
1864 ! Calculer Gamma * Cp * dz: (gamm est le gradient critique)
1865
1866 ✗ DO k = 2, klev
1867 ✗ DO i = 1, klon
1868 ✗ zdp = paprs(i, k) - paprs(i, k+1)
1869 ✗ zdpm = paprs(i, k-1) - paprs(i, k)
1870 gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)*(v_cptt(i,k-1)*zdpm+ &
1871 v_cptt(i,k)*zdp)+rlvtt/(zdpm+zdp)*(v_qs(i,k-1)*zdpm+ &
1872 v_qs(i,k)*zdp))*(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
1873 ✗ k-1)*zdpm+v_qsd(i,k)*zdp)/(zdpm+zdp))
1874 END DO
1875 END DO
1876
1877 ! ------------------------------------ modification des profils instables
1878 ✗ DO i = 1, klon
1879 ✗ itest(i) = .FALSE.
1880
1881 k1 = 0
1882 ✗ k2 = 1
1883
1884 810 CONTINUE ! chercher k1, le bas de la colonne
1885 ✗ k2 = k2 + 1
1886 ✗ IF (k2>klev) GO TO 9999
1887 ✗ zflo = v_cptt(i, k2-1) - v_cptt(i, k2) - gamcpdz(i, k2)
1888 zsat = (local_q(i,k2-1)-v_qs(i,k2-1))*(paprs(i,k2-1)-paprs(i,k2)) + &
1889 ✗ (local_q(i,k2)-v_qs(i,k2))*(paprs(i,k2)-paprs(i,k2+1))
1890 ✗ IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 810
1891 k1 = k2 - 1
1892 ✗ itest(i) = .TRUE.
1893
1894 820 CONTINUE ! chercher k2, le haut de la colonne
1895 ✗ IF (k2==klev) GO TO 821
1896 ✗ k2p = k2 + 1
1897 ✗ zsat = zsat + (paprs(i,k2p)-paprs(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
1898 ✗ zflo = v_cptt(i, k2p-1) - v_cptt(i, k2p) - gamcpdz(i, k2p)
1899 ✗ IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 821
1900 k2 = k2p
1901 ✗ GO TO 820
1902 821 CONTINUE
1903
1904 ! ------------------------------------------------------ ajustement local
1905 830 CONTINUE ! ajustement proprement dit
1906 ✗ v_cptj(k1) = 0.0
1907 ✗ zdp = paprs(i, k1) - paprs(i, k1+1)
1908 v_cptjk1 = ((1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1))+rlvtt*(local_q(i, &
1909 ✗ k1)-v_qs(i,k1)))*zdp
1910 ✗ v_ssig = zdp*(1.0+v_qsd(i,k1))
1911
1912 k1p = k1 + 1
1913 ✗ DO k = k1p, k2
1914 ✗ zdp = paprs(i, k) - paprs(i, k+1)
1915 ✗ v_cptj(k) = v_cptj(k-1) + gamcpdz(i, k)
1916 v_cptjk1 = v_cptjk1 + zdp*((1.0+v_qsd(i,k))*(v_cptt(i, &
1917 ✗ k)+v_cptj(k))+rlvtt*(local_q(i,k)-v_qs(i,k)))
1918 ✗ v_ssig = v_ssig + zdp*(1.0+v_qsd(i,k))
1919 END DO
1920
1921 ✗ DO k = k1, k2
1922 ✗ cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
1923 ✗ cp_delta_t(k) = cp_new_t(k) - v_cptt(i, k)
1924 ✗ new_qb(k) = v_qs(i, k) + v_qsd(i, k)*cp_delta_t(k)/rlvtt
1925 ✗ local_q(i, k) = new_qb(k)
1926 ✗ local_t(i, k) = cp_new_t(k)/rcpd
1927 END DO
1928
1929 ! --------------------------------------------------- sondage vers le bas
1930 ! -- on redefinit les variables prognostiques dans
1931 ! -- la colonne qui vient d'etre ajustee
1932
1933 ✗ DO k = k1, k2
1934 ✗ v_cptt(i, k) = rcpd*local_t(i, k)
1935 v_t = local_t(i, k)
1936 ✗ v_p = pplay(i, k)
1937
1938 ✗ IF (thermcep) THEN
1939 ✗ zdelta = max(0., sign(1.,rtt-v_t))
1940 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
1941 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
1942 ✗ v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
1943 ✗ v_qs(i, k) = min(0.5, v_qs(i,k))
1944 ✗ zcor = 1./(1.-retv*v_qs(i,k))
1945 ✗ v_qs(i, k) = v_qs(i, k)*zcor
1946 ✗ v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
1947 ELSE
1948 IF (v_t<t_coup) THEN
1949 v_qs(i, k) = qsats(v_t)/v_p
1950 v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
1951 ELSE
1952 v_qs(i, k) = qsatl(v_t)/v_p
1953 v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
1954 END IF
1955 END IF
1956 END DO
1957 ✗ DO k = 2, klev
1958 ✗ zdpm = paprs(i, k-1) - paprs(i, k)
1959 ✗ zdp = paprs(i, k) - paprs(i, k+1)
1960 gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)*(v_cptt(i,k-1)*zdpm+ &
1961 v_cptt(i,k)*zdp)+rlvtt/(zdpm+zdp)*(v_qs(i,k-1)*zdpm+ &
1962 v_qs(i,k)*zdp))*(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
1963 ✗ k-1)*zdpm+v_qsd(i,k)*zdp)/(zdpm+zdp))
1964 END DO
1965
1966 ! Verifier si l'on peut etendre la colonne vers le bas
1967
1968 ✗ IF (k1==1) GO TO 841 ! extension echouee
1969 ✗ zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
1970 zsat = (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1)) + &
1971 ✗ (local_q(i,k1)-v_qs(i,k1))*(paprs(i,k1)-paprs(i,k1+1))
1972 ✗ IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 841 ! extension echouee
1973
1974 840 CONTINUE
1975 ✗ k1 = k1 - 1
1976 ✗ IF (k1==1) GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
1977 ✗ zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1))
1978 ✗ zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
1979 ✗ IF (zflo>0.0 .AND. zsat>0.0) THEN
1980 GO TO 840
1981 ELSE
1982 ✗ GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
1983 END IF
1984 841 CONTINUE
1985
1986 ✗ GO TO 810 ! chercher d'autres blocks en haut
1987
1988 ✗ 9999 END DO ! boucle sur tous les points
1989 ! -----------------------------------------------------------------------
1990
1991 ! Determiner la fraction nuageuse (hypothese: la nebulosite a lieu
1992 ! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
1993
1994 ✗ DO k = 1, klev
1995 ✗ DO i = 1, klon
1996 ✗ IF (itest(i)) THEN
1997 ✗ delta_q(i, k) = local_q(i, k) - q(i, k)
1998 ✗ IF (delta_q(i,k)<0.) rneb(i, k) = 1.0
1999 END IF
2000 END DO
2001 END DO
2002
2003 ! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
2004 ! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
2005 ! diminue et d'une maniere proportionnelle a cet diminution):
2006
2007 ✗ DO i = 1, klon
2008 ✗ IF (itest(i)) THEN
2009 ✗ zq1(i) = 0.0
2010 ✗ zq2(i) = 0.0
2011 END IF
2012 END DO
2013 ✗ DO k = 1, klev
2014 ✗ DO i = 1, klon
2015 ✗ IF (itest(i)) THEN
2016 ✗ zdp = paprs(i, k) - paprs(i, k+1)
2017 ✗ zq1(i) = zq1(i) - delta_q(i, k)*zdp
2018 ✗ zq2(i) = zq2(i) - min(0.0, delta_q(i,k))*zdp
2019 END IF
2020 END DO
2021 END DO
2022 ✗ DO k = 1, klev
2023 ✗ DO i = 1, klon
2024 ✗ IF (itest(i)) THEN
2025 ✗ IF (zq2(i)/=0.0) d_ql(i, k) = -min(0.0, delta_q(i,k))*zq1(i)/zq2(i)
2026 END IF
2027 END DO
2028 END DO
2029
2030 ✗ DO k = 1, klev
2031 ✗ DO i = 1, klon
2032 ✗ local_q(i, k) = max(local_q(i,k), seuil_vap)
2033 END DO
2034 END DO
2035
2036 ✗ DO k = 1, klev
2037 ✗ DO i = 1, klon
2038 ✗ d_t(i, k) = local_t(i, k) - t(i, k)
2039 ✗ d_q(i, k) = local_q(i, k) - q(i, k)
2040 END DO
2041 END DO
2042
2043 ✗ RETURN
2044 END SUBROUTINE fiajh
2045 ✗ SUBROUTINE fiajc(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
2046 ✗ rain, snow, ibas, itop)
2047 USE dimphy
2048 IMPLICIT NONE
2049
2050 include "YOMCST.h"
2051
2052 ! Options:
2053
2054 INTEGER plb ! niveau de depart pour la convection
2055 PARAMETER (plb=4)
2056
2057 ! Mystere: cette option n'est pas innocente pour les resultats !
2058 ! Qui peut resoudre ce mystere ? (Z.X.Li mars 1995)
2059 LOGICAL vector ! calcul vectorise
2060 PARAMETER (vector=.FALSE.)
2061
2062 REAL t_coup
2063 PARAMETER (t_coup=234.0)
2064
2065 ! Arguments:
2066
2067 REAL q(klon, klev) ! humidite specifique (kg/kg)
2068 REAL t(klon, klev) ! temperature (K)
2069 REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
2070 REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
2071 REAL dtime ! intervalle du temps (s)
2072 REAL conv_q(klon, klev) ! taux de convergence de l'humidite
2073 REAL rneb(klon, klev) ! fraction nuageuse
2074 REAL d_q(klon, klev) ! incrementaion pour la vapeur d'eau
2075 REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
2076 REAL d_t(klon, klev) ! incrementation pour la temperature
2077 REAL rain(klon) ! variable non-utilisee
2078 REAL snow(klon) ! variable non-utilisee
2079 INTEGER itop(klon) ! variable non-utilisee
2080 INTEGER ibas(klon) ! variable non-utilisee
2081
2082 ✗ INTEGER kh(klon), i, k
2083 ✗ LOGICAL nuage(klon), test(klon, klev)
2084 ✗ REAL zconv(klon), zdeh(klon, klev), zvirt(klon)
2085 ✗ REAL zdqs(klon, klev), zqs(klon, klev)
2086 ✗ REAL ztt, zvar, zfrac(klon)
2087 ✗ REAL zq1(klon), zq2(klon)
2088 REAL zdelta, zcor, zcvm5
2089
2090 include "YOETHF.h"
2091 include "FCTTRE.h"
2092
2093 ! Initialiser les sorties:
2094
2095 ✗ DO k = 1, klev
2096 ✗ DO i = 1, klon
2097 ✗ rneb(i, k) = 0.0
2098 ✗ d_ql(i, k) = 0.0
2099 ✗ d_t(i, k) = 0.0
2100 ✗ d_q(i, k) = 0.0
2101 END DO
2102 END DO
2103 ✗ DO i = 1, klon
2104 ✗ itop(i) = 0
2105 ✗ ibas(i) = 0
2106 ✗ rain(i) = 0.0
2107 ✗ snow(i) = 0.0
2108 END DO
2109
2110 ! Calculer Qs et L/Cp * dQs/dT:
2111
2112 ✗ DO k = 1, klev
2113 ✗ DO i = 1, klon
2114 ✗ ztt = t(i, k)
2115 ✗ IF (thermcep) THEN
2116 ✗ zdelta = max(0., sign(1.,rtt-ztt))
2117 ✗ zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
2118 ✗ zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
2119 ✗ zqs(i, k) = r2es*foeew(ztt, zdelta)/pplay(i, k)
2120 ✗ zqs(i, k) = min(0.5, zqs(i,k))
2121 ✗ zcor = 1./(1.-retv*zqs(i,k))
2122 ✗ zqs(i, k) = zqs(i, k)*zcor
2123 ✗ zdqs(i, k) = foede(ztt, zdelta, zcvm5, zqs(i,k), zcor)
2124 ELSE
2125 IF (ztt<t_coup) THEN
2126 zqs(i, k) = qsats(ztt)/pplay(i, k)
2127 zdqs(i, k) = dqsats(ztt, zqs(i,k))
2128 ELSE
2129 zqs(i, k) = qsatl(ztt)/pplay(i, k)
2130 zdqs(i, k) = dqsatl(ztt, zqs(i,k))
2131 END IF
2132 END IF
2133 END DO
2134 END DO
2135
2136 ! Determiner la difference de l'energie totale saturee:
2137
2138 ✗ DO i = 1, klon
2139 k = plb
2140 zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k &
2141 ✗ )*(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
2142 ✗ zdeh(i, k) = zdeh(i, k)*0.5 ! on prend la moitie
2143 END DO
2144 ✗ DO k = plb + 1, klev
2145 ✗ DO i = 1, klon
2146 zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
2147 rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
2148 ✗ rlvtt*(zqs(i,k-1)-zqs(i,k))
2149 END DO
2150 END DO
2151
2152 ! Determiner le sommet du nuage selon l'instabilite
2153 ! Calculer les convergences d'humidite (reelle et virtuelle)
2154
2155 ✗ DO i = 1, klon
2156 ✗ nuage(i) = .TRUE.
2157 ✗ zconv(i) = 0.0
2158 ✗ zvirt(i) = 0.0
2159 ✗ kh(i) = -999
2160 END DO
2161 ✗ DO k = plb, klev
2162 ✗ DO i = 1, klon
2163 ✗ nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
2164 ✗ IF (nuage(i)) THEN
2165 ✗ kh(i) = k
2166 ✗ zconv(i) = zconv(i) + conv_q(i, k)*dtime*(paprs(i,k)-paprs(i,k+1))
2167 zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k)- &
2168 ✗ paprs(i,k+1))
2169 END IF
2170 END DO
2171 END DO
2172
2173 IF (vector) THEN
2174
2175
2176 DO k = plb, klev
2177 DO i = 1, klon
2178 IF (k<=kh(i) .AND. kh(i)>plb .AND. zconv(i)>0.0) THEN
2179 test(i, k) = .TRUE.
2180 zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
2181 ELSE
2182 test(i, k) = .FALSE.
2183 END IF
2184 END DO
2185 END DO
2186
2187 DO k = plb, klev
2188 DO i = 1, klon
2189 IF (test(i,k)) THEN
2190 zvar = zdeh(i, k)/(1.0+zdqs(i,k))
2191 d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
2192 conv_q(i, k)*dtime
2193 d_t(i, k) = zvar*zfrac(i)/rcpd
2194 END IF
2195 END DO
2196 END DO
2197
2198 DO i = 1, klon
2199 zq1(i) = 0.0
2200 zq2(i) = 0.0
2201 END DO
2202 DO k = plb, klev
2203 DO i = 1, klon
2204 IF (test(i,k)) THEN
2205 IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
2206 zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
2207 zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
2208 END IF
2209 END DO
2210 END DO
2211
2212 DO k = plb, klev
2213 DO i = 1, klon
2214 IF (test(i,k)) THEN
2215 IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
2216 END IF
2217 END DO
2218 END DO
2219
2220 ELSE ! (.NOT. vector)
2221
2222 ✗ DO i = 1, klon
2223 ✗ IF (kh(i)>plb .AND. zconv(i)>0.0) THEN
2224 ! cc IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite
2225 ! cc IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante
2226 ✗ zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
2227 ✗ DO k = plb, kh(i)
2228 ✗ zvar = zdeh(i, k)/(1.0+zdqs(i,k))
2229 d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
2230 ✗ conv_q(i, k)*dtime
2231 ✗ d_t(i, k) = zvar*zfrac(i)/rcpd
2232 END DO
2233
2234 ✗ zq1(i) = 0.0
2235 ✗ zq2(i) = 0.0
2236 ✗ DO k = plb, kh(i)
2237 ✗ IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
2238 ✗ zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
2239 ✗ zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
2240 END DO
2241 ✗ DO k = plb, kh(i)
2242 ✗ IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
2243 END DO
2244 END IF
2245 END DO
2246
2247 END IF ! fin de teste sur vector
2248
2249 ✗ RETURN
2250 END SUBROUTINE fiajc
2251