source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/phylmd/conflx.F @ 664

Last change on this file since 664 was 524, checked in by lmdzadmin, 21 years ago

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