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