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

Last change on this file since 41 was 41, checked in by lmdz, 24 years ago

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