source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/conflx.F @ 1338

Last change on this file since 1338 was 766, checked in by Laurent Fairhead, 17 years ago

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