GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/conflx.F90 Lines: 0 738 0.0 %
Date: 2023-06-30 12:56:34 Branches: 0 506 0.0 %

Line Branch Exec Source
1
2
! $Header$
3
4
SUBROUTINE 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
201
END SUBROUTINE conflx
202
! --------------------------------------------------------------------
203
SUBROUTINE 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
474
1000 CONTINUE
475
  RETURN
476
END SUBROUTINE flxmain
477
SUBROUTINE 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
599
END SUBROUTINE flxini
600
SUBROUTINE 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
676
290 END DO
677
678
  RETURN
679
END SUBROUTINE flxbase
680
SUBROUTINE 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
939
480 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
976
800 CONTINUE
977
  RETURN
978
END SUBROUTINE flxasc
979
SUBROUTINE 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
1224
END SUBROUTINE flxflux
1225
SUBROUTINE 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
1286
END SUBROUTINE flxdtdq
1287
SUBROUTINE 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
1390
290 END DO
1391
1392
  RETURN
1393
END SUBROUTINE flxdlfs
1394
SUBROUTINE 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
1505
180 END DO
1506
  RETURN
1507
END SUBROUTINE flxddraf
1508
SUBROUTINE 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
1580
230 CONTINUE
1581
  RETURN
1582
END SUBROUTINE flxadjtq
1583
SUBROUTINE 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
1608
END SUBROUTINE flxsetup