/[lmdze]/trunk/libf/phylmd/conflx.f
ViewVC logotype

Contents of /trunk/libf/phylmd/conflx.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21