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