source: lmdz_wrf/trunk/WRFV3/lmdz/conflx.F90 @ 951

Last change on this file since 951 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 60.3 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE conflx (dtime,pres_h,pres_f,                                        &
5       &                   t, q, con_t, con_q, pqhfl, w,                             &
6       &                   d_t, d_q, rain, snow,                                     &
7       &                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,                   &
8       &                   kcbot, kctop, kdtop, pmflxr, pmflxs)
9!c
10      USE dimphy
11      IMPLICIT none
12!c======================================================================
13!c Auteur(s): Z.X. Li (LMD/CNRS) date: 19941014
14!c Objet: Schema flux de masse pour la convection
15!c        (schema de Tiedtke avec qqs modifications mineures)
16!c Dec.97: Prise en compte des modifications introduites par
17!c         Olivier Boucher et Alexandre Armengaud pour melange
18!c         et lessivage des traceurs passifs.
19!c======================================================================
20!cym#include "dimensions.h"
21!cym#include "dimphy.h"
22#include "YOMCST.h"
23#include "YOETHF.h"
24!c 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
34!c 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
50!c 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)
56!c
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)
64!cAA
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)
73!cAA
74
75!c
76      INTEGER i, k
77      REAL zdelta, zqsat
78!c
79#include "FCTTRE.h"
80!c
81!c 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
113!c
114!c calculer la nature du sol (pour l'instant, ocean partout)
115      DO i = 1, klon
116         land(i) = .FALSE.
117      ENDDO
118!c
119!c preparer les variables d'entree (attention: l'ordre des niveaux
120!c 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)
130!c
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
151!c
152!c appeler la routine principale
153!c
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)
159!C
160!cAA--------------------------------------------------------
161!cAA rem : De la meme facon que l'on effectue le reindicage
162!cAA       pour la temperature t et le champ q
163!cAA       on reindice les flux necessaires a la convection
164!cAA       des traceurs
165!cAA--------------------------------------------------------
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
172!c
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      DO k = 2, klev
180      DO i = 1, klon
181         pmfu(i,klev+2-k)= zmfu(i,k)
182         pmfd(i,klev+2-k)= zmfd(i,k)
183      ENDDO
184      ENDDO
185!c
186      DO k = 1, klev
187      DO i = 1, klon
188         pen_u(i,klev+1-k)=  zen_u(i,k)
189         pde_u(i,klev+1-k)=  zde_u(i,k)
190      ENDDO
191      ENDDO
192!c
193      DO k = 1, klev-1
194      DO i = 1, klon
195         pen_d(i,klev+1-k)= -zen_d(i,k+1)
196         pde_d(i,klev+1-k)= -zde_d(i,k+1)
197      ENDDO
198      ENDDO
199
200      DO k = 1, klev+1
201      DO i = 1, klon
202         pmflxr(i,klev+2-k)= zmflxr(i,k)
203         pmflxs(i,klev+2-k)= zmflxs(i,k)
204      ENDDO
205      ENDDO
206
207      RETURN
208      END
209!c--------------------------------------------------------------------
210      SUBROUTINE flxmain(pdtime, pten, pqen, pqsen, pqhfl, pap, paph,                &
211       &                   pgeo, ldland, ptte, pqte, pvervel,                        &
212       &                   prsfc, pssfc, kcbot, kctop, kdtop,                        &
213!c     *                   ldcum, ktype,                                             &
214       &                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,                   &
215       &                   dt_con, dq_con, pmflxr, pmflxs)
216      USE dimphy
217      IMPLICIT none
218!C     ------------------------------------------------------------------
219!cym#include "dimensions.h"
220!cym#include "dimphy.h"
221#include "YOMCST.h"
222#include "YOETHF.h"
223#include "YOECUMF.h"
224!C     ----------------------------------------------------------------
225      REAL pten(klon,klev), pqen(klon,klev), pqsen(klon,klev)
226      REAL ptte(klon,klev)
227      REAL pqte(klon,klev)
228      REAL pvervel(klon,klev)
229      REAL pgeo(klon,klev), pap(klon,klev), paph(klon,klev+1)
230      REAL pqhfl(klon)
231!c
232      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
233      REAL plude(klon,klev)
234      REAL pmfu(klon,klev)
235      REAL prsfc(klon), pssfc(klon)
236      INTEGER  kcbot(klon), kctop(klon), ktype(klon)
237      LOGICAL  ldland(klon), ldcum(klon)
238!c
239      REAL ztenh(klon,klev), zqenh(klon,klev), zqsenh(klon,klev)
240      REAL zgeoh(klon,klev)
241      REAL zmfub(klon), zmfub1(klon)
242      REAL zmfus(klon,klev), zmfuq(klon,klev), zmful(klon,klev)
243      REAL zdmfup(klon,klev), zdpmel(klon,klev)
244      REAL zentr(klon), zhcbase(klon)
245      REAL zdqpbl(klon), zdqcv(klon), zdhpbl(klon)
246      REAL zrfl(klon)
247      REAL pmflxr(klon,klev+1)
248      REAL pmflxs(klon,klev+1)
249      INTEGER  ilab(klon,klev), ictop0(klon)
250      LOGICAL  llo1
251      REAL dt_con(klon,klev), dq_con(klon,klev)
252      REAL zmfmax, zdh
253      REAL pdtime, zqumqe, zdqmin, zalvdcp, zhsat, zzz
254      REAL zhhat, zpbmpt, zgam, zeps, zfac
255      INTEGER i, k, ikb, itopm2, kcum
256!c
257      REAL pen_u(klon,klev), pde_u(klon,klev)
258      REAL pen_d(klon,klev), pde_d(klon,klev)
259!c
260      REAL ptd(klon,klev), pqd(klon,klev), pmfd(klon,klev)
261      REAL zmfds(klon,klev), zmfdq(klon,klev), zdmfdp(klon,klev)
262      INTEGER kdtop(klon)
263      LOGICAL lddraf(klon)
264!C---------------------------------------------------------------------
265      LOGICAL firstcal
266      SAVE firstcal
267      DATA firstcal / .TRUE. /
268!$OMP THREADPRIVATE(firstcal)
269!C---------------------------------------------------------------------
270      IF (firstcal) THEN
271         CALL flxsetup
272         firstcal = .FALSE.
273      ENDIF
274!C---------------------------------------------------------------------
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
284!c----------------------------------------------------------------------
285!c initialiser les variables et faire l'interpolation verticale
286!c----------------------------------------------------------------------
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)
292!c---------------------------------------------------------------------
293!c determiner les valeurs au niveau de base de la tour convective
294!c---------------------------------------------------------------------
295      CALL flxbase(ztenh, zqenh, zgeoh, paph,                                        &
296       &            ptu, pqu, plu, ldcum, kcbot, ilab)
297!c---------------------------------------------------------------------
298!c calculer la convergence totale de l'humidite et celle en provenance
299!c de la couche limite, plus precisement, la convergence integree entre
300!c le sol et la base de la convection. Cette derniere convergence est
301!c comparee avec l'evaporation obtenue dans la couche limite pour
302!c determiner le type de la convection
303!c---------------------------------------------------------------------
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
310!c
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
321!c
322      DO i = 1, klon
323         ktype(i) = 2
324         if (zdqcv(i).GT.MAX(0.,-1.5*pqhfl(i)*RG)) ktype(i) = 1
325!ccc         if (zdqcv(i).GT.MAX(0.,-1.1*pqhfl(i)*RG)) ktype(i) = 1
326      ENDDO
327!c
328!c---------------------------------------------------------------------
329!c determiner le flux de masse entrant a travers la base.
330!c on ignore, pour l'instant, l'effet du panache descendant
331!c---------------------------------------------------------------------
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
352!C-----------------------------------------------------------------------
353!C DETERMINE CLOUD ASCENT FOR ENTRAINING PLUME
354!C-----------------------------------------------------------------------
355!c (A) calculer d'abord la hauteur "theorique" de la tour convective sans
356!c     considerer l'entrainement ni le detrainement du panache, sachant
357!c     ces derniers peuvent abaisser la hauteur theorique.
358!c
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
364!c
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
377!c
378!c (B) calculer le panache ascendant
379!c
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
387!C
388!C verifier l'epaisseur de la convection et changer eventuellement
389!c le taux d'entrainement/detrainement
390!C
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
397!c
398      IF (lmfdd) THEN  ! si l'on considere le panache descendant
399!c
400!c calculer la precipitation issue du panache ascendant pour
401!c 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
410!c
411!c 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)
417!c
418!c 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)
424!c
425!c calculer de nouveau le flux de masse entrant a travers la base
426!c de la convection, sachant qu'il a ete modifie par le panache
427!c 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
471!c
472      ENDIF   ! fin de test sur lmfdd
473!c
474!c-----------------------------------------------------------------------
475!c calculer de nouveau le panache ascendant
476!c-----------------------------------------------------------------------
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)
483!c
484!c-----------------------------------------------------------------------
485!c determiner les flux convectifs en forme finale, ainsi que
486!c la quantite des precipitations
487!c-----------------------------------------------------------------------
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)
493!c
494!c----------------------------------------------------------------------
495!c calculer les tendances pour T et Q
496!c----------------------------------------------------------------------
497      CALL flxdtdq(pdtime, itopm2, paph, ldcum, pten,                                &
498       &     zmfus, zmfds, zmfuq, zmfdq, zmful, zdmfup, zdmfdp, zdpmel,              &
499       &     dt_con,dq_con)
500!c
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
510!C----------------------------------------------------------------------
511!C THIS ROUTINE INTERPOLATES LARGE-SCALE FIELDS OF T,Q ETC.
512!C TO HALF LEVELS (I.E. GRID FOR MASSFLUX SCHEME),
513!C AND INITIALIZES VALUES FOR UPDRAFTS
514!C----------------------------------------------------------------------
515!cym#include "dimensions.h"
516!cym#include "dimphy.h"
517#include "YOMCST.h"
518#include "YOETHF.h"
519!C
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
529!C
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
538!c                              p-a a l'environnement
539      REAL pdpmel(klon,klev) ! quantite de neige fondue
540!c
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
547!c
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
552!C
553      INTEGER  klab(klon,klev)
554      LOGICAL  llflag(klon)
555      INTEGER k, i, icall
556      REAL zzs
557!C----------------------------------------------------------------------
558!C SPECIFY LARGE SCALE PARAMETERS AT HALF LEVELS
559!C ADJUST TEMPERATURE FIELDS IF STATICLY UNSTABLE
560!C----------------------------------------------------------------------
561      DO 130 k = 2, klev
562!c
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
570!c
571      icall=0
572      CALL flxadjtq(paph(1,k),ptenh(1,k),pqsenh(1,k),llflag,icall)
573!c
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
579!c
580  130 CONTINUE
581!C
582      DO 140 i = 1, klon
583         ptenh(i,klev)=(RCPD*pten(i,klev)+pgeo(i,klev)-                              &
584       &                   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
590!c
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
598!C
599!C-----------------------------------------------------------------------
600!C INITIALIZE VALUES FOR UPDRAFTS AND DOWNDRAFTS
601!C-----------------------------------------------------------------------
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.
613!c
614         klab(i,k) = 0
615!c
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
622!c
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
629!C
630      RETURN
631      END
632      SUBROUTINE flxbase(ptenh, pqenh, pgeoh, paph,                                  &
633       &     ptu, pqu, plu, ldcum, kcbot, klab)
634      USE dimphy
635      IMPLICIT none
636!C----------------------------------------------------------------------
637!C THIS ROUTINE CALCULATES CLOUD BASE VALUES (T AND Q)
638!C
639!C INPUT ARE ENVIRONM. VALUES OF T,Q,P,PHI AT HALF LEVELS.
640!C IT RETURNS CLOUD BASE VALUES AND FLAGS AS FOLLOWS;
641!C   klab=1 FOR SUBCLOUD LEVELS
642!C   klab=2 FOR CONDENSATION LEVEL
643!C
644!C LIFT SURFACE AIR DRY-ADIABATICALLY TO CLOUD BASE
645!C (NON ENTRAINING PLUME,I.E.CONSTANT MASSFLUX)
646!C----------------------------------------------------------------------
647!cym#include "dimensions.h"
648!cym#include "dimphy.h"
649#include "YOMCST.h"
650#include "YOETHF.h"
651!C       ----------------------------------------------------------------
652      REAL ptenh(klon,klev), pqenh(klon,klev)
653      REAL pgeoh(klon,klev), paph(klon,klev+1)
654!C
655      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
656      INTEGER  klab(klon,klev), kcbot(klon)
657!C
658      LOGICAL llflag(klon), ldcum(klon)
659      INTEGER i, k, icall, is
660      REAL zbuo, zqold(klon)
661!C----------------------------------------------------------------------
662!C INITIALIZE VALUES AT LIFTING LEVEL
663!C----------------------------------------------------------------------
664      DO i = 1, klon
665         klab(i,klev)=1
666         kcbot(i)=klev-1
667         ldcum(i)=.FALSE.
668      ENDDO
669!C----------------------------------------------------------------------
670!C DO ASCENT IN SUBCLOUD LAYER,
671!C CHECK FOR EXISTENCE OF CONDENSATION LEVEL,
672!C ADJUST T,Q AND L ACCORDINGLY
673!C CHECK FOR BUOYANCY AND SET FLAGS
674!C----------------------------------------------------------------------
675      DO 290 k = klev-1, 2, -1
676!c
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
684!c
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
695!c
696      icall=1
697      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
698!c
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
709!c
710  290 CONTINUE
711!c
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
722!C----------------------------------------------------------------------
723!C THIS ROUTINE DOES THE CALCULATIONS FOR CLOUD ASCENTS
724!C FOR CUMULUS PARAMETERIZATION
725!C----------------------------------------------------------------------
726!cym#include "dimensions.h"
727!cym#include "dimphy.h"
728#include "YOMCST.h"
729#include "YOETHF.h"
730#include "YOECUMF.h"
731!C
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
739!C
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)
749!C
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
757!c
758      REAL zpbot(klon), zptop(klon), zrho(klon)
759      REAL zdprho, zentr, zpmid, zmftest, zmfmax
760      LOGICAL llo1, llo2
761!c
762      REAL zwmax(klon), zzzmb
763      INTEGER klwmin(klon) ! level of maximum vertical velocity
764!C----------------------------------------------------------------------
765      ztglace = RTT - 13.
766!c
767!c 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
780!C----------------------------------------------------------------------
781!C SET DEFAULT VALUES
782!C----------------------------------------------------------------------
783      DO i = 1, klon
784         IF (.NOT.ldcum(i)) ktype(i)=0
785      ENDDO
786!c
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
800!c
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
810!C
811!C Initialiser les valeurs au niveau d'ascendance
812!C
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
824!c
825      DO i = 1, klon
826         ldcum(i) = .FALSE.
827      ENDDO
828!C----------------------------------------------------------------------
829!C  DO ASCENT: SUBCLOUD LAYER (klab=1) ,CLOUDS (klab=2)
830!C  BY DOING FIRST DRY-ADIABATIC ASCENT AND THEN
831!C  BY ADJUSTING T,Q AND L ACCORDINGLY IN *flxadjtq*,
832!C  THEN CHECK FOR BUOYANCY AND SET FLAGS ACCORDINGLY
833!C----------------------------------------------------------------------
834      DO 480 k = klev-1,3,-1
835!c
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
858!c
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
867!c
868!c calculer le taux d'entrainement et de detrainement
869!c
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
877!c
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
905!c
906!C----------------------------------------------------------------------
907!c DO ADIABATIC ASCENT FOR ENTRAINING/DETRAINING PLUME
908!C----------------------------------------------------------------------
909!c
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))
918!c 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)
920!c 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       &               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
940!c
941!C----------------------------------------------------------------------
942!c DO CORRECTIONS FOR MOIST ASCENT BY ADJUSTING T,Q AND L
943!C----------------------------------------------------------------------
944!c
945      icall = 1
946      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
947!C
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
978!C
979  480 CONTINUE
980!C----------------------------------------------------------------------
981!C DETERMINE CONVECTIVE FLUXES ABOVE NON-BUOYANCY LEVEL
982!C    (NOTE: CLOUD VARIABLES LIKE T,Q AND L ARE NOT
983!C           AFFECTED BY DETRAINMENT AND ARE ALREADY KNOWN
984!C           FROM PREVIOUS CALCULATIONS ABOVE)
985!C----------------------------------------------------------------------
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
990!c
991      ldcum(1)=ldcum(1)
992!c
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
999!c
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
1015!C
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
1027!C----------------------------------------------------------------------
1028!C THIS ROUTINE DOES THE FINAL CALCULATION OF CONVECTIVE
1029!C FLUXES IN THE CLOUD LAYER AND IN THE SUBCLOUD LAYER
1030!C----------------------------------------------------------------------
1031!cym#include "dimensions.h"
1032!cym#include "dimphy.h"
1033#include "YOMCST.h"
1034#include "YOETHF.h"
1035#include "YOECUMF.h"
1036!C
1037      REAL cevapcu(klon,klev)
1038!C     -----------------------------------------------------------------
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)
1042!c
1043      REAL pap(klon,klev)
1044      REAL ztmsmlt, zdelta, zqsat
1045!C
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)
1052!cjq The variable maxpdmfdp(klon) has been introduced by Olivier Boucher
1053!cjq 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
1065!c
1066      LOGICAL lddraf(klon)
1067      INTEGER kdtop(klon)
1068!c
1069#include "FCTTRE.h"
1070!c
1071      DO 101 k=1,klev
1072        DO i=1,klon
1073      CEVAPCU(i,k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293)                               &
1074       & *SQRT(0.5*(paph(i,k)+paph(i,k+1))/paph(i,klev+1)) ) * 0.5/RG
1075        ENDDO
1076 101  CONTINUE
1077!c
1078!c SPECIFY CONSTANTS
1079!c
1080      zcons1 = RCPD/(RLMLT*RG*pdtime)
1081      zcons2 = 1./(RG*pdtime)
1082      zcucov = 0.05
1083      ztmelp2 = RTT + 2.
1084!c
1085!c DETERMINE FINAL CONVECTIVE FLUXES
1086!c
1087      itop=klev
1088      DO 110 i = 1, klon
1089         itop=MIN(itop,kctop(i))
1090         IF (.NOT.ldcum(i) .OR. kdtop(i).LT.kctop(i)) lddraf(i)=.FALSE.
1091         IF(.NOT.ldcum(i)) ktype(i)=0
1092  110 CONTINUE
1093!c
1094      ktopm2=itop-2
1095      DO 120 k=ktopm2,klev
1096      DO 115 i = 1, klon
1097      IF(ldcum(i).AND.k.GE.kctop(i)-1) THEN
1098         pmfus(i,k)=pmfus(i,k)-pmfu(i,k)*                                            &
1099       &                (RCPD*ptenh(i,k)+pgeoh(i,k))
1100         pmfuq(i,k)=pmfuq(i,k)-pmfu(i,k)*pqenh(i,k)
1101         zdp = 1.5E4
1102         IF ( ldland(i) ) zdp = 3.E4
1103!c
1104!c        l'eau liquide detrainee est precipitee quand certaines
1105!c        conditions sont reunies (sinon, elle est consideree
1106!c        evaporee dans l'environnement)
1107!c
1108         IF(paph(i,kcbot(i))-paph(i,kctop(i)).GE.zdp.AND.                            &
1109       &      pqen(i,k-1).GT.0.8*pqsen(i,k-1))                                       &
1110       &      pdmfup(i,k-1)=pdmfup(i,k-1)+plude(i,k-1)
1111!c
1112         IF(lddraf(i).AND.k.GE.kdtop(i)) THEN
1113            pmfds(i,k)=pmfds(i,k)-pmfd(i,k)*                                         &
1114       &                   (RCPD*ptenh(i,k)+pgeoh(i,k))
1115            pmfdq(i,k)=pmfdq(i,k)-pmfd(i,k)*pqenh(i,k)
1116         ELSE
1117            pmfd(i,k)=0.
1118            pmfds(i,k)=0.
1119            pmfdq(i,k)=0.
1120            pdmfdp(i,k-1)=0.
1121         END IF
1122      ELSE
1123         pmfu(i,k)=0.
1124         pmfus(i,k)=0.
1125         pmfuq(i,k)=0.
1126         pmful(i,k)=0.
1127         pdmfup(i,k-1)=0.
1128         plude(i,k-1)=0.
1129         pmfd(i,k)=0.
1130         pmfds(i,k)=0.
1131         pmfdq(i,k)=0.
1132         pdmfdp(i,k-1)=0.
1133      ENDIF
1134  115 CONTINUE
1135  120 CONTINUE
1136!c
1137      DO 130 k=ktopm2,klev
1138      DO 125 i = 1, klon
1139      IF(ldcum(i).AND.k.GT.kcbot(i)) THEN
1140         ikb=kcbot(i)
1141         zzp=((paph(i,klev+1)-paph(i,k))/                                            &
1142       &        (paph(i,klev+1)-paph(i,ikb)))
1143         IF (ktype(i).EQ.3) zzp = zzp**2
1144         pmfu(i,k)=pmfu(i,ikb)*zzp
1145         pmfus(i,k)=pmfus(i,ikb)*zzp
1146         pmfuq(i,k)=pmfuq(i,ikb)*zzp
1147         pmful(i,k)=pmful(i,ikb)*zzp
1148      ENDIF
1149  125 CONTINUE
1150  130 CONTINUE
1151!c
1152!c CALCULATE RAIN/SNOW FALL RATES
1153!c CALCULATE MELTING OF SNOW
1154!c CALCULATE EVAPORATION OF PRECIP
1155!c
1156      DO k = 1, klev+1
1157      DO i = 1, klon
1158         pmflxr(i,k) = 0.0
1159         pmflxs(i,k) = 0.0
1160      ENDDO
1161      ENDDO
1162      DO k = ktopm2, klev
1163      DO i = 1, klon
1164      IF (ldcum(i)) THEN
1165         IF (pmflxs(i,k).GT.0.0 .AND. pten(i,k).GT.ztmelp2) THEN
1166            zfac=zcons1*(paph(i,k+1)-paph(i,k))
1167            zsnmlt=MIN(pmflxs(i,k),zfac*(pten(i,k)-ztmelp2))
1168            pdpmel(i,k)=zsnmlt
1169            ztmsmlt=pten(i,k)-zsnmlt/zfac
1170            zdelta=MAX(0.,SIGN(1.,RTT-ztmsmlt))
1171            zqsat=R2ES*FOEEW(ztmsmlt, zdelta) / pap(i,k)
1172            zqsat=MIN(0.5,zqsat)
1173            zqsat=zqsat/(1.-RETV  *zqsat)
1174            pqsen(i,k) = zqsat
1175         ENDIF
1176         IF (pten(i,k).GT.RTT) THEN
1177         pmflxr(i,k+1)=pmflxr(i,k)+pdmfup(i,k)+pdmfdp(i,k)+pdpmel(i,k)
1178         pmflxs(i,k+1)=pmflxs(i,k)-pdpmel(i,k)
1179         ELSE
1180           pmflxs(i,k+1)=pmflxs(i,k)+pdmfup(i,k)+pdmfdp(i,k)
1181           pmflxr(i,k+1)=pmflxr(i,k)
1182         ENDIF
1183!c        si la precipitation est negative, on ajuste le plux du
1184!c        panache descendant pour eliminer la negativite
1185         IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN
1186            pdmfdp(i,k) = -pmflxr(i,k)-pmflxs(i,k)-pdmfup(i,k)
1187            pmflxr(i,k+1) = 0.0
1188            pmflxs(i,k+1) = 0.0
1189            pdpmel(i,k) = 0.0
1190         ENDIF
1191      ENDIF
1192      ENDDO
1193      ENDDO
1194!c
1195!cjq The new variable is initialized here.
1196!cjq It contains the humidity which is fed to the downdraft
1197!cjq by evaporation of precipitation in the column below the base
1198!cjq of convection.
1199!cjq
1200!cjq In the former version, this term has been subtracted from precip
1201!cjq as well as the evaporation.
1202!cjq     
1203      DO k = 1, klev
1204      DO i = 1, klon
1205         maxpdmfdp(i,k)=0.0
1206      ENDDO
1207      ENDDO
1208      DO k = 1, klev
1209       DO kp = k, klev
1210        DO i = 1, klon
1211         maxpdmfdp(i,k)=maxpdmfdp(i,k)+pdmfdp(i,kp)
1212        ENDDO
1213       ENDDO
1214      ENDDO
1215!cjq End of initialization
1216!c     
1217      DO k = ktopm2, klev
1218      DO i = 1, klon
1219      IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN
1220         zrfl = pmflxr(i,k) + pmflxs(i,k)
1221         IF (zrfl.GT.1.0E-20) THEN
1222            zrnew=(MAX(0.,SQRT(zrfl/zcucov)-                                         &
1223       &            CEVAPCU(i,k)*(paph(i,k+1)-paph(i,k))*                            &
1224       &            MAX(0.,pqsen(i,k)-pqen(i,k))))**2*zcucov
1225            zrmin=zrfl-zcucov*MAX(0.,0.8*pqsen(i,k)-pqen(i,k))                       &
1226       &            *zcons2*(paph(i,k+1)-paph(i,k))
1227            zrnew=MAX(zrnew,zrmin)
1228            zrfln=MAX(zrnew,0.)
1229            zdrfl=MIN(0.,zrfln-zrfl)
1230!cjq At least the amount of precipiation needed to feed the downdraft
1231!cjq with humidity below the base of convection has to be left and can't
1232!cjq be evaporated (surely the evaporation can't be positive):           
1233            zdrfl=MAX(zdrfl,                                                         &
1234       &            MIN(-pmflxr(i,k)-pmflxs(i,k)-maxpdmfdp(i,k),0.0))
1235!cjq End of insertion
1236!c           
1237            zdenom=1.0/MAX(1.0E-20,pmflxr(i,k)+pmflxs(i,k))
1238            IF (pten(i,k).GT.RTT) THEN
1239               zpdr = pdmfdp(i,k)
1240               zpds = 0.0
1241            ELSE
1242               zpdr = 0.0
1243               zpds = pdmfdp(i,k)
1244            ENDIF
1245            pmflxr(i,k+1) = pmflxr(i,k) + zpdr + pdpmel(i,k)                         &
1246       &                    + zdrfl*pmflxr(i,k)*zdenom
1247            pmflxs(i,k+1) = pmflxs(i,k) + zpds - pdpmel(i,k)                         &
1248       &                    + zdrfl*pmflxs(i,k)*zdenom
1249            pdmfup(i,k) = pdmfup(i,k) + zdrfl
1250         ELSE
1251            pmflxr(i,k+1) = 0.0
1252            pmflxs(i,k+1) = 0.0
1253            pdmfdp(i,k) = 0.0
1254            pdpmel(i,k) = 0.0
1255         ENDIF         
1256         if (pmflxr(i,k) + pmflxs(i,k).lt.-1.e-26)                                   &
1257       &    write(*,*) 'precip. < 1e-16 ',pmflxr(i,k) + pmflxs(i,k)
1258      ENDIF
1259      ENDDO
1260      ENDDO
1261!c
1262      DO 210 i = 1, klon
1263         prfl(i) = pmflxr(i,klev+1)
1264         psfl(i) = pmflxs(i,klev+1)
1265  210 CONTINUE
1266!c
1267      RETURN
1268      END
1269      SUBROUTINE flxdtdq(pdtime, ktopm2, paph, ldcum, pten                           &
1270       &  ,  pmfus, pmfds, pmfuq, pmfdq, pmful, pdmfup, pdmfdp                       &
1271       &  ,  pdpmel, dt_con, dq_con)
1272      USE dimphy
1273      IMPLICIT none
1274!c----------------------------------------------------------------------
1275!c calculer les tendances T et Q
1276!c----------------------------------------------------------------------
1277!cym#include "dimensions.h"
1278!cym#include "dimphy.h"
1279#include "YOMCST.h"
1280#include "YOETHF.h"
1281#include "YOECUMF.h"
1282!C     -----------------------------------------------------------------
1283      LOGICAL  llo1
1284!C
1285      REAL pten(klon,klev), paph(klon,klev+1)
1286      REAL pmfus(klon,klev), pmfuq(klon,klev), pmful(klon,klev)
1287      REAL pmfds(klon,klev), pmfdq(klon,klev)
1288      REAL pdmfup(klon,klev)
1289      REAL pdmfdp(klon,klev)
1290      REAL pdpmel(klon,klev)
1291      LOGICAL ldcum(klon)
1292      REAL dt_con(klon,klev), dq_con(klon,klev)
1293!c
1294      INTEGER ktopm2
1295      REAL pdtime
1296!c
1297      INTEGER i, k
1298      REAL zalv, zdtdt, zdqdt
1299!c
1300      DO 210 k=ktopm2,klev-1
1301      DO 220 i = 1, klon
1302      IF (ldcum(i)) THEN
1303         llo1 = (pten(i,k)-RTT).GT.0.
1304         zalv = RLSTT
1305         IF (llo1) zalv = RLVTT
1306         zdtdt=RG/(paph(i,k+1)-paph(i,k))/RCPD                                       &
1307       &        *(pmfus(i,k+1)-pmfus(i,k)                                            &
1308       &         +pmfds(i,k+1)-pmfds(i,k)                                            &
1309       &          -RLMLT*pdpmel(i,k)                                                 &
1310       &          -zalv*(pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))            &
1311       &         )
1312         dt_con(i,k)=zdtdt
1313         zdqdt=RG/(paph(i,k+1)-paph(i,k))                                            &
1314       &        *(pmfuq(i,k+1)-pmfuq(i,k)                                            &
1315       &         +pmfdq(i,k+1)-pmfdq(i,k)                                            &
1316       &          +pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1317         dq_con(i,k)=zdqdt
1318      ENDIF
1319  220 CONTINUE
1320  210 CONTINUE
1321!C
1322      k = klev
1323      DO 230 i = 1, klon
1324      IF (ldcum(i)) THEN
1325         llo1 = (pten(i,k)-RTT).GT.0.
1326         zalv = RLSTT
1327         IF (llo1) zalv = RLVTT
1328         zdtdt=-RG/(paph(i,k+1)-paph(i,k))/RCPD                                      &
1329       &         *(pmfus(i,k)+pmfds(i,k)+RLMLT*pdpmel(i,k)                           &
1330       &           -zalv*(pmful(i,k)+pdmfup(i,k)+pdmfdp(i,k)))
1331         dt_con(i,k)=zdtdt
1332         zdqdt=-RG/(paph(i,k+1)-paph(i,k))                                           &
1333       &            *(pmfuq(i,k)+pmfdq(i,k)+pmful(i,k)                               &
1334       &             +pdmfup(i,k)+pdmfdp(i,k))
1335         dq_con(i,k)=zdqdt
1336      ENDIF
1337  230 CONTINUE
1338!C
1339      RETURN
1340      END
1341      SUBROUTINE flxdlfs(ptenh, pqenh, pgeoh, paph, ptu, pqu,                        &
1342       &     ldcum, kcbot, kctop, pmfub, prfl, ptd, pqd,                             &
1343       &     pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)
1344      USE dimphy
1345      IMPLICIT none
1346!C
1347!C----------------------------------------------------------------------
1348!C THIS ROUTINE CALCULATES LEVEL OF FREE SINKING FOR
1349!C CUMULUS DOWNDRAFTS AND SPECIFIES T,Q,U AND V VALUES
1350!C
1351!C TO PRODUCE LFS-VALUES FOR CUMULUS DOWNDRAFTS
1352!C FOR MASSFLUX CUMULUS PARAMETERIZATION
1353!C
1354!C INPUT ARE ENVIRONMENTAL VALUES OF T,Q,U,V,P,PHI
1355!C AND UPDRAFT VALUES T,Q,U AND V AND ALSO
1356!C CLOUD BASE MASSFLUX AND CU-PRECIPITATION RATE.
1357!C IT RETURNS T,Q,U AND V VALUES AND MASSFLUX AT LFS.
1358!C
1359!C CHECK FOR NEGATIVE BUOYANCY OF AIR OF EQUAL PARTS OF
1360!C MOIST ENVIRONMENTAL AIR AND CLOUD AIR.
1361!C----------------------------------------------------------------------
1362!cym#include "dimensions.h"
1363!cym#include "dimphy.h"
1364#include "YOMCST.h"
1365#include "YOETHF.h"
1366#include "YOECUMF.h"
1367!C
1368      REAL ptenh(klon,klev)
1369      REAL pqenh(klon,klev)
1370      REAL pgeoh(klon,klev), paph(klon,klev+1)
1371      REAL ptu(klon,klev), pqu(klon,klev)
1372      REAL pmfub(klon)
1373      REAL prfl(klon)
1374!C
1375      REAL ptd(klon,klev), pqd(klon,klev)
1376      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1377      REAL pdmfdp(klon,klev)
1378      INTEGER  kcbot(klon), kctop(klon), kdtop(klon)
1379      LOGICAL  ldcum(klon), lddraf(klon)
1380!C
1381      REAL ztenwb(klon,klev), zqenwb(klon,klev), zcond(klon)
1382      REAL zttest, zqtest, zbuo, zmftop
1383      LOGICAL  llo2(klon)
1384      INTEGER i, k, is, icall
1385!C----------------------------------------------------------------------
1386      DO i= 1, klon
1387         lddraf(i)=.FALSE.
1388         kdtop(i)=klev+1
1389      ENDDO
1390!C
1391!C----------------------------------------------------------------------
1392!C DETERMINE LEVEL OF FREE SINKING BY
1393!C DOING A SCAN FROM TOP TO BASE OF CUMULUS CLOUDS
1394!C
1395!C FOR EVERY POINT AND PROCEED AS FOLLOWS:
1396!C     (1) DETEMINE WET BULB ENVIRONMENTAL T AND Q
1397!C     (2) DO MIXING WITH CUMULUS CLOUD AIR
1398!C     (3) CHECK FOR NEGATIVE BUOYANCY
1399!C
1400!C THE ASSUMPTION IS THAT AIR OF DOWNDRAFTS IS MIXTURE
1401!C OF 50% CLOUD AIR + 50% ENVIRONMENTAL AIR AT WET BULB
1402!C TEMPERATURE (I.E. WHICH BECAME SATURATED DUE TO
1403!C EVAPORATION OF RAIN AND CLOUD WATER)
1404!C----------------------------------------------------------------------
1405!C
1406      DO 290 k = 3, klev-3
1407!C
1408      is=0
1409      DO 212 i= 1, klon
1410         ztenwb(i,k)=ptenh(i,k)
1411         zqenwb(i,k)=pqenh(i,k)
1412         llo2(i) = ldcum(i).AND.prfl(i).GT.0.                                        &
1413       &             .AND..NOT.lddraf(i)                                             &
1414       &             .AND.(k.LT.kcbot(i).AND.k.GT.kctop(i))
1415         IF ( llo2(i) ) is = is + 1
1416  212 CONTINUE
1417      IF(is.EQ.0) GO TO 290
1418!C
1419      icall=2
1420      CALL flxadjtq(paph(1,k), ztenwb(1,k), zqenwb(1,k), llo2, icall)
1421!C
1422!C----------------------------------------------------------------------
1423!C DO MIXING OF CUMULUS AND ENVIRONMENTAL AIR
1424!C AND CHECK FOR NEGATIVE BUOYANCY.
1425!C THEN SET VALUES FOR DOWNDRAFT AT LFS.
1426!C----------------------------------------------------------------------
1427      DO 222 i= 1, klon
1428      IF (llo2(i)) THEN
1429         zttest=0.5*(ptu(i,k)+ztenwb(i,k))
1430         zqtest=0.5*(pqu(i,k)+zqenwb(i,k))
1431         zbuo=zttest*(1.+RETV*zqtest)-                                               &
1432       &        ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1433         zcond(i)=pqenh(i,k)-zqenwb(i,k)
1434         zmftop=-CMFDEPS*pmfub(i)
1435         IF (zbuo.LT.0..AND.prfl(i).GT.10.*zmftop*zcond(i)) THEN
1436            kdtop(i)=k
1437            lddraf(i)=.TRUE.
1438            ptd(i,k)=zttest
1439            pqd(i,k)=zqtest
1440            pmfd(i,k)=zmftop
1441            pmfds(i,k)=pmfd(i,k)*(RCPD*ptd(i,k)+pgeoh(i,k))
1442            pmfdq(i,k)=pmfd(i,k)*pqd(i,k)
1443            pdmfdp(i,k-1)=-0.5*pmfd(i,k)*zcond(i)
1444            prfl(i)=prfl(i)+pdmfdp(i,k-1)
1445         ENDIF
1446      ENDIF
1447  222 CONTINUE
1448!c
1449  290 CONTINUE
1450!C
1451      RETURN
1452      END
1453      SUBROUTINE flxddraf(ptenh, pqenh, pgeoh, paph, prfl,                           &
1454       &           ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp,                             &
1455       &           lddraf, pen_d, pde_d)
1456      USE dimphy
1457      IMPLICIT none
1458!C
1459!C----------------------------------------------------------------------
1460!C          THIS ROUTINE CALCULATES CUMULUS DOWNDRAFT DESCENT
1461!C
1462!C          TO PRODUCE THE VERTICAL PROFILES FOR CUMULUS DOWNDRAFTS
1463!C          (I.E. T,Q,U AND V AND FLUXES)
1464!C
1465!C          INPUT IS T,Q,P,PHI,U,V AT HALF LEVELS.
1466!C          IT RETURNS FLUXES OF S,Q AND EVAPORATION RATE
1467!C          AND U,V AT LEVELS WHERE DOWNDRAFT OCCURS
1468!C
1469!C          CALCULATE MOIST DESCENT FOR ENTRAINING/DETRAINING PLUME BY
1470!C          A) MOVING AIR DRY-ADIABATICALLY TO NEXT LEVEL BELOW AND
1471!C          B) CORRECTING FOR EVAPORATION TO OBTAIN SATURATED STATE.
1472!C
1473!C----------------------------------------------------------------------
1474!cym#include "dimensions.h"
1475!cym#include "dimphy.h"
1476#include "YOMCST.h"
1477#include "YOETHF.h"
1478#include "YOECUMF.h"
1479!C
1480      REAL ptenh(klon,klev), pqenh(klon,klev)
1481      REAL pgeoh(klon,klev), paph(klon,klev+1)
1482!C
1483      REAL ptd(klon,klev), pqd(klon,klev)
1484      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1485      REAL pdmfdp(klon,klev)
1486      REAL prfl(klon)
1487      LOGICAL lddraf(klon)
1488!C
1489      REAL pen_d(klon,klev), pde_d(klon,klev), zcond(klon)
1490      LOGICAL llo2(klon), llo1
1491      INTEGER i, k, is, icall, itopde
1492      REAL zentr, zseen, zqeen, zsdde, zqdde, zmfdsk, zmfdqk, zdmfdp
1493      REAL zbuo
1494!C----------------------------------------------------------------------
1495!C CALCULATE MOIST DESCENT FOR CUMULUS DOWNDRAFT BY
1496!C       (A) CALCULATING ENTRAINMENT RATES, ASSUMING
1497!C           LINEAR DECREASE OF MASSFLUX IN PBL
1498!C       (B) DOING MOIST DESCENT - EVAPORATIVE COOLING
1499!C           AND MOISTENING IS CALCULATED IN *flxadjtq*
1500!C       (C) CHECKING FOR NEGATIVE BUOYANCY AND
1501!C           SPECIFYING FINAL T,Q,U,V AND DOWNWARD FLUXES
1502!C
1503      DO 180 k = 3, klev
1504!c
1505      is = 0
1506      DO i = 1, klon
1507         llo2(i)=lddraf(i).AND.pmfd(i,k-1).LT.0.
1508         IF (llo2(i)) is = is + 1
1509      ENDDO
1510      IF (is.EQ.0) GOTO 180
1511!c
1512      DO i = 1, klon
1513      IF (llo2(i)) THEN
1514         zentr = ENTRDD*pmfd(i,k-1)*RD*ptenh(i,k-1)/                                 &
1515       &           (RG*paph(i,k-1))*(paph(i,k)-paph(i,k-1))
1516         pen_d(i,k) = zentr
1517         pde_d(i,k) = zentr
1518      ENDIF
1519      ENDDO
1520!c
1521      itopde = klev-2
1522      IF (k.GT.itopde) THEN
1523         DO i = 1, klon
1524         IF (llo2(i)) THEN
1525            pen_d(i,k)=0.
1526            pde_d(i,k)=pmfd(i,itopde)*                                               &
1527       &      (paph(i,k)-paph(i,k-1))/(paph(i,klev+1)-paph(i,itopde))
1528         ENDIF
1529         ENDDO
1530      ENDIF
1531!C
1532      DO i = 1, klon
1533      IF (llo2(i)) THEN
1534         pmfd(i,k) = pmfd(i,k-1)+pen_d(i,k)-pde_d(i,k)
1535         zseen = (RCPD*ptenh(i,k-1)+pgeoh(i,k-1))*pen_d(i,k)
1536         zqeen = pqenh(i,k-1)*pen_d(i,k)
1537         zsdde = (RCPD*ptd(i,k-1)+pgeoh(i,k-1))*pde_d(i,k)
1538         zqdde = pqd(i,k-1)*pde_d(i,k)
1539         zmfdsk = pmfds(i,k-1)+zseen-zsdde
1540         zmfdqk = pmfdq(i,k-1)+zqeen-zqdde
1541         pqd(i,k) = zmfdqk*(1./MIN(-CMFCMIN,pmfd(i,k)))
1542         ptd(i,k) = (zmfdsk*(1./MIN(-CMFCMIN,pmfd(i,k)))-                            &
1543       &               pgeoh(i,k))/RCPD
1544         ptd(i,k) = MIN(400.,ptd(i,k))
1545         ptd(i,k) = MAX(100.,ptd(i,k))
1546         zcond(i) = pqd(i,k)
1547      ENDIF
1548      ENDDO
1549!C
1550      icall = 2
1551      CALL flxadjtq(paph(1,k), ptd(1,k), pqd(1,k), llo2, icall)
1552!C
1553      DO i = 1, klon
1554      IF (llo2(i)) THEN
1555         zcond(i) = zcond(i)-pqd(i,k)
1556         zbuo = ptd(i,k)*(1.+RETV  *pqd(i,k))-                                       &
1557       &          ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1558         llo1 = zbuo.LT.0..AND.(prfl(i)-pmfd(i,k)*zcond(i).GT.0.)
1559         IF (.not.llo1) pmfd(i,k) = 0.0
1560         pmfds(i,k) = (RCPD*ptd(i,k)+pgeoh(i,k))*pmfd(i,k)
1561         pmfdq(i,k) = pqd(i,k)*pmfd(i,k)
1562         zdmfdp = -pmfd(i,k)*zcond(i)
1563         pdmfdp(i,k-1) = zdmfdp
1564         prfl(i) = prfl(i)+zdmfdp
1565      ENDIF
1566      ENDDO
1567!c
1568  180 CONTINUE
1569      RETURN
1570      END
1571      SUBROUTINE flxadjtq(pp, pt, pq, ldflag, kcall)
1572      USE dimphy
1573      IMPLICIT none
1574!c======================================================================
1575!c Objet: ajustement entre T et Q
1576!c======================================================================
1577!C NOTE: INPUT PARAMETER kcall DEFINES CALCULATION AS
1578!C        kcall=0    ENV. T AND QS IN*CUINI*
1579!C        kcall=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
1580!C        kcall=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1581!C
1582!cym#include "dimensions.h"
1583!cym#include "dimphy.h"
1584#include "YOMCST.h"
1585!C
1586      REAL pt(klon), pq(klon), pp(klon)
1587      LOGICAL ldflag(klon)
1588      INTEGER kcall
1589!c
1590      REAL zcond(klon), zcond1
1591      REAL Z5alvcp, z5alscp, zalvdcp, zalsdcp
1592      REAL zdelta, zcvm5, zldcp, zqsat, zcor
1593      INTEGER is, i
1594#include "YOETHF.h"
1595#include "FCTTRE.h"
1596!C
1597      z5alvcp = r5les*RLVTT/RCPD
1598      z5alscp = r5ies*RLSTT/RCPD
1599      zalvdcp = rlvtt/RCPD
1600      zalsdcp = rlstt/RCPD
1601!C
1602
1603      DO i = 1, klon
1604         zcond(i) = 0.0
1605      ENDDO
1606
1607      DO 210 i =1, klon
1608      IF (ldflag(i)) THEN
1609         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1610         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1611         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1612         zqsat = R2ES*FOEEW(pt(i),zdelta) / pp(i)
1613         zqsat = MIN(0.5,zqsat)
1614         zcor = 1./(1.-RETV*zqsat)
1615         zqsat = zqsat*zcor
1616         zcond(i) = (pq(i)-zqsat)                                                    &
1617       &     / (1. + FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor))
1618         IF (kcall.EQ.1) zcond(i) = MAX(zcond(i),0.)
1619         IF (kcall.EQ.2) zcond(i) = MIN(zcond(i),0.)
1620         pt(i) = pt(i) + zldcp*zcond(i)
1621         pq(i) = pq(i) - zcond(i)
1622      ENDIF
1623  210 CONTINUE
1624!C
1625      is = 0
1626      DO i =1, klon
1627         IF (zcond(i).NE.0.) is = is + 1
1628      ENDDO
1629      IF (is.EQ.0) GOTO 230
1630!C
1631      DO 220 i = 1, klon
1632      IF(ldflag(i).AND.zcond(i).NE.0.) THEN
1633         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1634         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1635         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1636         zqsat = R2ES* FOEEW(pt(i),zdelta) / pp(i)
1637         zqsat = MIN(0.5,zqsat)
1638         zcor = 1./(1.-RETV*zqsat)
1639         zqsat = zqsat*zcor
1640         zcond1 = (pq(i)-zqsat)                                                      &
1641       &     / (1. + FOEDE(pt(i),zdelta,zcvm5,zqsat,zcor))
1642         pt(i) = pt(i) + zldcp*zcond1
1643         pq(i) = pq(i) - zcond1
1644      ENDIF
1645  220 CONTINUE
1646!C
1647  230 CONTINUE
1648      RETURN
1649      END
1650      SUBROUTINE flxsetup
1651      IMPLICIT none
1652!C
1653!C     THIS ROUTINE DEFINES DISPOSABLE PARAMETERS FOR MASSFLUX SCHEME
1654!C
1655#include "YOECUMF.h"
1656!C
1657      ENTRPEN=1.0E-4  ! ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
1658      ENTRSCV=3.0E-4  ! ENTRAINMENT RATE FOR SHALLOW CONVECTION
1659      ENTRMID=1.0E-4  ! ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
1660      ENTRDD =2.0E-4  ! ENTRAINMENT RATE FOR DOWNDRAFTS
1661      CMFCTOP=0.33  ! RELATIVE CLOUD MASSFLUX AT LEVEL ABOVE NONBUO LEVEL
1662      CMFCMAX=1.0  ! MAXIMUM MASSFLUX VALUE ALLOWED FOR UPDRAFTS ETC
1663      CMFCMIN=1.E-10  ! MINIMUM MASSFLUX VALUE (FOR SAFETY)
1664      CMFDEPS=0.3  ! FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
1665      CPRCON =2.0E-4  ! CONVERSION FROM CLOUD WATER TO RAIN
1666      RHCDD=1.  ! RELATIVE SATURATION IN DOWNDRAFRS (NO LONGER USED)
1667!c                 (FORMULATION IMPLIES SATURATION)
1668      LMFPEN = .TRUE.
1669      LMFSCV = .TRUE.
1670      LMFMID = .TRUE.
1671      LMFDD = .TRUE.
1672      LMFDUDV = .TRUE.
1673!c
1674      RETURN
1675      END
Note: See TracBrowser for help on using the repository browser.