source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/phylmd/conlmd.F @ 1671

Last change on this file since 1671 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

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