source: LMDZ6/branches/Amaury_dev/libf/phylmd/conflx.F90 @ 5157

Last change on this file since 5157 was 5153, checked in by abarral, 7 weeks ago

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