source: LMDZ.3.3/trunk/libf/phylmd/conflx.F @ 5415

Last change on this file since 5415 was 343, checked in by lmdz, 23 years ago

On corrige le double bug sur la vitesse verticale
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 55.2 KB
Line 
1      SUBROUTINE conflx (dtime,pres_h,pres_f,
2     e                   t, q, con_t, con_q, pqhfl, w,
3     s                   d_t, d_q, rain, snow,
4     s                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
5     s                   kcbot, kctop, kdtop, pmflxr, pmflxs)
6c
7      IMPLICIT none
8c======================================================================
9c Auteur(s): Z.X. Li (LMD/CNRS) date: 19941014
10c Objet: Schema flux de masse pour la convection
11c        (schema de Tiedtke avec qqs modifications mineures)
12c Dec.97: Prise en compte des modifications introduites par
13c         Olivier Boucher et Alexandre Armengaud pour melange
14c         et lessivage des traceurs passifs.
15c======================================================================
16#include "dimensions.h"
17#include "dimphy.h"
18#include "YOMCST.h"
19#include "YOETHF.h"
20c Entree:
21      REAL dtime            ! pas d'integration (s)
22      REAL pres_h(klon,klev+1) ! pression half-level (Pa)
23      REAL pres_f(klon,klev)! pression full-level (Pa)
24      REAL t(klon,klev)     ! temperature (K)
25      REAL q(klon,klev)     ! humidite specifique (g/g)
26      REAL w(klon,klev)     ! vitesse verticale (Pa/s)
27      REAL con_t(klon,klev) ! convergence de temperature (K/s)
28      REAL con_q(klon,klev) ! convergence de l'eau vapeur (g/g/s)
29      REAL pqhfl(klon)      ! evaporation (negative vers haut) mm/s
30c Sortie:
31      REAL d_t(klon,klev)   ! incrementation de temperature
32      REAL d_q(klon,klev)   ! incrementation d'humidite
33      REAL pmfu(klon,klev)  ! flux masse (kg/m2/s) panache ascendant
34      REAL pmfd(klon,klev)  ! flux masse (kg/m2/s) panache descendant
35      REAL pen_u(klon,klev)
36      REAL pen_d(klon,klev)
37      REAL pde_u(klon,klev)
38      REAL pde_d(klon,klev)
39      REAL rain(klon)       ! pluie (mm/s)
40      REAL snow(klon)       ! neige (mm/s)
41      REAL pmflxr(klon,klev+1)
42      REAL pmflxs(klon,klev+1)
43      INTEGER kcbot(klon)  ! niveau du bas de la convection
44      INTEGER kctop(klon)  ! niveau du haut de la convection
45      INTEGER kdtop(klon)  ! niveau du haut des downdrafts
46c Local:
47      REAL pt(klon,klev)
48      REAL pq(klon,klev)
49      REAL pqs(klon,klev)
50      REAL pvervel(klon,klev)
51      LOGICAL land(klon)
52c
53      REAL d_t_bis(klon,klev)
54      REAL d_q_bis(klon,klev)
55      REAL paprs(klon,klev+1)
56      REAL paprsf(klon,klev)
57      REAL zgeom(klon,klev)
58      REAL zcvgq(klon,klev)
59      REAL zcvgt(klon,klev)
60cAA
61      REAL zmfu(klon,klev)
62      REAL zmfd(klon,klev)
63      REAL zen_u(klon,klev)
64      REAL zen_d(klon,klev)
65      REAL zde_u(klon,klev)
66      REAL zde_d(klon,klev)
67      REAL zmflxr(klon,klev+1)
68      REAL zmflxs(klon,klev+1)
69cAA
70
71c
72      INTEGER i, k
73      REAL zdelta, zqsat
74c
75#include "FCTTRE.h"
76c
77c initialiser les variables de sortie (pour securite)
78      DO i = 1, klon
79         rain(i) = 0.0
80         snow(i) = 0.0
81         kcbot(i) = 0
82         kctop(i) = 0
83         kdtop(i) = 0
84      ENDDO
85      DO k = 1, klev
86      DO i = 1, klon
87         d_t(i,k) = 0.0
88         d_q(i,k) = 0.0
89         pmfu(i,k) = 0.0
90         pmfd(i,k) = 0.0
91         pen_u(i,k) = 0.0
92         pde_u(i,k) = 0.0
93         pen_d(i,k) = 0.0
94         pde_d(i,k) = 0.0
95         zmfu(i,k) = 0.0
96         zmfd(i,k) = 0.0
97         zen_u(i,k) = 0.0
98         zde_u(i,k) = 0.0
99         zen_d(i,k) = 0.0
100         zde_d(i,k) = 0.0
101      ENDDO
102      ENDDO
103      DO k = 1, klev+1
104      DO i = 1, klon
105         zmflxr(i,k) = 0.0
106         zmflxs(i,k) = 0.0
107      ENDDO
108      ENDDO
109c
110c calculer la nature du sol (pour l'instant, ocean partout)
111      DO i = 1, klon
112         land(i) = .FALSE.
113      ENDDO
114c
115c preparer les variables d'entree (attention: l'ordre des niveaux
116c verticaux augmente du haut vers le bas)
117      DO k = 1, klev
118      DO i = 1, klon
119         pt(i,k) = t(i,klev-k+1)
120         pq(i,k) = q(i,klev-k+1)
121         paprsf(i,k) = pres_f(i,klev-k+1)
122         paprs(i,k) = pres_h(i,klev+1-k+1)
123         pvervel(i,k) = w(i,klev+1-k)
124         zcvgt(i,k) = con_t(i,klev-k+1)
125         zcvgq(i,k) = con_q(i,klev-k+1)
126c
127         zdelta=MAX(0.,SIGN(1.,RTT-pt(i,k)))
128         zqsat=R2ES*FOEEW ( pt(i,k), zdelta ) / paprsf(i,k)
129         zqsat=MIN(0.5,zqsat)
130         zqsat=zqsat/(1.-RETV  *zqsat)
131         pqs(i,k) = zqsat
132      ENDDO
133      ENDDO
134      DO i = 1, klon
135         paprs(i,klev+1) = pres_h(i,1)
136         zgeom(i,klev) = RD * pt(i,klev)
137     .                   / (0.5*(paprs(i,klev+1)+paprsf(i,klev)))
138     .                   * (paprs(i,klev+1)-paprsf(i,klev))
139      ENDDO
140      DO k = klev-1, 1, -1
141      DO i = 1, klon
142         zgeom(i,k) = zgeom(i,k+1)
143     .              + RD * 0.5*(pt(i,k+1)+pt(i,k)) / paprs(i,k+1)
144     .                   * (paprsf(i,k+1)-paprsf(i,k))
145      ENDDO
146      ENDDO
147c
148c appeler la routine principale
149c
150      CALL flxmain(dtime, pt, pq, pqs, pqhfl,
151     .             paprsf, paprs, zgeom, land, zcvgt, zcvgq, pvervel,
152     .             rain, snow, kcbot, kctop, kdtop,
153     .             zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,
154     .             d_t_bis, d_q_bis, zmflxr, zmflxs)
155C
156cAA--------------------------------------------------------
157cAA rem : De la meme facon que l'on effectue le reindicage
158cAA       pour la temperature t et le champ q
159cAA       on reindice les flux necessaires a la convection
160cAA       des traceurs
161cAA--------------------------------------------------------
162      DO k = 1, klev
163      DO i = 1, klon
164         d_q(i,klev+1-k) = dtime*d_q_bis(i,k)
165         d_t(i,klev+1-k) = dtime*d_t_bis(i,k)
166      ENDDO
167      ENDDO
168c
169      DO i = 1, klon
170         pmfu(i,1)= 0.
171         pmfd(i,1)= 0.
172         pen_d(i,1)= 0.
173         pde_d(i,1)= 0.
174      ENDDO
175     
176      DO k = 2, klev
177      DO i = 1, klon
178         pmfu(i,klev+2-k)= zmfu(i,k)
179         pmfd(i,klev+2-k)= zmfd(i,k)
180      ENDDO
181      ENDDO
182c
183      DO k = 1, klev
184      DO i = 1, klon
185         pen_u(i,klev+1-k)=  zen_u(i,k)
186         pde_u(i,klev+1-k)=  zde_u(i,k)
187      ENDDO
188      ENDDO
189c
190      DO k = 1, klev-1
191      DO i = 1, klon
192         pen_d(i,klev+1-k)= -zen_d(i,k+1)
193         pde_d(i,klev+1-k)= -zde_d(i,k+1)
194      ENDDO
195      ENDDO
196
197      DO k = 1, klev+1
198      DO i = 1, klon
199         pmflxr(i,klev+2-k)= zmflxr(i,k)
200         pmflxs(i,klev+2-k)= zmflxs(i,k)
201      ENDDO
202      ENDDO
203
204      RETURN
205      END
206c--------------------------------------------------------------------
207      SUBROUTINE flxmain(pdtime, pten, pqen, pqsen, pqhfl, pap, paph,
208     .                   pgeo, ldland, ptte, pqte, pvervel,
209     .                   prsfc, pssfc, kcbot, kctop, kdtop,
210c     *                   ldcum, ktype,
211     .                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
212     .                   dt_con, dq_con, pmflxr, pmflxs)
213      IMPLICIT none
214C     ------------------------------------------------------------------
215#include "dimensions.h"
216#include "dimphy.h"
217#include "YOMCST.h"
218#include "YOETHF.h"
219#include "YOECUMF.h"
220C     ----------------------------------------------------------------
221      REAL pten(klon,klev), pqen(klon,klev), pqsen(klon,klev)
222      REAL ptte(klon,klev)
223      REAL pqte(klon,klev)
224      REAL pvervel(klon,klev)
225      REAL pgeo(klon,klev), pap(klon,klev), paph(klon,klev+1)
226      REAL pqhfl(klon)
227c
228      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
229      REAL plude(klon,klev)
230      REAL pmfu(klon,klev)
231      REAL prsfc(klon), pssfc(klon)
232      INTEGER  kcbot(klon), kctop(klon), ktype(klon)
233      LOGICAL  ldland(klon), ldcum(klon)
234c
235      REAL ztenh(klon,klev), zqenh(klon,klev), zqsenh(klon,klev)
236      REAL zgeoh(klon,klev)
237      REAL zmfub(klon), zmfub1(klon)
238      REAL zmfus(klon,klev), zmfuq(klon,klev), zmful(klon,klev)
239      REAL zdmfup(klon,klev), zdpmel(klon,klev)
240      REAL zentr(klon), zhcbase(klon)
241      REAL zdqpbl(klon), zdqcv(klon), zdhpbl(klon)
242      REAL zrfl(klon)
243      REAL pmflxr(klon,klev+1)
244      REAL pmflxs(klon,klev+1)
245      INTEGER  ilab(klon,klev), ictop0(klon)
246      LOGICAL  llo1
247      REAL dt_con(klon,klev), dq_con(klon,klev)
248      REAL zmfmax, zdh
249      REAL pdtime, zqumqe, zdqmin, zalvdcp, zhsat, zzz
250      REAL zhhat, zpbmpt, zgam, zeps, zfac
251      INTEGER i, k, ikb, itopm2, kcum
252c
253      REAL pen_u(klon,klev), pde_u(klon,klev)
254      REAL pen_d(klon,klev), pde_d(klon,klev)
255c
256      REAL ptd(klon,klev), pqd(klon,klev), pmfd(klon,klev)
257      REAL zmfds(klon,klev), zmfdq(klon,klev), zdmfdp(klon,klev)
258      INTEGER kdtop(klon)
259      LOGICAL lddraf(klon)
260C---------------------------------------------------------------------
261      LOGICAL firstcal
262      SAVE firstcal
263      DATA firstcal / .TRUE. /
264C---------------------------------------------------------------------
265      IF (firstcal) THEN
266         CALL flxsetup
267         firstcal = .FALSE.
268      ENDIF
269C---------------------------------------------------------------------
270      DO i = 1, klon
271         ldcum(i) = .FALSE.
272      ENDDO
273      DO k = 1, klev
274      DO i = 1, klon
275         dt_con(i,k) = 0.0
276         dq_con(i,k) = 0.0
277      ENDDO
278      ENDDO
279c----------------------------------------------------------------------
280c initialiser les variables et faire l'interpolation verticale
281c----------------------------------------------------------------------
282      CALL flxini(pten, pqen, pqsen, pgeo,
283     .     paph, zgeoh, ztenh, zqenh, zqsenh,
284     .     ptu, pqu, ptd, pqd, pmfd, zmfds, zmfdq, zdmfdp,
285     .     pmfu, zmfus, zmfuq, zdmfup,
286     .     zdpmel, plu, plude, ilab, pen_u, pde_u, pen_d, pde_d)
287c---------------------------------------------------------------------
288c determiner les valeurs au niveau de base de la tour convective
289c---------------------------------------------------------------------
290      CALL flxbase(ztenh, zqenh, zgeoh, paph,
291     *            ptu, pqu, plu, ldcum, kcbot, ilab)
292c---------------------------------------------------------------------
293c calculer la convergence totale de l'humidite et celle en provenance
294c de la couche limite, plus precisement, la convergence integree entre
295c le sol et la base de la convection. Cette derniere convergence est
296c comparee avec l'evaporation obtenue dans la couche limite pour
297c determiner le type de la convection
298c---------------------------------------------------------------------
299      k=1
300      DO i = 1, klon
301         zdqcv(i) = pqte(i,k)*(paph(i,k+1)-paph(i,k))
302         zdhpbl(i) = 0.0
303         zdqpbl(i) = 0.0
304      ENDDO
305c
306      DO k=2,klev
307      DO i = 1, klon
308          zdqcv(i)=zdqcv(i)+pqte(i,k)*(paph(i,k+1)-paph(i,k))
309          IF (k.GE.kcbot(i)) THEN
310             zdqpbl(i)=zdqpbl(i)+pqte(i,k)*(paph(i,k+1)-paph(i,k))
311             zdhpbl(i)=zdhpbl(i)+(RCPD*ptte(i,k)+RLVTT*pqte(i,k))
312     .                          *(paph(i,k+1)-paph(i,k))
313          ENDIF
314      ENDDO
315      ENDDO
316c
317      DO i = 1, klon
318         ktype(i) = 2
319         if (zdqcv(i).GT.MAX(0.,-1.5*pqhfl(i)*RG)) ktype(i) = 1
320ccc         if (zdqcv(i).GT.MAX(0.,-1.1*pqhfl(i)*RG)) ktype(i) = 1
321      ENDDO
322c
323c---------------------------------------------------------------------
324c determiner le flux de masse entrant a travers la base.
325c on ignore, pour l'instant, l'effet du panache descendant
326c---------------------------------------------------------------------
327      DO i = 1, klon
328         ikb=kcbot(i)
329         zqumqe=pqu(i,ikb)+plu(i,ikb)-zqenh(i,ikb)
330         zdqmin=MAX(0.01*zqenh(i,ikb),1.E-10)
331         IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)) THEN
332            zmfub(i) = zdqpbl(i)/(RG*MAX(zqumqe,zdqmin))
333         ELSE
334            zmfub(i) = 0.01
335            ldcum(i)=.FALSE.
336         ENDIF
337         IF (ktype(i).EQ.2) THEN
338            zdh = RCPD*(ptu(i,ikb)-ztenh(i,ikb)) + RLVTT*zqumqe
339            zdh = RG * MAX(zdh,1.0E5*zdqmin)
340            IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub(i)=zdhpbl(i)/zdh
341         ENDIF
342         zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)
343         zmfub(i) = MIN(zmfub(i),zmfmax)
344         zentr(i) = ENTRSCV
345         IF (ktype(i).EQ.1) zentr(i) = ENTRPEN
346      ENDDO
347C-----------------------------------------------------------------------
348C DETERMINE CLOUD ASCENT FOR ENTRAINING PLUME
349C-----------------------------------------------------------------------
350c (A) calculer d'abord la hauteur "theorique" de la tour convective sans
351c     considerer l'entrainement ni le detrainement du panache, sachant
352c     ces derniers peuvent abaisser la hauteur theorique.
353c
354      DO i = 1, klon
355         ikb=kcbot(i)
356         zhcbase(i)=RCPD*ptu(i,ikb)+zgeoh(i,ikb)+RLVTT*pqu(i,ikb)
357         ictop0(i)=kcbot(i)-1
358      ENDDO
359c
360      zalvdcp=RLVTT/RCPD
361      DO k=klev-1,3,-1
362      DO i = 1, klon
363         zhsat=RCPD*ztenh(i,k)+zgeoh(i,k)+RLVTT*zqsenh(i,k)
364         zgam=R5LES*zalvdcp*zqsenh(i,k)/
365     .        ((1.-RETV  *zqsenh(i,k))*(ztenh(i,k)-R4LES)**2)
366         zzz=RCPD*ztenh(i,k)*0.608
367         zhhat=zhsat-(zzz+zgam*zzz)/(1.+zgam*zzz/RLVTT)*
368     .               MAX(zqsenh(i,k)-zqenh(i,k),0.)
369         IF(k.LT.ictop0(i).AND.zhcbase(i).GT.zhhat) ictop0(i)=k
370      ENDDO
371      ENDDO
372c
373c (B) calculer le panache ascendant
374c
375      CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,
376     .     pgeo, zgeoh, pap, paph, pqte, pvervel,
377     .     ldland, ldcum, ktype, ilab,
378     .     ptu, pqu, plu, pmfu, zmfub, zentr,
379     .     zmfus, zmfuq, zmful, plude, zdmfup,
380     .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)
381      IF (kcum.EQ.0) GO TO 1000
382C
383C verifier l'epaisseur de la convection et changer eventuellement
384c le taux d'entrainement/detrainement
385C
386      DO i = 1, klon
387         zpbmpt=paph(i,kcbot(i))-paph(i,kctop(i))
388         IF(ldcum(i).AND.ktype(i).EQ.1.AND.zpbmpt.LT.2.E4)ktype(i)=2
389         IF(ldcum(i)) ictop0(i)=kctop(i)
390         IF(ktype(i).EQ.2) zentr(i)=ENTRSCV
391      ENDDO
392c
393      IF (lmfdd) THEN  ! si l'on considere le panache descendant
394c
395c calculer la precipitation issue du panache ascendant pour
396c determiner l'existence du panache descendant dans la convection
397      DO i = 1, klon
398         zrfl(i)=zdmfup(i,1)
399      ENDDO
400      DO k=2,klev
401      DO i = 1, klon
402         zrfl(i)=zrfl(i)+zdmfup(i,k)
403      ENDDO
404      ENDDO
405c
406c determiner le LFS (level of free sinking: niveau de plonge libre)
407      CALL flxdlfs(ztenh, zqenh, zgeoh, paph, ptu, pqu,
408     *     ldcum,    kcbot,    kctop,    zmfub,    zrfl,
409     *     ptd,      pqd,
410     *     pmfd,     zmfds,    zmfdq,    zdmfdp,
411     *     kdtop,    lddraf)
412c
413c calculer le panache descendant
414      CALL flxddraf(ztenh,    zqenh,
415     *     zgeoh,    paph,     zrfl,
416     *     ptd,      pqd,
417     *     pmfd,     zmfds,    zmfdq,    zdmfdp,
418     *     lddraf, pen_d, pde_d)
419c
420c calculer de nouveau le flux de masse entrant a travers la base
421c de la convection, sachant qu'il a ete modifie par le panache
422c descendant
423      DO i = 1, klon
424      IF (lddraf(i)) THEN
425         ikb = kcbot(i)
426         llo1 = PMFD(i,ikb).LT.0.
427         zeps = 0.
428         IF ( llo1 ) zeps = CMFDEPS
429         zqumqe = pqu(i,ikb)+plu(i,ikb)-
430     .            zeps*pqd(i,ikb)-(1.-zeps)*zqenh(i,ikb)
431         zdqmin = MAX(0.01*zqenh(i,ikb),1.E-10)
432         zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)
433         IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)
434     .       .AND.zmfub(i).LT.zmfmax) THEN
435            zmfub1(i) = zdqpbl(i) / (RG*MAX(zqumqe,zdqmin))
436         ELSE
437            zmfub1(i) = zmfub(i)
438         ENDIF
439         IF (ktype(i).EQ.2) THEN
440            zdh = RCPD*(ptu(i,ikb)-zeps*ptd(i,ikb)-
441     .            (1.-zeps)*ztenh(i,ikb))+RLVTT*zqumqe
442            zdh = RG * MAX(zdh,1.0E5*zdqmin)
443            IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub1(i)=zdhpbl(i)/zdh
444         ENDIF
445         IF ( .NOT.((ktype(i).EQ.1.OR.ktype(i).EQ.2).AND.
446     .              ABS(zmfub1(i)-zmfub(i)).LT.0.2*zmfub(i)) )
447     .      zmfub1(i) = zmfub(i)
448      ENDIF
449      ENDDO
450      DO k = 1, klev
451      DO i = 1, klon
452      IF (lddraf(i)) THEN
453         zfac = zmfub1(i)/MAX(zmfub(i),1.E-10)
454         pmfd(i,k) = pmfd(i,k)*zfac
455         zmfds(i,k) = zmfds(i,k)*zfac
456         zmfdq(i,k) = zmfdq(i,k)*zfac
457         zdmfdp(i,k) = zdmfdp(i,k)*zfac
458         pen_d(i,k) = pen_d(i,k)*zfac
459         pde_d(i,k) = pde_d(i,k)*zfac
460      ENDIF
461      ENDDO
462      ENDDO
463      DO i = 1, klon
464         IF (lddraf(i)) zmfub(i)=zmfub1(i)
465      ENDDO
466c
467      ENDIF   ! fin de test sur lmfdd
468c
469c-----------------------------------------------------------------------
470c calculer de nouveau le panache ascendant
471c-----------------------------------------------------------------------
472      CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,
473     .     pgeo, zgeoh, pap, paph, pqte, pvervel,
474     .     ldland, ldcum, ktype, ilab,
475     .     ptu, pqu, plu, pmfu, zmfub, zentr,
476     .     zmfus, zmfuq, zmful, plude, zdmfup,
477     .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)
478c
479c-----------------------------------------------------------------------
480c determiner les flux convectifs en forme finale, ainsi que
481c la quantite des precipitations
482c-----------------------------------------------------------------------
483      CALL flxflux(pdtime, pqen, pqsen, ztenh, zqenh, pap, paph,
484     .     ldland, zgeoh, kcbot, kctop, lddraf, kdtop, ktype, ldcum,
485     .     pmfu, pmfd, zmfus, zmfds, zmfuq, zmfdq, zmful, plude,
486     .     zdmfup, zdmfdp, pten, prsfc, pssfc, zdpmel, itopm2,
487     .     pmflxr, pmflxs)
488c
489c----------------------------------------------------------------------
490c calculer les tendances pour T et Q
491c----------------------------------------------------------------------
492      CALL flxdtdq(pdtime, itopm2, paph, ldcum, pten,
493     e     zmfus, zmfds, zmfuq, zmfdq, zmful, zdmfup, zdmfdp, zdpmel,
494     s     dt_con,dq_con)
495c
496 1000 CONTINUE
497      RETURN
498      END
499      SUBROUTINE flxini(pten, pqen, pqsen, pgeo, paph, pgeoh, ptenh,
500     .           pqenh, pqsenh, ptu, pqu, ptd, pqd, pmfd, pmfds, pmfdq,
501     .           pdmfdp, pmfu, pmfus, pmfuq, pdmfup, pdpmel, plu, plude,
502     .           klab,pen_u, pde_u, pen_d, pde_d)
503      IMPLICIT none
504C----------------------------------------------------------------------
505C THIS ROUTINE INTERPOLATES LARGE-SCALE FIELDS OF T,Q ETC.
506C TO HALF LEVELS (I.E. GRID FOR MASSFLUX SCHEME),
507C AND INITIALIZES VALUES FOR UPDRAFTS
508C----------------------------------------------------------------------
509#include "dimensions.h"
510#include "dimphy.h"
511#include "YOMCST.h"
512#include "YOETHF.h"
513C
514      REAL pten(klon,klev)   ! temperature (environnement)
515      REAL pqen(klon,klev)   ! humidite (environnement)
516      REAL pqsen(klon,klev)  ! humidite saturante (environnement)
517      REAL pgeo(klon,klev)   ! geopotentiel (g * metre)
518      REAL pgeoh(klon,klev)  ! geopotentiel aux demi-niveaux
519      REAL paph(klon,klev+1) ! pression aux demi-niveaux
520      REAL ptenh(klon,klev)  ! temperature aux demi-niveaux
521      REAL pqenh(klon,klev)  ! humidite aux demi-niveaux
522      REAL pqsenh(klon,klev) ! humidite saturante aux demi-niveaux
523C
524      REAL ptu(klon,klev)    ! temperature du panache ascendant (p-a)
525      REAL pqu(klon,klev)    ! humidite du p-a
526      REAL plu(klon,klev)    ! eau liquide du p-a
527      REAL pmfu(klon,klev)   ! flux de masse du p-a
528      REAL pmfus(klon,klev)  ! flux de l'energie seche dans le p-a
529      REAL pmfuq(klon,klev)  ! flux de l'humidite dans le p-a
530      REAL pdmfup(klon,klev) ! quantite de l'eau precipitee dans p-a
531      REAL plude(klon,klev)  ! quantite de l'eau liquide jetee du
532c                              p-a a l'environnement
533      REAL pdpmel(klon,klev) ! quantite de neige fondue
534c
535      REAL ptd(klon,klev)    ! temperature du panache descendant (p-d)
536      REAL pqd(klon,klev)    ! humidite du p-d
537      REAL pmfd(klon,klev)   ! flux de masse du p-d
538      REAL pmfds(klon,klev)  ! flux de l'energie seche dans le p-d
539      REAL pmfdq(klon,klev)  ! flux de l'humidite dans le p-d
540      REAL pdmfdp(klon,klev) ! quantite de precipitation dans p-d
541c
542      REAL pen_u(klon,klev) ! quantite de masse entrainee pour p-a
543      REAL pde_u(klon,klev) ! quantite de masse detrainee pour p-a
544      REAL pen_d(klon,klev) ! quantite de masse entrainee pour p-d
545      REAL pde_d(klon,klev) ! quantite de masse detrainee pour p-d
546C
547      INTEGER  klab(klon,klev)
548      LOGICAL  llflag(klon)
549      INTEGER k, i, icall
550      REAL zzs
551C----------------------------------------------------------------------
552C SPECIFY LARGE SCALE PARAMETERS AT HALF LEVELS
553C ADJUST TEMPERATURE FIELDS IF STATICLY UNSTABLE
554C----------------------------------------------------------------------
555      DO 130 k = 2, klev
556c
557      DO i = 1, klon
558         pgeoh(i,k)=pgeo(i,k)+(pgeo(i,k-1)-pgeo(i,k))*0.5
559         ptenh(i,k)=(MAX(RCPD*pten(i,k-1)+pgeo(i,k-1),
560     .             RCPD*pten(i,k)+pgeo(i,k))-pgeoh(i,k))/RCPD
561         pqsenh(i,k)=pqsen(i,k-1)
562         llflag(i)=.TRUE.
563      ENDDO
564c
565      icall=0
566      CALL flxadjtq(paph(1,k),ptenh(1,k),pqsenh(1,k),llflag,icall)
567c
568      DO i = 1, klon
569         pqenh(i,k)=MIN(pqen(i,k-1),pqsen(i,k-1))
570     .               +(pqsenh(i,k)-pqsen(i,k-1))
571         pqenh(i,k)=MAX(pqenh(i,k),0.)
572      ENDDO
573c
574  130 CONTINUE
575C
576      DO 140 i = 1, klon
577         ptenh(i,klev)=(RCPD*pten(i,klev)+pgeo(i,klev)-
578     1                   pgeoh(i,klev))/RCPD
579         pqenh(i,klev)=pqen(i,klev)
580         ptenh(i,1)=pten(i,1)
581         pqenh(i,1)=pqen(i,1)
582         pgeoh(i,1)=pgeo(i,1)
583  140 CONTINUE
584c
585      DO 160 k = klev-1, 2, -1
586      DO 150 i = 1, klon
587         zzs = MAX(RCPD*ptenh(i,k)+pgeoh(i,k),
588     .             RCPD*ptenh(i,k+1)+pgeoh(i,k+1))
589         ptenh(i,k) = (zzs-pgeoh(i,k))/RCPD
590  150 CONTINUE
591  160 CONTINUE
592C
593C-----------------------------------------------------------------------
594C INITIALIZE VALUES FOR UPDRAFTS AND DOWNDRAFTS
595C-----------------------------------------------------------------------
596      DO k = 1, klev
597      DO i = 1, klon
598         ptu(i,k) = ptenh(i,k)
599         pqu(i,k) = pqenh(i,k)
600         plu(i,k) = 0.
601         pmfu(i,k) = 0.
602         pmfus(i,k) = 0.
603         pmfuq(i,k) = 0.
604         pdmfup(i,k) = 0.
605         pdpmel(i,k) = 0.
606         plude(i,k) = 0.
607c
608         klab(i,k) = 0
609c
610         ptd(i,k) = ptenh(i,k)
611         pqd(i,k) = pqenh(i,k)
612         pmfd(i,k) = 0.0
613         pmfds(i,k) = 0.0
614         pmfdq(i,k) = 0.0
615         pdmfdp(i,k) = 0.0
616c
617         pen_u(i,k) = 0.0
618         pde_u(i,k) = 0.0
619         pen_d(i,k) = 0.0
620         pde_d(i,k) = 0.0
621      ENDDO
622      ENDDO
623C
624      RETURN
625      END
626      SUBROUTINE flxbase(ptenh, pqenh, pgeoh, paph,
627     *     ptu, pqu, plu, ldcum, kcbot, klab)
628      IMPLICIT none
629C----------------------------------------------------------------------
630C THIS ROUTINE CALCULATES CLOUD BASE VALUES (T AND Q)
631C
632C INPUT ARE ENVIRONM. VALUES OF T,Q,P,PHI AT HALF LEVELS.
633C IT RETURNS CLOUD BASE VALUES AND FLAGS AS FOLLOWS;
634C   klab=1 FOR SUBCLOUD LEVELS
635C   klab=2 FOR CONDENSATION LEVEL
636C
637C LIFT SURFACE AIR DRY-ADIABATICALLY TO CLOUD BASE
638C (NON ENTRAINING PLUME,I.E.CONSTANT MASSFLUX)
639C----------------------------------------------------------------------
640#include "dimensions.h"
641#include "dimphy.h"
642#include "YOMCST.h"
643#include "YOETHF.h"
644C       ----------------------------------------------------------------
645      REAL ptenh(klon,klev), pqenh(klon,klev)
646      REAL pgeoh(klon,klev), paph(klon,klev+1)
647C
648      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
649      INTEGER  klab(klon,klev), kcbot(klon)
650C
651      LOGICAL llflag(klon), ldcum(klon)
652      INTEGER i, k, icall, is
653      REAL zbuo, zqold(klon)
654C----------------------------------------------------------------------
655C INITIALIZE VALUES AT LIFTING LEVEL
656C----------------------------------------------------------------------
657      DO i = 1, klon
658         klab(i,klev)=1
659         kcbot(i)=klev-1
660         ldcum(i)=.FALSE.
661      ENDDO
662C----------------------------------------------------------------------
663C DO ASCENT IN SUBCLOUD LAYER,
664C CHECK FOR EXISTENCE OF CONDENSATION LEVEL,
665C ADJUST T,Q AND L ACCORDINGLY
666C CHECK FOR BUOYANCY AND SET FLAGS
667C----------------------------------------------------------------------
668      DO 290 k = klev-1, 2, -1
669c
670      is = 0
671      DO i = 1, klon
672         IF (klab(i,k+1).EQ.1) is = is + 1
673         llflag(i) = .FALSE.
674         IF (klab(i,k+1).EQ.1) llflag(i) = .TRUE.
675      ENDDO
676      IF (is.EQ.0) GOTO 290
677c
678      DO i = 1, klon
679      IF(llflag(i)) THEN
680         pqu(i,k) = pqu(i,k+1)
681         ptu(i,k) = ptu(i,k+1)+(pgeoh(i,k+1)-pgeoh(i,k))/RCPD
682         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
683     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5
684         IF (zbuo.GT.0.) klab(i,k) = 1
685         zqold(i) = pqu(i,k)
686      ENDIF
687      ENDDO
688c
689      icall=1
690      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
691c
692      DO i = 1, klon
693      IF (llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN
694         klab(i,k) = 2
695         plu(i,k) = plu(i,k) + zqold(i)-pqu(i,k)
696         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
697     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5
698         IF (zbuo.GT.0.) kcbot(i) = k
699         IF (zbuo.GT.0.) ldcum(i) = .TRUE.
700      ENDIF
701      ENDDO
702c
703  290 CONTINUE
704c
705      RETURN
706      END
707      SUBROUTINE flxasc(pdtime, ptenh, pqenh, pten, pqen, pqsen,
708     .     pgeo, pgeoh, pap, paph, pqte, pvervel,
709     .     ldland, ldcum, ktype, klab, ptu, pqu, plu,
710     .     pmfu, pmfub, pentr, pmfus, pmfuq,
711     .     pmful, plude, pdmfup, kcbot, kctop, kctop0, kcum,
712     .     pen_u, pde_u)
713      IMPLICIT none
714C----------------------------------------------------------------------
715C THIS ROUTINE DOES THE CALCULATIONS FOR CLOUD ASCENTS
716C FOR CUMULUS PARAMETERIZATION
717C----------------------------------------------------------------------
718#include "dimensions.h"
719#include "dimphy.h"
720#include "YOMCST.h"
721#include "YOETHF.h"
722#include "YOECUMF.h"
723C
724      REAL pdtime
725      REAL pten(klon,klev), ptenh(klon,klev)
726      REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)
727      REAL pgeo(klon,klev), pgeoh(klon,klev)
728      REAL pap(klon,klev), paph(klon,klev+1)
729      REAL pqte(klon,klev)
730      REAL pvervel(klon,klev) ! vitesse verticale en Pa/s
731C
732      REAL pmfub(klon), pentr(klon)
733      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
734      REAL plude(klon,klev)
735      REAL pmfu(klon,klev), pmfus(klon,klev)
736      REAL pmfuq(klon,klev), pmful(klon,klev)
737      REAL pdmfup(klon,klev)
738      INTEGER ktype(klon), klab(klon,klev), kcbot(klon), kctop(klon)
739      INTEGER kctop0(klon)
740      LOGICAL ldland(klon), ldcum(klon)
741C
742      REAL pen_u(klon,klev), pde_u(klon,klev)
743      REAL zqold(klon)
744      REAL zdland(klon)
745      LOGICAL llflag(klon)
746      INTEGER k, i, is, icall, kcum
747      REAL ztglace, zdphi, zqeen, zseen, zscde, zqude
748      REAL zmfusk, zmfuqk, zmfulk, zbuo, zdnoprc, zprcon, zlnew
749c
750      REAL zpbot(klon), zptop(klon), zrho(klon)
751      REAL zdprho, zentr, zpmid, zmftest, zmfmax
752      LOGICAL llo1, llo2
753c
754      REAL zwmax(klon), zzzmb
755      INTEGER klwmin(klon) ! level of maximum vertical velocity
756C----------------------------------------------------------------------
757      ztglace = RTT - 13.
758c
759c Chercher le niveau ou la vitesse verticale est maximale:
760      DO i = 1, klon
761         klwmin(i) = klev
762         zwmax(i) = 0.0
763      ENDDO
764      DO k = klev, 3, -1
765      DO i = 1, klon
766      IF (pvervel(i,k).LT.zwmax(i)) THEN
767         zwmax(i) = pvervel(i,k)
768         klwmin(i) = k
769      ENDIF
770      ENDDO
771      ENDDO
772C----------------------------------------------------------------------
773C SET DEFAULT VALUES
774C----------------------------------------------------------------------
775      DO i = 1, klon
776         IF (.NOT.ldcum(i)) ktype(i)=0
777      ENDDO
778c
779      DO k=1,klev
780      DO i = 1, klon
781         plu(i,k)=0.
782         pmfu(i,k)=0.
783         pmfus(i,k)=0.
784         pmfuq(i,k)=0.
785         pmful(i,k)=0.
786         plude(i,k)=0.
787         pdmfup(i,k)=0.
788         IF(.NOT.ldcum(i).OR.ktype(i).EQ.3) klab(i,k)=0
789         IF(.NOT.ldcum(i).AND.paph(i,k).LT.4.E4) kctop0(i)=k
790      ENDDO
791      ENDDO
792c
793      DO i = 1, klon
794      IF (ldland(i)) THEN
795         zdland(i)=3.0E4
796         zdphi=pgeoh(i,kctop0(i))-pgeoh(i,kcbot(i))
797         IF (ptu(i,kctop0(i)).GE.ztglace) zdland(i)=zdphi
798         zdland(i)=MAX(3.0E4,zdland(i))
799         zdland(i)=MIN(5.0E4,zdland(i))
800      ENDIF
801      ENDDO
802C
803C Initialiser les valeurs au niveau d'ascendance
804C
805      DO i = 1, klon
806         kctop(i) = klev-1
807         IF (.NOT.ldcum(i)) THEN
808            kcbot(i) = klev-1
809            pmfub(i) = 0.
810            pqu(i,klev) = 0.
811         ENDIF
812         pmfu(i,klev) = pmfub(i)
813         pmfus(i,klev) = pmfub(i)*(RCPD*ptu(i,klev)+pgeoh(i,klev))
814         pmfuq(i,klev) = pmfub(i)*pqu(i,klev)
815      ENDDO
816c
817      DO i = 1, klon
818         ldcum(i) = .FALSE.
819      ENDDO
820C----------------------------------------------------------------------
821C  DO ASCENT: SUBCLOUD LAYER (klab=1) ,CLOUDS (klab=2)
822C  BY DOING FIRST DRY-ADIABATIC ASCENT AND THEN
823C  BY ADJUSTING T,Q AND L ACCORDINGLY IN *flxadjtq*,
824C  THEN CHECK FOR BUOYANCY AND SET FLAGS ACCORDINGLY
825C----------------------------------------------------------------------
826      DO 480 k = klev-1,3,-1
827c
828      IF (LMFMID .AND. k.LT.klev-1 .AND. k.GT.klev/2) THEN
829         DO i = 1, klon
830         IF (.NOT.ldcum(i) .AND. klab(i,k+1).EQ.0 .AND.
831     .       pqen(i,k).GT.0.9*pqsen(i,k)) THEN
832            ptu(i,k+1) = pten(i,k) +(pgeo(i,k)-pgeoh(i,k+1))/RCPD
833            pqu(i,k+1) = pqen(i,k)
834            plu(i,k+1) = 0.0
835            zzzmb = MAX(CMFCMIN, -pvervel(i,k))
836            zmfmax = (paph(i,k)-paph(i,k-1))/(RG*pdtime)
837            pmfub(i) = MIN(zzzmb,zmfmax)
838            pmfu(i,k+1) = pmfub(i)
839            pmfus(i,k+1) = pmfub(i)*(RCPD*ptu(i,k+1)+pgeoh(i,k+1))
840            pmfuq(i,k+1) = pmfub(i)*pqu(i,k+1)
841            pmful(i,k+1) = 0.0
842            pdmfup(i,k+1) = 0.0
843            kcbot(i) = k
844            klab(i,k+1) = 1
845            ktype(i) = 3
846            pentr(i) = ENTRMID
847         ENDIF
848         ENDDO
849      ENDIF
850c
851      is = 0
852      DO i = 1, klon
853         is = is + klab(i,k+1)
854         IF (klab(i,k+1) .EQ. 0) klab(i,k) = 0
855         llflag(i) = .FALSE.
856         IF (klab(i,k+1) .GT. 0) llflag(i) = .TRUE.
857      ENDDO
858      IF (is .EQ. 0) GOTO 480
859c
860c calculer le taux d'entrainement et de detrainement
861c
862      DO i = 1, klon
863         pen_u(i,k) = 0.0
864         pde_u(i,k) = 0.0
865         zrho(i)=paph(i,k+1)/(RD*ptenh(i,k+1))
866         zpbot(i)=paph(i,kcbot(i))
867         zptop(i)=paph(i,kctop0(i))
868      ENDDO
869c
870      DO 125 i = 1, klon
871      IF(ldcum(i)) THEN
872         zdprho=(paph(i,k+1)-paph(i,k))/(RG*zrho(i))
873         zentr=pentr(i)*pmfu(i,k+1)*zdprho
874         llo1=k.LT.kcbot(i)
875         IF(llo1) pde_u(i,k)=zentr
876         zpmid=0.5*(zpbot(i)+zptop(i))
877         llo2=llo1.AND.ktype(i).EQ.2.AND.
878     .        (zpbot(i)-paph(i,k).LT.0.2E5.OR.
879     .         paph(i,k).GT.zpmid)
880         IF(llo2) pen_u(i,k)=zentr
881         llo2=llo1.AND.(ktype(i).EQ.1.OR.ktype(i).EQ.3).AND.
882     .        (k.GE.MAX(klwmin(i),kctop0(i)+2).OR.pap(i,k).GT.zpmid)
883         IF(llo2) pen_u(i,k)=zentr
884         llo1=pen_u(i,k).GT.0..AND.(ktype(i).EQ.1.OR.ktype(i).EQ.2)
885         IF(llo1) THEN
886            zentr=zentr*(1.+3.*(1.-MIN(1.,(zpbot(i)-pap(i,k))/1.5E4)))
887            pen_u(i,k)=pen_u(i,k)*(1.+3.*(1.-MIN(1.,
888     .                 (zpbot(i)-pap(i,k))/1.5E4)))
889            pde_u(i,k)=pde_u(i,k)*(1.+3.*(1.-MIN(1.,
890     .                 (zpbot(i)-pap(i,k))/1.5E4)))
891         ENDIF
892         IF(llo2.AND.pqenh(i,k+1).GT.1.E-5)
893     .   pen_u(i,k)=zentr+MAX(pqte(i,k),0.)/pqenh(i,k+1)*
894     .              zrho(i)*zdprho
895      ENDIF
896  125 CONTINUE
897c
898C----------------------------------------------------------------------
899c DO ADIABATIC ASCENT FOR ENTRAINING/DETRAINING PLUME
900C----------------------------------------------------------------------
901c
902      DO 420 i = 1, klon
903      IF (llflag(i)) THEN
904         IF (k.LT.kcbot(i)) THEN
905            zmftest = pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
906            zmfmax = MIN(zmftest,(paph(i,k)-paph(i,k-1))/(RG*pdtime))
907            pen_u(i,k)=MAX(pen_u(i,k)-MAX(0.0,zmftest-zmfmax),0.0)
908         ENDIF
909         pde_u(i,k)=MIN(pde_u(i,k),0.75*pmfu(i,k+1))
910c calculer le flux de masse du niveau k a partir de celui du k+1
911         pmfu(i,k)=pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
912c calculer les valeurs Su, Qu et l du niveau k dans le panache montant
913         zqeen=pqenh(i,k+1)*pen_u(i,k)
914         zseen=(RCPD*ptenh(i,k+1)+pgeoh(i,k+1))*pen_u(i,k)
915         zscde=(RCPD*ptu(i,k+1)+pgeoh(i,k+1))*pde_u(i,k)
916         zqude=pqu(i,k+1)*pde_u(i,k)
917         plude(i,k)=plu(i,k+1)*pde_u(i,k)
918         zmfusk=pmfus(i,k+1)+zseen-zscde
919         zmfuqk=pmfuq(i,k+1)+zqeen-zqude
920         zmfulk=pmful(i,k+1)    -plude(i,k)
921         plu(i,k)=zmfulk*(1./MAX(CMFCMIN,pmfu(i,k)))
922         pqu(i,k)=zmfuqk*(1./MAX(CMFCMIN,pmfu(i,k)))
923         ptu(i,k)=(zmfusk*(1./MAX(CMFCMIN,pmfu(i,k)))-
924     1               pgeoh(i,k))/RCPD
925         ptu(i,k)=MAX(100.,ptu(i,k))
926         ptu(i,k)=MIN(400.,ptu(i,k))
927         zqold(i)=pqu(i,k)
928      ELSE
929         zqold(i)=0.0
930      ENDIF
931  420 CONTINUE
932c
933C----------------------------------------------------------------------
934c DO CORRECTIONS FOR MOIST ASCENT BY ADJUSTING T,Q AND L
935C----------------------------------------------------------------------
936c
937      icall = 1
938      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
939C
940      DO 440 i = 1, klon
941      IF(llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN
942         klab(i,k) = 2
943         plu(i,k) = plu(i,k)+zqold(i)-pqu(i,k)
944         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
945     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))
946         IF (klab(i,k+1).EQ.1) zbuo=zbuo+0.5
947         IF (zbuo.GT.0..AND.pmfu(i,k).GE.0.1*pmfub(i)) THEN
948            kctop(i) = k
949            ldcum(i) = .TRUE.
950            zdnoprc = 1.5E4
951            IF (ldland(i)) zdnoprc = zdland(i)
952            zprcon = CPRCON
953            IF ((zpbot(i)-paph(i,k)).LT.zdnoprc) zprcon = 0.0
954            zlnew=plu(i,k)/(1.+zprcon*(pgeoh(i,k)-pgeoh(i,k+1)))
955            pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
956            plu(i,k)=zlnew
957         ELSE
958            klab(i,k)=0
959            pmfu(i,k)=0.
960         ENDIF
961      ENDIF
962  440 CONTINUE
963      DO 455 i = 1, klon
964      IF (llflag(i)) THEN
965         pmful(i,k)=plu(i,k)*pmfu(i,k)
966         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
967         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
968      ENDIF
969  455 CONTINUE
970C
971  480 CONTINUE
972C----------------------------------------------------------------------
973C DETERMINE CONVECTIVE FLUXES ABOVE NON-BUOYANCY LEVEL
974C    (NOTE: CLOUD VARIABLES LIKE T,Q AND L ARE NOT
975C           AFFECTED BY DETRAINMENT AND ARE ALREADY KNOWN
976C           FROM PREVIOUS CALCULATIONS ABOVE)
977C----------------------------------------------------------------------
978      DO i = 1, klon
979         IF (kctop(i).EQ.klev-1) ldcum(i) = .FALSE.
980         kcbot(i) = MAX(kcbot(i),kctop(i))
981      ENDDO
982c
983      ldcum(1)=ldcum(1)
984c
985      is = 0
986      DO i = 1, klon
987         if (ldcum(i)) is = is + 1
988      ENDDO
989      kcum = is
990      IF (is.EQ.0) GOTO 800
991c
992      DO 530 i = 1, klon
993      IF (ldcum(i)) THEN
994         k=kctop(i)-1
995         pde_u(i,k)=(1.-CMFCTOP)*pmfu(i,k+1)
996         plude(i,k)=pde_u(i,k)*plu(i,k+1)
997         pmfu(i,k)=pmfu(i,k+1)-pde_u(i,k)
998         zlnew=plu(i,k)
999         pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
1000         plu(i,k)=zlnew
1001         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
1002         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
1003         pmful(i,k)=plu(i,k)*pmfu(i,k)
1004         plude(i,k-1)=pmful(i,k)
1005      ENDIF
1006  530 CONTINUE
1007C
1008  800 CONTINUE
1009      RETURN
1010      END
1011      SUBROUTINE flxflux(pdtime, pqen, pqsen, ptenh, pqenh, pap
1012     .  ,  paph, ldland, pgeoh, kcbot, kctop, lddraf, kdtop
1013     .  ,  ktype, ldcum, pmfu, pmfd, pmfus, pmfds
1014     .  ,  pmfuq, pmfdq, pmful, plude, pdmfup, pdmfdp
1015     .  ,  pten, prfl, psfl, pdpmel, ktopm2
1016     .  ,  pmflxr, pmflxs)
1017      IMPLICIT none
1018C----------------------------------------------------------------------
1019C THIS ROUTINE DOES THE FINAL CALCULATION OF CONVECTIVE
1020C FLUXES IN THE CLOUD LAYER AND IN THE SUBCLOUD LAYER
1021C----------------------------------------------------------------------
1022#include "dimensions.h"
1023#include "dimphy.h"
1024#include "YOMCST.h"
1025#include "YOETHF.h"
1026#include "YOECUMF.h"
1027C
1028      REAL cevapcu(klev)
1029C     -----------------------------------------------------------------
1030      REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)
1031      REAL pten(klon,klev), ptenh(klon,klev)
1032      REAL paph(klon,klev+1), pgeoh(klon,klev)
1033c
1034      REAL pap(klon,klev)
1035      REAL ztmsmlt, zdelta, zqsat
1036C
1037      REAL pmfu(klon,klev), pmfus(klon,klev)
1038      REAL pmfd(klon,klev), pmfds(klon,klev)
1039      REAL pmfuq(klon,klev), pmful(klon,klev)
1040      REAL pmfdq(klon,klev)
1041      REAL plude(klon,klev)
1042      REAL pdmfup(klon,klev), pdpmel(klon,klev)
1043cjq The variable maxpdmfdp(klon) has been introduced by Olivier Boucher
1044cjq 14/11/00 to fix the problem with the negative precipitation.     
1045      REAL pdmfdp(klon,klev), maxpdmfdp(klon,klev)
1046      REAL prfl(klon), psfl(klon)
1047      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
1048      INTEGER  kcbot(klon), kctop(klon), ktype(klon)
1049      LOGICAL  ldland(klon), ldcum(klon)
1050      INTEGER k, kp, i
1051      REAL zcons1, zcons2, zcucov, ztmelp2
1052      REAL pdtime, zdp, zzp, zfac, zsnmlt, zrfl, zrnew
1053      REAL zrmin, zrfln, zdrfl
1054      REAL zpds, zpdr, zdenom
1055      INTEGER ktopm2, itop, ikb
1056c
1057      LOGICAL lddraf(klon)
1058      INTEGER kdtop(klon)
1059c
1060#include "FCTTRE.h"
1061c
1062      DO 101 k=1,klev
1063      CEVAPCU(k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293)
1064     1 *SQRT(0.5*(paph(1,k)+paph(1,k+1))/paph(1,klev+1)) ) * 0.5/RG
1065 101  CONTINUE
1066c
1067c SPECIFY CONSTANTS
1068c
1069      zcons1 = RCPD/(RLMLT*RG*pdtime)
1070      zcons2 = 1./(RG*pdtime)
1071      zcucov = 0.05
1072      ztmelp2 = RTT + 2.
1073c
1074c DETERMINE FINAL CONVECTIVE FLUXES
1075c
1076      itop=klev
1077      DO 110 i = 1, klon
1078         itop=MIN(itop,kctop(i))
1079         IF (.NOT.ldcum(i) .OR. kdtop(i).LT.kctop(i)) lddraf(i)=.FALSE.
1080         IF(.NOT.ldcum(i)) ktype(i)=0
1081  110 CONTINUE
1082c
1083      ktopm2=itop-2
1084      DO 120 k=ktopm2,klev
1085      DO 115 i = 1, klon
1086      IF(ldcum(i).AND.k.GE.kctop(i)-1) THEN
1087         pmfus(i,k)=pmfus(i,k)-pmfu(i,k)*
1088     .                (RCPD*ptenh(i,k)+pgeoh(i,k))
1089         pmfuq(i,k)=pmfuq(i,k)-pmfu(i,k)*pqenh(i,k)
1090         zdp = 1.5E4
1091         IF ( ldland(i) ) zdp = 3.E4
1092c
1093c        l'eau liquide detrainee est precipitee quand certaines
1094c        conditions sont reunies (sinon, elle est consideree
1095c        evaporee dans l'environnement)
1096c
1097         IF(paph(i,kcbot(i))-paph(i,kctop(i)).GE.zdp.AND.
1098     .      pqen(i,k-1).GT.0.8*pqsen(i,k-1))
1099     .      pdmfup(i,k-1)=pdmfup(i,k-1)+plude(i,k-1)
1100c
1101         IF(lddraf(i).AND.k.GE.kdtop(i)) THEN
1102            pmfds(i,k)=pmfds(i,k)-pmfd(i,k)*
1103     .                   (RCPD*ptenh(i,k)+pgeoh(i,k))
1104            pmfdq(i,k)=pmfdq(i,k)-pmfd(i,k)*pqenh(i,k)
1105         ELSE
1106            pmfd(i,k)=0.
1107            pmfds(i,k)=0.
1108            pmfdq(i,k)=0.
1109            pdmfdp(i,k-1)=0.
1110         END IF
1111      ELSE
1112         pmfu(i,k)=0.
1113         pmfus(i,k)=0.
1114         pmfuq(i,k)=0.
1115         pmful(i,k)=0.
1116         pdmfup(i,k-1)=0.
1117         plude(i,k-1)=0.
1118         pmfd(i,k)=0.
1119         pmfds(i,k)=0.
1120         pmfdq(i,k)=0.
1121         pdmfdp(i,k-1)=0.
1122      ENDIF
1123  115 CONTINUE
1124  120 CONTINUE
1125c
1126      DO 130 k=ktopm2,klev
1127      DO 125 i = 1, klon
1128      IF(ldcum(i).AND.k.GT.kcbot(i)) THEN
1129         ikb=kcbot(i)
1130         zzp=((paph(i,klev+1)-paph(i,k))/
1131     .        (paph(i,klev+1)-paph(i,ikb)))
1132         IF (ktype(i).EQ.3) zzp = zzp**2
1133         pmfu(i,k)=pmfu(i,ikb)*zzp
1134         pmfus(i,k)=pmfus(i,ikb)*zzp
1135         pmfuq(i,k)=pmfuq(i,ikb)*zzp
1136         pmful(i,k)=pmful(i,ikb)*zzp
1137      ENDIF
1138  125 CONTINUE
1139  130 CONTINUE
1140c
1141c CALCULATE RAIN/SNOW FALL RATES
1142c CALCULATE MELTING OF SNOW
1143c CALCULATE EVAPORATION OF PRECIP
1144c
1145      DO k = 1, klev+1
1146      DO i = 1, klon
1147         pmflxr(i,k) = 0.0
1148         pmflxs(i,k) = 0.0
1149      ENDDO
1150      ENDDO
1151      DO k = ktopm2, klev
1152      DO i = 1, klon
1153      IF (ldcum(i)) THEN
1154         IF (pmflxs(i,k).GT.0.0 .AND. pten(i,k).GT.ztmelp2) THEN
1155            zfac=zcons1*(paph(i,k+1)-paph(i,k))
1156            zsnmlt=MIN(pmflxs(i,k),zfac*(pten(i,k)-ztmelp2))
1157            pdpmel(i,k)=zsnmlt
1158            ztmsmlt=pten(i,k)-zsnmlt/zfac
1159            zdelta=MAX(0.,SIGN(1.,RTT-ztmsmlt))
1160            zqsat=R2ES*FOEEW(ztmsmlt, zdelta) / pap(i,k)
1161            zqsat=MIN(0.5,zqsat)
1162            zqsat=zqsat/(1.-RETV  *zqsat)
1163            pqsen(i,k) = zqsat
1164         ENDIF
1165         IF (pten(i,k).GT.RTT) THEN
1166         pmflxr(i,k+1)=pmflxr(i,k)+pdmfup(i,k)+pdmfdp(i,k)+pdpmel(i,k)
1167         pmflxs(i,k+1)=pmflxs(i,k)-pdpmel(i,k)
1168         ELSE
1169           pmflxs(i,k+1)=pmflxs(i,k)+pdmfup(i,k)+pdmfdp(i,k)
1170           pmflxr(i,k+1)=pmflxr(i,k)
1171         ENDIF
1172c        si la precipitation est negative, on ajuste le plux du
1173c        panache descendant pour eliminer la negativite
1174         IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN
1175            pdmfdp(i,k) = -pmflxr(i,k)-pmflxs(i,k)-pdmfup(i,k)
1176            pmflxr(i,k+1) = 0.0
1177            pmflxs(i,k+1) = 0.0
1178            pdpmel(i,k) = 0.0
1179         ENDIF
1180      ENDIF
1181      ENDDO
1182      ENDDO
1183c
1184cjq The new variable is initialized here.
1185cjq It contains the humidity which is fed to the downdraft
1186cjq by evaporation of precipitation in the column below the base
1187cjq of convection.
1188cjq
1189cjq In the former version, this term has been subtracted from precip
1190cjq as well as the evaporation.
1191cjq     
1192      DO k = 1, klev
1193      DO i = 1, klon
1194         maxpdmfdp(i,k)=0.0
1195      ENDDO
1196      ENDDO
1197      DO k = 1, klev
1198       DO kp = k, klev
1199        DO i = 1, klon
1200         maxpdmfdp(i,k)=maxpdmfdp(i,k)+pdmfdp(i,kp)
1201        ENDDO
1202       ENDDO
1203      ENDDO
1204cjq End of initialization
1205c     
1206      DO k = ktopm2, klev
1207      DO i = 1, klon
1208      IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN
1209         zrfl = pmflxr(i,k) + pmflxs(i,k)
1210         IF (zrfl.GT.1.0E-20) THEN
1211            zrnew=(MAX(0.,SQRT(zrfl/zcucov)-
1212     .            CEVAPCU(k)*(paph(i,k+1)-paph(i,k))*
1213     .            MAX(0.,pqsen(i,k)-pqen(i,k))))**2*zcucov
1214            zrmin=zrfl-zcucov*MAX(0.,0.8*pqsen(i,k)-pqen(i,k))
1215     .            *zcons2*(paph(i,k+1)-paph(i,k))
1216            zrnew=MAX(zrnew,zrmin)
1217            zrfln=MAX(zrnew,0.)
1218            zdrfl=MIN(0.,zrfln-zrfl)
1219cjq At least the amount of precipiation needed to feed the downdraft
1220cjq with humidity below the base of convection has to be left and can't
1221cjq be evaporated (surely the evaporation can't be positive):           
1222            zdrfl=MAX(zdrfl,
1223     .            MIN(-pmflxr(i,k)-pmflxs(i,k)-maxpdmfdp(i,k),0.0))
1224cjq End of insertion
1225c           
1226            zdenom=1.0/MAX(1.0E-20,pmflxr(i,k)+pmflxs(i,k))
1227            IF (pten(i,k).GT.RTT) THEN
1228               zpdr = pdmfdp(i,k)
1229               zpds = 0.0
1230            ELSE
1231               zpdr = 0.0
1232               zpds = pdmfdp(i,k)
1233            ENDIF
1234            pmflxr(i,k+1) = pmflxr(i,k) + zpdr + pdpmel(i,k)
1235     .                    + zdrfl*pmflxr(i,k)*zdenom
1236            pmflxs(i,k+1) = pmflxs(i,k) + zpds - pdpmel(i,k)
1237     .                    + zdrfl*pmflxs(i,k)*zdenom
1238            pdmfup(i,k) = pdmfup(i,k) + zdrfl
1239         ELSE
1240            pmflxr(i,k+1) = 0.0
1241            pmflxs(i,k+1) = 0.0
1242            pdmfdp(i,k) = 0.0
1243            pdpmel(i,k) = 0.0
1244         ENDIF         
1245         if (pmflxr(i,k) + pmflxs(i,k).lt.-1.e-26)
1246     .    write(*,*) 'precip. < 1e-16 ',pmflxr(i,k) + pmflxs(i,k)
1247      ENDIF
1248      ENDDO
1249      ENDDO
1250c
1251      DO 210 i = 1, klon
1252         prfl(i) = pmflxr(i,klev+1)
1253         psfl(i) = pmflxs(i,klev+1)
1254  210 CONTINUE
1255c
1256      RETURN
1257      END
1258      SUBROUTINE flxdtdq(pdtime, ktopm2, paph, ldcum, pten
1259     .  ,  pmfus, pmfds, pmfuq, pmfdq, pmful, pdmfup, pdmfdp
1260     .  ,  pdpmel, dt_con, dq_con)
1261      IMPLICIT none
1262c----------------------------------------------------------------------
1263c calculer les tendances T et Q
1264c----------------------------------------------------------------------
1265#include "dimensions.h"
1266#include "dimphy.h"
1267#include "YOMCST.h"
1268#include "YOETHF.h"
1269#include "YOECUMF.h"
1270C     -----------------------------------------------------------------
1271      LOGICAL  llo1
1272C
1273      REAL pten(klon,klev), paph(klon,klev+1)
1274      REAL pmfus(klon,klev), pmfuq(klon,klev), pmful(klon,klev)
1275      REAL pmfds(klon,klev), pmfdq(klon,klev)
1276      REAL pdmfup(klon,klev)
1277      REAL pdmfdp(klon,klev)
1278      REAL pdpmel(klon,klev)
1279      LOGICAL ldcum(klon)
1280      REAL dt_con(klon,klev), dq_con(klon,klev)
1281c
1282      INTEGER ktopm2
1283      REAL pdtime
1284c
1285      INTEGER i, k
1286      REAL zalv, zdtdt, zdqdt
1287c
1288      DO 210 k=ktopm2,klev-1
1289      DO 220 i = 1, klon
1290      IF (ldcum(i)) THEN
1291         llo1 = (pten(i,k)-RTT).GT.0.
1292         zalv = RLSTT
1293         IF (llo1) zalv = RLVTT
1294         zdtdt=RG/(paph(i,k+1)-paph(i,k))/RCPD
1295     .        *(pmfus(i,k+1)-pmfus(i,k)
1296     .         +pmfds(i,k+1)-pmfds(i,k)
1297     .          -RLMLT*pdpmel(i,k)
1298     .          -zalv*(pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1299     .         )
1300         dt_con(i,k)=zdtdt
1301         zdqdt=RG/(paph(i,k+1)-paph(i,k))
1302     .        *(pmfuq(i,k+1)-pmfuq(i,k)
1303     .         +pmfdq(i,k+1)-pmfdq(i,k)
1304     .          +pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1305         dq_con(i,k)=zdqdt
1306      ENDIF
1307  220 CONTINUE
1308  210 CONTINUE
1309C
1310      k = klev
1311      DO 230 i = 1, klon
1312      IF (ldcum(i)) THEN
1313         llo1 = (pten(i,k)-RTT).GT.0.
1314         zalv = RLSTT
1315         IF (llo1) zalv = RLVTT
1316         zdtdt=-RG/(paph(i,k+1)-paph(i,k))/RCPD
1317     .         *(pmfus(i,k)+pmfds(i,k)+RLMLT*pdpmel(i,k)
1318     .           -zalv*(pmful(i,k)+pdmfup(i,k)+pdmfdp(i,k)))
1319         dt_con(i,k)=zdtdt
1320         zdqdt=-RG/(paph(i,k+1)-paph(i,k))
1321     .            *(pmfuq(i,k)+pmfdq(i,k)+pmful(i,k)
1322     .             +pdmfup(i,k)+pdmfdp(i,k))
1323         dq_con(i,k)=zdqdt
1324      ENDIF
1325  230 CONTINUE
1326C
1327      RETURN
1328      END
1329      SUBROUTINE flxdlfs(ptenh, pqenh, pgeoh, paph, ptu, pqu,
1330     .     ldcum, kcbot, kctop, pmfub, prfl, ptd, pqd,
1331     .     pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)
1332      IMPLICIT none
1333C
1334C----------------------------------------------------------------------
1335C THIS ROUTINE CALCULATES LEVEL OF FREE SINKING FOR
1336C CUMULUS DOWNDRAFTS AND SPECIFIES T,Q,U AND V VALUES
1337C
1338C TO PRODUCE LFS-VALUES FOR CUMULUS DOWNDRAFTS
1339C FOR MASSFLUX CUMULUS PARAMETERIZATION
1340C
1341C INPUT ARE ENVIRONMENTAL VALUES OF T,Q,U,V,P,PHI
1342C AND UPDRAFT VALUES T,Q,U AND V AND ALSO
1343C CLOUD BASE MASSFLUX AND CU-PRECIPITATION RATE.
1344C IT RETURNS T,Q,U AND V VALUES AND MASSFLUX AT LFS.
1345C
1346C CHECK FOR NEGATIVE BUOYANCY OF AIR OF EQUAL PARTS OF
1347C MOIST ENVIRONMENTAL AIR AND CLOUD AIR.
1348C----------------------------------------------------------------------
1349#include "dimensions.h"
1350#include "dimphy.h"
1351#include "YOMCST.h"
1352#include "YOETHF.h"
1353#include "YOECUMF.h"
1354C
1355      REAL ptenh(klon,klev)
1356      REAL pqenh(klon,klev)
1357      REAL pgeoh(klon,klev), paph(klon,klev+1)
1358      REAL ptu(klon,klev), pqu(klon,klev)
1359      REAL pmfub(klon)
1360      REAL prfl(klon)
1361C
1362      REAL ptd(klon,klev), pqd(klon,klev)
1363      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1364      REAL pdmfdp(klon,klev)
1365      INTEGER  kcbot(klon), kctop(klon), kdtop(klon)
1366      LOGICAL  ldcum(klon), lddraf(klon)
1367C
1368      REAL ztenwb(klon,klev), zqenwb(klon,klev), zcond(klon)
1369      REAL zttest, zqtest, zbuo, zmftop
1370      LOGICAL  llo2(klon)
1371      INTEGER i, k, is, icall
1372C----------------------------------------------------------------------
1373      DO i= 1, klon
1374         lddraf(i)=.FALSE.
1375         kdtop(i)=klev+1
1376      ENDDO
1377C
1378C----------------------------------------------------------------------
1379C DETERMINE LEVEL OF FREE SINKING BY
1380C DOING A SCAN FROM TOP TO BASE OF CUMULUS CLOUDS
1381C
1382C FOR EVERY POINT AND PROCEED AS FOLLOWS:
1383C     (1) DETEMINE WET BULB ENVIRONMENTAL T AND Q
1384C     (2) DO MIXING WITH CUMULUS CLOUD AIR
1385C     (3) CHECK FOR NEGATIVE BUOYANCY
1386C
1387C THE ASSUMPTION IS THAT AIR OF DOWNDRAFTS IS MIXTURE
1388C OF 50% CLOUD AIR + 50% ENVIRONMENTAL AIR AT WET BULB
1389C TEMPERATURE (I.E. WHICH BECAME SATURATED DUE TO
1390C EVAPORATION OF RAIN AND CLOUD WATER)
1391C----------------------------------------------------------------------
1392C
1393      DO 290 k = 3, klev-3
1394C
1395      is=0
1396      DO 212 i= 1, klon
1397         ztenwb(i,k)=ptenh(i,k)
1398         zqenwb(i,k)=pqenh(i,k)
1399         llo2(i) = ldcum(i).AND.prfl(i).GT.0.
1400     .             .AND..NOT.lddraf(i)
1401     .             .AND.(k.LT.kcbot(i).AND.k.GT.kctop(i))
1402         IF ( llo2(i) ) is = is + 1
1403  212 CONTINUE
1404      IF(is.EQ.0) GO TO 290
1405C
1406      icall=2
1407      CALL flxadjtq(paph(1,k), ztenwb(1,k), zqenwb(1,k), llo2, icall)
1408C
1409C----------------------------------------------------------------------
1410C DO MIXING OF CUMULUS AND ENVIRONMENTAL AIR
1411C AND CHECK FOR NEGATIVE BUOYANCY.
1412C THEN SET VALUES FOR DOWNDRAFT AT LFS.
1413C----------------------------------------------------------------------
1414      DO 222 i= 1, klon
1415      IF (llo2(i)) THEN
1416         zttest=0.5*(ptu(i,k)+ztenwb(i,k))
1417         zqtest=0.5*(pqu(i,k)+zqenwb(i,k))
1418         zbuo=zttest*(1.+RETV*zqtest)-
1419     .        ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1420         zcond(i)=pqenh(i,k)-zqenwb(i,k)
1421         zmftop=-CMFDEPS*pmfub(i)
1422         IF (zbuo.LT.0..AND.prfl(i).GT.10.*zmftop*zcond(i)) THEN
1423            kdtop(i)=k
1424            lddraf(i)=.TRUE.
1425            ptd(i,k)=zttest
1426            pqd(i,k)=zqtest
1427            pmfd(i,k)=zmftop
1428            pmfds(i,k)=pmfd(i,k)*(RCPD*ptd(i,k)+pgeoh(i,k))
1429            pmfdq(i,k)=pmfd(i,k)*pqd(i,k)
1430            pdmfdp(i,k-1)=-0.5*pmfd(i,k)*zcond(i)
1431            prfl(i)=prfl(i)+pdmfdp(i,k-1)
1432         ENDIF
1433      ENDIF
1434  222 CONTINUE
1435c
1436  290 CONTINUE
1437C
1438      RETURN
1439      END
1440      SUBROUTINE flxddraf(ptenh, pqenh, pgeoh, paph, prfl,
1441     .           ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp,
1442     .           lddraf, pen_d, pde_d)
1443      IMPLICIT none
1444C
1445C----------------------------------------------------------------------
1446C          THIS ROUTINE CALCULATES CUMULUS DOWNDRAFT DESCENT
1447C
1448C          TO PRODUCE THE VERTICAL PROFILES FOR CUMULUS DOWNDRAFTS
1449C          (I.E. T,Q,U AND V AND FLUXES)
1450C
1451C          INPUT IS T,Q,P,PHI,U,V AT HALF LEVELS.
1452C          IT RETURNS FLUXES OF S,Q AND EVAPORATION RATE
1453C          AND U,V AT LEVELS WHERE DOWNDRAFT OCCURS
1454C
1455C          CALCULATE MOIST DESCENT FOR ENTRAINING/DETRAINING PLUME BY
1456C          A) MOVING AIR DRY-ADIABATICALLY TO NEXT LEVEL BELOW AND
1457C          B) CORRECTING FOR EVAPORATION TO OBTAIN SATURATED STATE.
1458C
1459C----------------------------------------------------------------------
1460#include "dimensions.h"
1461#include "dimphy.h"
1462#include "YOMCST.h"
1463#include "YOETHF.h"
1464#include "YOECUMF.h"
1465C
1466      REAL ptenh(klon,klev), pqenh(klon,klev)
1467      REAL pgeoh(klon,klev), paph(klon,klev+1)
1468C
1469      REAL ptd(klon,klev), pqd(klon,klev)
1470      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1471      REAL pdmfdp(klon,klev)
1472      REAL prfl(klon)
1473      LOGICAL lddraf(klon)
1474C
1475      REAL pen_d(klon,klev), pde_d(klon,klev), zcond(klon)
1476      LOGICAL llo2(klon), llo1
1477      INTEGER i, k, is, icall, itopde
1478      REAL zentr, zseen, zqeen, zsdde, zqdde, zmfdsk, zmfdqk, zdmfdp
1479      REAL zbuo
1480C----------------------------------------------------------------------
1481C CALCULATE MOIST DESCENT FOR CUMULUS DOWNDRAFT BY
1482C       (A) CALCULATING ENTRAINMENT RATES, ASSUMING
1483C           LINEAR DECREASE OF MASSFLUX IN PBL
1484C       (B) DOING MOIST DESCENT - EVAPORATIVE COOLING
1485C           AND MOISTENING IS CALCULATED IN *flxadjtq*
1486C       (C) CHECKING FOR NEGATIVE BUOYANCY AND
1487C           SPECIFYING FINAL T,Q,U,V AND DOWNWARD FLUXES
1488C
1489      DO 180 k = 3, klev
1490c
1491      is = 0
1492      DO i = 1, klon
1493         llo2(i)=lddraf(i).AND.pmfd(i,k-1).LT.0.
1494         IF (llo2(i)) is = is + 1
1495      ENDDO
1496      IF (is.EQ.0) GOTO 180
1497c
1498      DO i = 1, klon
1499      IF (llo2(i)) THEN
1500         zentr = ENTRDD*pmfd(i,k-1)*RD*ptenh(i,k-1)/
1501     .           (RG*paph(i,k-1))*(paph(i,k)-paph(i,k-1))
1502         pen_d(i,k) = zentr
1503         pde_d(i,k) = zentr
1504      ENDIF
1505      ENDDO
1506c
1507      itopde = klev-2
1508      IF (k.GT.itopde) THEN
1509         DO i = 1, klon
1510         IF (llo2(i)) THEN
1511            pen_d(i,k)=0.
1512            pde_d(i,k)=pmfd(i,itopde)*
1513     .      (paph(i,k)-paph(i,k-1))/(paph(i,klev+1)-paph(i,itopde))
1514         ENDIF
1515         ENDDO
1516      ENDIF
1517C
1518      DO i = 1, klon
1519      IF (llo2(i)) THEN
1520         pmfd(i,k) = pmfd(i,k-1)+pen_d(i,k)-pde_d(i,k)
1521         zseen = (RCPD*ptenh(i,k-1)+pgeoh(i,k-1))*pen_d(i,k)
1522         zqeen = pqenh(i,k-1)*pen_d(i,k)
1523         zsdde = (RCPD*ptd(i,k-1)+pgeoh(i,k-1))*pde_d(i,k)
1524         zqdde = pqd(i,k-1)*pde_d(i,k)
1525         zmfdsk = pmfds(i,k-1)+zseen-zsdde
1526         zmfdqk = pmfdq(i,k-1)+zqeen-zqdde
1527         pqd(i,k) = zmfdqk*(1./MIN(-CMFCMIN,pmfd(i,k)))
1528         ptd(i,k) = (zmfdsk*(1./MIN(-CMFCMIN,pmfd(i,k)))-
1529     .               pgeoh(i,k))/RCPD
1530         ptd(i,k) = MIN(400.,ptd(i,k))
1531         ptd(i,k) = MAX(100.,ptd(i,k))
1532         zcond(i) = pqd(i,k)
1533      ENDIF
1534      ENDDO
1535C
1536      icall = 2
1537      CALL flxadjtq(paph(1,k), ptd(1,k), pqd(1,k), llo2, icall)
1538C
1539      DO i = 1, klon
1540      IF (llo2(i)) THEN
1541         zcond(i) = zcond(i)-pqd(i,k)
1542         zbuo = ptd(i,k)*(1.+RETV  *pqd(i,k))-
1543     .          ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1544         llo1 = zbuo.LT.0..AND.(prfl(i)-pmfd(i,k)*zcond(i).GT.0.)
1545         IF (.not.llo1) pmfd(i,k) = 0.0
1546         pmfds(i,k) = (RCPD*ptd(i,k)+pgeoh(i,k))*pmfd(i,k)
1547         pmfdq(i,k) = pqd(i,k)*pmfd(i,k)
1548         zdmfdp = -pmfd(i,k)*zcond(i)
1549         pdmfdp(i,k-1) = zdmfdp
1550         prfl(i) = prfl(i)+zdmfdp
1551      ENDIF
1552      ENDDO
1553c
1554  180 CONTINUE
1555      RETURN
1556      END
1557      SUBROUTINE flxadjtq(pp, pt, pq, ldflag, kcall)
1558      IMPLICIT none
1559c======================================================================
1560c Objet: ajustement entre T et Q
1561c======================================================================
1562C NOTE: INPUT PARAMETER kcall DEFINES CALCULATION AS
1563C        kcall=0    ENV. T AND QS IN*CUINI*
1564C        kcall=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
1565C        kcall=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1566C
1567#include "dimensions.h"
1568#include "dimphy.h"
1569#include "YOMCST.h"
1570C
1571      REAL pt(klon), pq(klon), pp(klon)
1572      LOGICAL ldflag(klon)
1573      INTEGER kcall
1574c
1575      REAL zcond(klon), zcond1
1576      REAL Z5alvcp, z5alscp, zalvdcp, zalsdcp
1577      REAL zdelta, zcvm5, zldcp, zqsat, zcor
1578      INTEGER is, i
1579#include "YOETHF.h"
1580#include "FCTTRE.h"
1581C
1582      z5alvcp = r5les*RLVTT/RCPD
1583      z5alscp = r5ies*RLSTT/RCPD
1584      zalvdcp = rlvtt/RCPD
1585      zalsdcp = rlstt/RCPD
1586C
1587
1588      DO i = 1, klon
1589         zcond(i) = 0.0
1590      ENDDO
1591
1592      DO 210 i =1, klon
1593      IF (ldflag(i)) THEN
1594         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1595         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1596         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1597         zqsat = R2ES*FOEEW(pt(i),zdelta) / pp(i)
1598         zqsat = MIN(0.5,zqsat)
1599         zcor = 1./(1.-RETV*zqsat)
1600         zqsat = zqsat*zcor
1601         zcond(i) = (pq(i)-zqsat)
1602     .     / (1. + FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor))
1603         IF (kcall.EQ.1) zcond(i) = MAX(zcond(i),0.)
1604         IF (kcall.EQ.2) zcond(i) = MIN(zcond(i),0.)
1605         pt(i) = pt(i) + zldcp*zcond(i)
1606         pq(i) = pq(i) - zcond(i)
1607      ENDIF
1608  210 CONTINUE
1609C
1610      is = 0
1611      DO i =1, klon
1612         IF (zcond(i).NE.0.) is = is + 1
1613      ENDDO
1614      IF (is.EQ.0) GOTO 230
1615C
1616      DO 220 i = 1, klon
1617      IF(ldflag(i).AND.zcond(i).NE.0.) THEN
1618         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1619         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1620         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1621         zqsat = R2ES* FOEEW(pt(i),zdelta) / pp(i)
1622         zqsat = MIN(0.5,zqsat)
1623         zcor = 1./(1.-RETV*zqsat)
1624         zqsat = zqsat*zcor
1625         zcond1 = (pq(i)-zqsat)
1626     .     / (1. + FOEDE(pt(i),zdelta,zcvm5,zqsat,zcor))
1627         pt(i) = pt(i) + zldcp*zcond1
1628         pq(i) = pq(i) - zcond1
1629      ENDIF
1630  220 CONTINUE
1631C
1632  230 CONTINUE
1633      RETURN
1634      END
1635      SUBROUTINE flxsetup
1636      IMPLICIT none
1637C
1638C     THIS ROUTINE DEFINES DISPOSABLE PARAMETERS FOR MASSFLUX SCHEME
1639C
1640#include "YOECUMF.h"
1641C
1642      ENTRPEN=1.0E-4  ! ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
1643      ENTRSCV=3.0E-4  ! ENTRAINMENT RATE FOR SHALLOW CONVECTION
1644      ENTRMID=1.0E-4  ! ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
1645      ENTRDD =2.0E-4  ! ENTRAINMENT RATE FOR DOWNDRAFTS
1646      CMFCTOP=0.33  ! RELATIVE CLOUD MASSFLUX AT LEVEL ABOVE NONBUO LEVEL
1647      CMFCMAX=1.0  ! MAXIMUM MASSFLUX VALUE ALLOWED FOR UPDRAFTS ETC
1648      CMFCMIN=1.E-10  ! MINIMUM MASSFLUX VALUE (FOR SAFETY)
1649      CMFDEPS=0.3  ! FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
1650      CPRCON =2.0E-4  ! CONVERSION FROM CLOUD WATER TO RAIN
1651      RHCDD=1.  ! RELATIVE SATURATION IN DOWNDRAFRS (NO LONGER USED)
1652c                 (FORMULATION IMPLIES SATURATION)
1653      LMFPEN = .TRUE.
1654      LMFSCV = .TRUE.
1655      LMFMID = .TRUE.
1656      LMFDD = .TRUE.
1657      LMFDUDV = .TRUE.
1658c
1659      RETURN
1660      END
Note: See TracBrowser for help on using the repository browser.