source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/conflx.F @ 5080

Last change on this file since 5080 was 634, checked in by Laurent Fairhead, 20 years ago

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