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