source: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conflx.F90 @ 3222

Last change on this file since 3222 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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