source: LMDZ6/trunk/libf/phylmd/conflx.f90 @ 5322

Last change on this file since 5322 was 5289, checked in by abarral, 3 weeks ago

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