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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 51 - (show annotations)
Tue Sep 20 09:14:34 2011 UTC (12 years, 7 months ago) by guez
File size: 56613 byte(s)
Split "getincom.f90" into "getincom.f90" and "getincom2.f90". Split
"nuage.f" into "nuage.f90", "diagcld1.f90" and "diagcld2.f90". Created
module "chem" from included file "chem.h". Moved "YOEGWD.f90" to
directory "Orography".

In "physiq", for evaporation of water, "zlsdcp" was equal to
"zlvdc". Removed useless variables.

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

  ViewVC Help
Powered by ViewVC 1.1.21