source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/conflx.F90 @ 5080

Last change on this file since 5080 was 1999, checked in by Laurent Fairhead, 11 years ago

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