source: LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_cv.F90 @ 5157

Last change on this file since 5157 was 5142, checked in by abarral, 8 weeks ago

Put cvparam.h, fcg_gcssold.h, planete.h, tsoilnudge.h, YOECUMF.h into modules
  • Property svn:keywords set to Id
File size: 53.1 KB
Line 
1! $Id: lmdz_cv.F90 5142 2024-07-29 13:07:34Z abarral $
2
3MODULE lmdz_cv
4  !------------------------------------------------------------
5  ! Parameters for convectL:
6  ! (includes - microphysical parameters,
7  !                     - parameters that control the rate of approach
8  !               to quasi-equilibrium)
9  !                     - noff & minorig (previously in input of convect1)
10  !------------------------------------------------------------
11
12  IMPLICIT NONE; PRIVATE
13  PUBLIC elcrit, tlcrit, entp, sigs, sigd, omtrain, omtsnow, coeffr, coeffs &
14          , dtmax, cu, betad, alpha, damp, delta, noff, minorig, nl, nlp, nlm, &
15          cv_param, cv_prelim, cv_feed, cv_undilute1, cv_trigger, cv_compress, &
16          cv_undilute2, cv_closure, cv_mixing, cv_unsat, cv_yield, cv_uncompress
17
18  INTEGER noff, minorig, nl, nlp, nlm
19  REAL elcrit, tlcrit
20  REAL entp
21  REAL sigs, sigd
22  REAL omtrain, omtsnow, coeffr, coeffs
23  REAL dtmax
24  REAL cu
25  REAL betad
26  REAL alpha, damp
27  REAL delta
28
29  !$OMP THREADPRIVATE(elcrit, tlcrit, entp, sigs, sigd, omtrain, omtsnow, coeffr, coeffs &
30  !$OMP      , dtmax, cu, betad, alpha, damp, delta, noff, minorig, nl, nlp, nlm)
31
32CONTAINS
33
34SUBROUTINE cv_param(nd)
35  IMPLICIT NONE
36
37  ! ------------------------------------------------------------
38  ! Set parameters for convectL
39  ! (includes microphysical parameters and parameters that
40  ! control the rate of approach to quasi-equilibrium)
41  ! ------------------------------------------------------------
42
43  ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
44  ! ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
45  ! ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
46  ! ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
47  ! ***               BETWEEN 0 C AND TLCRIT)                        ***
48  ! ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
49  ! ***                       FORMULATION                            ***
50  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
51  ! ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
52  ! ***                        OF CLOUD                              ***
53  ! ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
54  ! ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
55  ! ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
56  ! ***                          OF RAIN                             ***
57  ! ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
58  ! ***                          OF SNOW                             ***
59  ! ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
60  ! ***                         TRANSPORT                            ***
61  ! ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
62  ! ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
63  ! ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
64  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
65  ! ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
66  ! ***                   (DAMP MUST BE LESS THAN 1)                 ***
67
68    INTEGER nd
69  CHARACTER (LEN = 20) :: modname = 'cv_routines'
70  CHARACTER (LEN = 80) :: abort_message
71
72  ! noff: integer limit for convection (nd-noff)
73  ! minorig: First level of convection
74
75  noff = 2
76  minorig = 2
77
78  nl = nd - noff
79  nlp = nl + 1
80  nlm = nl - 1
81
82  elcrit = 0.0011
83  tlcrit = -55.0
84  entp = 1.5
85  sigs = 0.12
86  sigd = 0.05
87  omtrain = 50.0
88  omtsnow = 5.5
89  coeffr = 1.0
90  coeffs = 0.8
91  dtmax = 0.9
92
93  cu = 0.70
94
95  betad = 10.0
96
97  damp = 0.1
98  alpha = 0.2
99
100  delta = 0.01 ! cld
101
102END SUBROUTINE cv_param
103
104SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
105  USE lmdz_cvthermo
106
107  IMPLICIT NONE
108
109  ! =====================================================================
110  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
111  ! =====================================================================
112
113  ! inputs:
114  INTEGER len, nd, ndp1
115  REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
116
117  ! outputs:
118  REAL lv(len, nd), cpn(len, nd), tv(len, nd)
119  REAL gz(len, nd), h(len, nd), hm(len, nd)
120
121  ! local variables:
122  INTEGER k, i
123  REAL cpx(len, nd)
124
125 
126  DO k = 1, nlp
127    DO i = 1, len
128      lv(i, k) = lv0 - clmcpv * (t(i, k) - t0)
129      cpn(i, k) = cpd * (1.0 - q(i, k)) + cpv * q(i, k)
130      cpx(i, k) = cpd * (1.0 - q(i, k)) + cl * q(i, k)
131      tv(i, k) = t(i, k) * (1.0 + q(i, k) * epsim1)
132    END DO
133  END DO
134
135  ! gz = phi at the full levels (same as p).
136
137  DO i = 1, len
138    gz(i, 1) = 0.0
139  END DO
140  DO k = 2, nlp
141    DO i = 1, len
142      gz(i, k) = gz(i, k - 1) + hrd * (tv(i, k - 1) + tv(i, k)) * (p(i, k - 1) - p(i, k)) / ph(i, &
143              k)
144    END DO
145  END DO
146
147  ! h  = phi + cpT (dry static energy).
148  ! hm = phi + cp(T-Tbase)+Lq
149
150  DO k = 1, nlp
151    DO i = 1, len
152      h(i, k) = gz(i, k) + cpn(i, k) * t(i, k)
153      hm(i, k) = gz(i, k) + cpx(i, k) * (t(i, k) - t(i, 1)) + lv(i, k) * q(i, k)
154    END DO
155  END DO
156
157END SUBROUTINE cv_prelim
158
159SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, &
160        qnk, gznk, plcl)
161  IMPLICIT NONE
162
163  ! ================================================================
164  ! Purpose: CONVECTIVE FEED
165  ! ================================================================
166
167 
168  ! inputs:
169  INTEGER len, nd
170  REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
171  REAL hm(len, nd), gz(len, nd)
172
173  ! outputs:
174  INTEGER iflag(len), nk(len), icb(len), icbmax
175  REAL tnk(len), qnk(len), gznk(len), plcl(len)
176
177  ! local variables:
178  INTEGER i, k
179  INTEGER ihmin(len)
180  REAL work(len)
181  REAL pnk(len), qsnk(len), rh(len), chi(len)
182
183  ! -------------------------------------------------------------------
184  ! --- Find level of minimum moist static energy
185  ! --- If level of minimum moist static energy coincides with
186  ! --- or is lower than minimum allowable parcel origin level,
187  ! --- set iflag to 6.
188  ! -------------------------------------------------------------------
189
190  DO i = 1, len
191    work(i) = 1.0E12
192    ihmin(i) = nl
193  END DO
194  DO k = 2, nlp
195    DO i = 1, len
196      IF ((hm(i, k)<work(i)) .AND. (hm(i, k)<hm(i, k - 1))) THEN
197        work(i) = hm(i, k)
198        ihmin(i) = k
199      END IF
200    END DO
201  END DO
202  DO i = 1, len
203    ihmin(i) = min(ihmin(i), nlm)
204    IF (ihmin(i)<=minorig) THEN
205      iflag(i) = 6
206    END IF
207  END DO
208
209  ! -------------------------------------------------------------------
210  ! --- Find that model level below the level of minimum moist static
211  ! --- energy that has the maximum value of moist static energy
212  ! -------------------------------------------------------------------
213
214  DO i = 1, len
215    work(i) = hm(i, minorig)
216    nk(i) = minorig
217  END DO
218  DO k = minorig + 1, nl
219    DO i = 1, len
220      IF ((hm(i, k)>work(i)) .AND. (k<=ihmin(i))) THEN
221        work(i) = hm(i, k)
222        nk(i) = k
223      END IF
224    END DO
225  END DO
226  ! -------------------------------------------------------------------
227  ! --- Check whether parcel level temperature and specific humidity
228  ! --- are reasonable
229  ! -------------------------------------------------------------------
230  DO i = 1, len
231    IF (((t(i, nk(i))<250.0) .OR. (q(i, nk(i))<=0.0) .OR. (p(i, ihmin(i))< &
232            400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
233  END DO
234  ! -------------------------------------------------------------------
235  ! --- Calculate lifted condensation level of air at parcel origin level
236  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
237  ! -------------------------------------------------------------------
238  DO i = 1, len
239    tnk(i) = t(i, nk(i))
240    qnk(i) = q(i, nk(i))
241    gznk(i) = gz(i, nk(i))
242    pnk(i) = p(i, nk(i))
243    qsnk(i) = qs(i, nk(i))
244
245    rh(i) = qnk(i) / qsnk(i)
246    rh(i) = min(1.0, rh(i))
247    chi(i) = tnk(i) / (1669.0 - 122.0 * rh(i) - tnk(i))
248    plcl(i) = pnk(i) * (rh(i)**chi(i))
249    IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
250            ) = 8
251  END DO
252  ! -------------------------------------------------------------------
253  ! --- Calculate first level above lcl (=icb)
254  ! -------------------------------------------------------------------
255  DO i = 1, len
256    icb(i) = nlm
257  END DO
258
259  DO k = minorig, nl
260    DO i = 1, len
261      IF ((k>=(nk(i) + 1)) .AND. (p(i, k)<plcl(i))) icb(i) = min(icb(i), k)
262    END DO
263  END DO
264
265  DO i = 1, len
266    IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
267  END DO
268
269  ! Compute icbmax.
270
271  icbmax = 2
272  DO i = 1, len
273    icbmax = max(icbmax, icb(i))
274  END DO
275
276END SUBROUTINE cv_feed
277
278SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, &
279        clw)
280  USE lmdz_cvthermo
281
282  IMPLICIT NONE
283
284 
285  ! inputs:
286  INTEGER len, nd
287  INTEGER nk(len), icb(len), icbmax
288  REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
289  REAL p(len, nd)
290
291  ! outputs:
292  REAL tp(len, nd), tvp(len, nd), clw(len, nd)
293
294  ! local variables:
295  INTEGER i, k
296  REAL tg, qg, alv, s, ahg, tc, denom, es, rg
297  REAL ah0(len), cpp(len)
298  REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
299
300  ! -------------------------------------------------------------------
301  ! --- Calculates the lifted parcel virtual temperature at nk,
302  ! --- the actual temperature, and the adiabatic
303  ! --- liquid water content. The procedure is to solve the equation.
304  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
305  ! -------------------------------------------------------------------
306
307  DO i = 1, len
308    tnk(i) = t(i, nk(i))
309    qnk(i) = q(i, nk(i))
310    gznk(i) = gz(i, nk(i))
311    ticb(i) = t(i, icb(i))
312    gzicb(i) = gz(i, icb(i))
313  END DO
314
315  ! ***  Calculate certain parcel quantities, including static energy   ***
316
317  DO i = 1, len
318    ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
319            273.15)) + gznk(i)
320    cpp(i) = cpd * (1. - qnk(i)) + qnk(i) * cpv
321  END DO
322
323  ! ***   Calculate lifted parcel quantities below cloud base   ***
324
325  DO k = minorig, icbmax - 1
326    DO i = 1, len
327      tp(i, k) = tnk(i) - (gz(i, k) - gznk(i)) / cpp(i)
328      tvp(i, k) = tp(i, k) * (1. + qnk(i) * epsi)
329    END DO
330  END DO
331
332  ! ***  Find lifted parcel quantities above cloud base    ***
333
334  DO i = 1, len
335    tg = ticb(i)
336    qg = qs(i, icb(i))
337    alv = lv0 - clmcpv * (ticb(i) - t0)
338
339    ! First iteration.
340
341    s = cpd + alv * alv * qg / (rrv * ticb(i) * ticb(i))
342    s = 1. / s
343    ahg = cpd * tg + (cl - cpd) * qnk(i) * ticb(i) + alv * qg + gzicb(i)
344    tg = tg + s * (ah0(i) - ahg)
345    tg = max(tg, 35.0)
346    tc = tg - t0
347    denom = 243.5 + tc
348    IF (tc>=0.0) THEN
349      es = 6.112 * exp(17.67 * tc / denom)
350    ELSE
351      es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
352    END IF
353    qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
354
355    ! Second iteration.
356
357    s = cpd + alv * alv * qg / (rrv * ticb(i) * ticb(i))
358    s = 1. / s
359    ahg = cpd * tg + (cl - cpd) * qnk(i) * ticb(i) + alv * qg + gzicb(i)
360    tg = tg + s * (ah0(i) - ahg)
361    tg = max(tg, 35.0)
362    tc = tg - t0
363    denom = 243.5 + tc
364    IF (tc>=0.0) THEN
365      es = 6.112 * exp(17.67 * tc / denom)
366    ELSE
367      es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
368    END IF
369    qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
370
371    alv = lv0 - clmcpv * (ticb(i) - 273.15)
372    tp(i, icb(i)) = (ah0(i) - (cl - cpd) * qnk(i) * ticb(i) - gz(i, icb(i)) - alv * qg) / cpd
373    clw(i, icb(i)) = qnk(i) - qg
374    clw(i, icb(i)) = max(0.0, clw(i, icb(i)))
375    rg = qg / (1. - qnk(i))
376    tvp(i, icb(i)) = tp(i, icb(i)) * (1. + rg * epsi)
377  END DO
378
379  DO k = minorig, icbmax
380    DO i = 1, len
381      tvp(i, k) = tvp(i, k) - tp(i, k) * qnk(i)
382    END DO
383  END DO
384
385END SUBROUTINE cv_undilute1
386
387SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
388  IMPLICIT NONE
389
390  ! -------------------------------------------------------------------
391  ! --- Test for instability.
392  ! --- If there was no convection at last time step and parcel
393  ! --- is stable at icb, then set iflag to 4.
394  ! -------------------------------------------------------------------
395
396 
397  ! inputs:
398  INTEGER len, nd, icb(len)
399  REAL cbmf(len), tv(len, nd), tvp(len, nd)
400
401  ! outputs:
402  INTEGER iflag(len) ! also an input
403
404  ! local variables:
405  INTEGER i
406
407  DO i = 1, len
408    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
409            icb(i))<=(tv(i, icb(i)) - dtmax))) iflag(i) = 4
410  END DO
411
412END SUBROUTINE cv_trigger
413
414SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
415        tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
416        tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
417        v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
418  USE lmdz_print_control, ONLY: lunout
419  USE lmdz_abort_physic, ONLY: abort_physic
420  IMPLICIT NONE
421
422 
423  ! inputs:
424  INTEGER len, ncum, nd, nloc
425  INTEGER iflag1(len), nk1(len), icb1(len)
426  REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
427  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
428  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
429  REAL p1(len, nd), ph1(len, nd + 1), tv1(len, nd), tp1(len, nd)
430  REAL tvp1(len, nd), clw1(len, nd)
431
432  ! outputs:
433  INTEGER iflag(nloc), nk(nloc), icb(nloc)
434  REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
435  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
436  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
437  REAL p(nloc, nd), ph(nloc, nd + 1), tv(nloc, nd), tp(nloc, nd)
438  REAL tvp(nloc, nd), clw(nloc, nd)
439  REAL dph(nloc, nd)
440
441  ! local variables:
442  INTEGER i, k, nn
443  CHARACTER (LEN = 20) :: modname = 'cv_compress'
444  CHARACTER (LEN = 80) :: abort_message
445
446  DO k = 1, nl + 1
447    nn = 0
448    DO i = 1, len
449      IF (iflag1(i)==0) THEN
450        nn = nn + 1
451        t(nn, k) = t1(i, k)
452        q(nn, k) = q1(i, k)
453        qs(nn, k) = qs1(i, k)
454        u(nn, k) = u1(i, k)
455        v(nn, k) = v1(i, k)
456        gz(nn, k) = gz1(i, k)
457        h(nn, k) = h1(i, k)
458        lv(nn, k) = lv1(i, k)
459        cpn(nn, k) = cpn1(i, k)
460        p(nn, k) = p1(i, k)
461        ph(nn, k) = ph1(i, k)
462        tv(nn, k) = tv1(i, k)
463        tp(nn, k) = tp1(i, k)
464        tvp(nn, k) = tvp1(i, k)
465        clw(nn, k) = clw1(i, k)
466      END IF
467    END DO
468  END DO
469
470  IF (nn/=ncum) THEN
471    WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
472    abort_message = ''
473    CALL abort_physic(modname, abort_message, 1)
474  END IF
475
476  nn = 0
477  DO i = 1, len
478    IF (iflag1(i)==0) THEN
479      nn = nn + 1
480      cbmf(nn) = cbmf1(i)
481      plcl(nn) = plcl1(i)
482      tnk(nn) = tnk1(i)
483      qnk(nn) = qnk1(i)
484      gznk(nn) = gznk1(i)
485      nk(nn) = nk1(i)
486      icb(nn) = icb1(i)
487      iflag(nn) = iflag1(i)
488    END IF
489  END DO
490
491  DO k = 1, nl
492    DO i = 1, ncum
493      dph(i, k) = ph(i, k) - ph(i, k + 1)
494    END DO
495  END DO
496
497END SUBROUTINE cv_compress
498
499SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
500        gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
501  USE lmdz_cvthermo
502
503  IMPLICIT NONE
504
505  ! ---------------------------------------------------------------------
506  ! Purpose:
507  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
508  ! &
509  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
510  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
511  ! &
512  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
513  ! ---------------------------------------------------------------------
514
515 
516  ! inputs:
517  INTEGER ncum, nd, nloc
518  INTEGER icb(nloc), nk(nloc)
519  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
520  REAL p(nloc, nd), dph(nloc, nd)
521  REAL tnk(nloc), qnk(nloc), gznk(nloc)
522  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
523
524  ! outputs:
525  INTEGER inb(nloc), inb1(nloc)
526  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
527  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
528  REAL frac(nloc)
529
530  ! local variables:
531  INTEGER i, k
532  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
533  REAL by, defrac
534  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
535  LOGICAL lcape(nloc)
536
537  ! =====================================================================
538  ! --- SOME INITIALIZATIONS
539  ! =====================================================================
540
541  DO k = 1, nl
542    DO i = 1, ncum
543      ep(i, k) = 0.0
544      sigp(i, k) = sigs
545    END DO
546  END DO
547
548  ! =====================================================================
549  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
550  ! =====================================================================
551
552  ! ---       The procedure is to solve the equation.
553  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
554
555  ! ***  Calculate certain parcel quantities, including static energy   ***
556
557  DO i = 1, ncum
558    ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
559            t0)) + gznk(i)
560  END DO
561
562
563  ! ***  Find lifted parcel quantities above cloud base    ***
564
565  DO k = minorig + 1, nl
566    DO i = 1, ncum
567      IF (k>=(icb(i) + 1)) THEN
568        tg = t(i, k)
569        qg = qs(i, k)
570        alv = lv0 - clmcpv * (t(i, k) - t0)
571
572        ! First iteration.
573
574        s = cpd + alv * alv * qg / (rrv * t(i, k) * t(i, k))
575        s = 1. / s
576        ahg = cpd * tg + (cl - cpd) * qnk(i) * t(i, k) + alv * qg + gz(i, k)
577        tg = tg + s * (ah0(i) - ahg)
578        tg = max(tg, 35.0)
579        tc = tg - t0
580        denom = 243.5 + tc
581        IF (tc>=0.0) THEN
582          es = 6.112 * exp(17.67 * tc / denom)
583        ELSE
584          es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
585        END IF
586        qg = eps * es / (p(i, k) - es * (1. - eps))
587
588        ! Second iteration.
589
590        s = cpd + alv * alv * qg / (rrv * t(i, k) * t(i, k))
591        s = 1. / s
592        ahg = cpd * tg + (cl - cpd) * qnk(i) * t(i, k) + alv * qg + gz(i, k)
593        tg = tg + s * (ah0(i) - ahg)
594        tg = max(tg, 35.0)
595        tc = tg - t0
596        denom = 243.5 + tc
597        IF (tc>=0.0) THEN
598          es = 6.112 * exp(17.67 * tc / denom)
599        ELSE
600          es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
601        END IF
602        qg = eps * es / (p(i, k) - es * (1. - eps))
603
604        alv = lv0 - clmcpv * (t(i, k) - t0)
605        ! PRINT*,'cpd dans convect2 ',cpd
606        ! PRINT*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
607        ! PRINT*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
608        tp(i, k) = (ah0(i) - (cl - cpd) * qnk(i) * t(i, k) - gz(i, k) - alv * qg) / cpd
609        ! if (.NOT.cpd.gt.1000.) THEN
610        ! PRINT*,'CPD=',cpd
611        ! stop
612        ! END IF
613        clw(i, k) = qnk(i) - qg
614        clw(i, k) = max(0.0, clw(i, k))
615        rg = qg / (1. - qnk(i))
616        tvp(i, k) = tp(i, k) * (1. + rg * epsi)
617      END IF
618    END DO
619  END DO
620
621  ! =====================================================================
622  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
623  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
624  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
625  ! =====================================================================
626
627  DO k = minorig + 1, nl
628    DO i = 1, ncum
629      IF (k>=(nk(i) + 1)) THEN
630        tca = tp(i, k) - t0
631        IF (tca>=0.0) THEN
632          elacrit = elcrit
633        ELSE
634          elacrit = elcrit * (1.0 - tca / tlcrit)
635        END IF
636        elacrit = max(elacrit, 0.0)
637        ep(i, k) = 1.0 - elacrit / max(clw(i, k), 1.0E-8)
638        ep(i, k) = max(ep(i, k), 0.0)
639        ep(i, k) = min(ep(i, k), 1.0)
640        sigp(i, k) = sigs
641      END IF
642    END DO
643  END DO
644
645  ! =====================================================================
646  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
647  ! --- VIRTUAL TEMPERATURE
648  ! =====================================================================
649
650  DO k = minorig + 1, nl
651    DO i = 1, ncum
652      IF (k>=(icb(i) + 1)) THEN
653        tvp(i, k) = tvp(i, k) * (1.0 - qnk(i) + ep(i, k) * clw(i, k))
654        ! PRINT*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
655        ! PRINT*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
656      END IF
657    END DO
658  END DO
659  DO i = 1, ncum
660    tvp(i, nlp) = tvp(i, nl) - (gz(i, nlp) - gz(i, nl)) / cpd
661  END DO
662
663  ! =====================================================================
664  ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
665  ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
666  ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
667  ! =====================================================================
668
669  DO i = 1, ncum
670    cape(i) = 0.0
671    capem(i) = 0.0
672    inb(i) = icb(i) + 1
673    inb1(i) = inb(i)
674  END DO
675
676  ! Originial Code
677
678  ! do 530 k=minorig+1,nl-1
679  ! do 520 i=1,ncum
680  ! IF(k.ge.(icb(i)+1))THEN
681  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
682  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
683  ! cape(i)=cape(i)+by
684  ! IF(by.ge.0.0)inb1(i)=k+1
685  ! IF(cape(i).gt.0.0)THEN
686  ! inb(i)=k+1
687  ! capem(i)=cape(i)
688  ! END IF
689  ! END IF
690  ! 520    continue
691  ! 530  continue
692  ! do 540 i=1,ncum
693  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
694  ! cape(i)=capem(i)+byp
695  ! defrac=capem(i)-cape(i)
696  ! defrac=max(defrac,0.001)
697  ! frac(i)=-cape(i)/defrac
698  ! frac(i)=min(frac(i),1.0)
699  ! frac(i)=max(frac(i),0.0)
700  ! 540   continue
701
702  ! K Emanuel fix
703
704  ! CALL zilch(byp,ncum)
705  ! do 530 k=minorig+1,nl-1
706  ! do 520 i=1,ncum
707  ! IF(k.ge.(icb(i)+1))THEN
708  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
709  ! cape(i)=cape(i)+by
710  ! IF(by.ge.0.0)inb1(i)=k+1
711  ! IF(cape(i).gt.0.0)THEN
712  ! inb(i)=k+1
713  ! capem(i)=cape(i)
714  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
715  ! END IF
716  ! END IF
717  ! 520    continue
718  ! 530  continue
719  ! do 540 i=1,ncum
720  ! inb(i)=max(inb(i),inb1(i))
721  ! cape(i)=capem(i)+byp(i)
722  ! defrac=capem(i)-cape(i)
723  ! defrac=max(defrac,0.001)
724  ! frac(i)=-cape(i)/defrac
725  ! frac(i)=min(frac(i),1.0)
726  ! frac(i)=max(frac(i),0.0)
727  ! 540   continue
728
729  ! J Teixeira fix
730
731  CALL zilch(byp, ncum)
732  DO i = 1, ncum
733    lcape(i) = .TRUE.
734  END DO
735  DO k = minorig + 1, nl - 1
736    DO i = 1, ncum
737      IF (cape(i)<0.0) lcape(i) = .FALSE.
738      IF ((k>=(icb(i) + 1)) .AND. lcape(i)) THEN
739        by = (tvp(i, k) - tv(i, k)) * dph(i, k) / p(i, k)
740        byp(i) = (tvp(i, k + 1) - tv(i, k + 1)) * dph(i, k + 1) / p(i, k + 1)
741        cape(i) = cape(i) + by
742        IF (by>=0.0) inb1(i) = k + 1
743        IF (cape(i)>0.0) THEN
744          inb(i) = k + 1
745          capem(i) = cape(i)
746        END IF
747      END IF
748    END DO
749  END DO
750  DO i = 1, ncum
751    cape(i) = capem(i) + byp(i)
752    defrac = capem(i) - cape(i)
753    defrac = max(defrac, 0.001)
754    frac(i) = -cape(i) / defrac
755    frac(i) = min(frac(i), 1.0)
756    frac(i) = max(frac(i), 0.0)
757  END DO
758
759  ! =====================================================================
760  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
761  ! =====================================================================
762
763  ! initialization:
764  DO i = 1, ncum * nlp
765    hp(i, 1) = h(i, 1)
766  END DO
767
768  DO k = minorig + 1, nl
769    DO i = 1, ncum
770      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
771        hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k &
772                )
773      END IF
774    END DO
775  END DO
776
777END SUBROUTINE cv_undilute2
778
779SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
780        cpn, iflag, cbmf)
781  USE lmdz_cvthermo
782
783  IMPLICIT NONE
784
785  ! inputs:
786  INTEGER ncum, nd, nloc
787  INTEGER nk(nloc), icb(nloc)
788  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
789  REAL ph(nloc, nd + 1) ! caution nd instead ndp1 to be consistent...
790  REAL plcl(nloc), cpn(nloc, nd)
791
792  ! outputs:
793  INTEGER iflag(nloc)
794  REAL cbmf(nloc) ! also an input
795
796  ! local variables:
797  INTEGER i, k, icbmax
798  REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
799  REAL work(nloc)
800
801 
802  ! -------------------------------------------------------------------
803  ! Compute icbmax.
804  ! -------------------------------------------------------------------
805
806  icbmax = 2
807  DO i = 1, ncum
808    icbmax = max(icbmax, icb(i))
809  END DO
810
811  ! =====================================================================
812  ! ---  CALCULATE CLOUD BASE MASS FLUX
813  ! =====================================================================
814
815  ! tvpplcl = parcel temperature lifted adiabatically from level
816  ! icb-1 to the LCL.
817  ! tvaplcl = virtual temperature at the LCL.
818
819  DO i = 1, ncum
820    dtpbl(i) = 0.0
821    tvpplcl(i) = tvp(i, icb(i) - 1) - rrd * tvp(i, icb(i) - 1) * (p(i, icb(i) - 1) - plcl(&
822            i)) / (cpn(i, icb(i) - 1) * p(i, icb(i) - 1))
823    tvaplcl(i) = tv(i, icb(i)) + (tvp(i, icb(i)) - tvp(i, icb(i) + 1)) * (plcl(i) - p(i &
824            , icb(i))) / (p(i, icb(i)) - p(i, icb(i) + 1))
825  END DO
826
827  ! -------------------------------------------------------------------
828  ! --- Interpolate difference between lifted parcel and
829  ! --- environmental temperatures to lifted condensation level
830  ! -------------------------------------------------------------------
831
832  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
833
834  DO k = minorig, icbmax
835    DO i = 1, ncum
836      IF ((k>=nk(i)) .AND. (k<=(icb(i) - 1))) THEN
837        dtpbl(i) = dtpbl(i) + (tvp(i, k) - tv(i, k)) * dph(i, k)
838      END IF
839    END DO
840  END DO
841  DO i = 1, ncum
842    dtpbl(i) = dtpbl(i) / (ph(i, nk(i)) - ph(i, icb(i)))
843    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
844  END DO
845
846  ! -------------------------------------------------------------------
847  ! --- Adjust cloud base mass flux
848  ! -------------------------------------------------------------------
849
850  DO i = 1, ncum
851    work(i) = cbmf(i)
852    cbmf(i) = max(0.0, (1.0 - damp) * cbmf(i) + 0.1 * alpha * dtmin(i))
853    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
854      iflag(i) = 3
855    END IF
856  END DO
857
858END SUBROUTINE cv_closure
859
860SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
861        h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
862        sij, elij)
863  USE lmdz_cvthermo
864
865  IMPLICIT NONE
866
867 
868  ! inputs:
869  INTEGER ncum, nd, nloc
870  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
871  REAL cbmf(nloc), qnk(nloc)
872  REAL ph(nloc, nd + 1)
873  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
874  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
875  REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
876
877  ! outputs:
878  INTEGER nent(nloc, nd)
879  REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
880  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
881  REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
882
883  ! local variables:
884  INTEGER i, j, k, ij
885  INTEGER num1, num2
886  REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
887  REAL alt, qp1, smid, sjmin, sjmax, delp, delm
888  REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
889  REAL bsum(nloc, nd)
890  LOGICAL lwork(nloc)
891
892  ! =====================================================================
893  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
894  ! =====================================================================
895
896  DO i = 1, ncum * nlp
897    nent(i, 1) = 0
898    m(i, 1) = 0.0
899  END DO
900
901  DO k = 1, nlp
902    DO j = 1, nlp
903      DO i = 1, ncum
904        qent(i, k, j) = q(i, j)
905        uent(i, k, j) = u(i, j)
906        vent(i, k, j) = v(i, j)
907        elij(i, k, j) = 0.0
908        ment(i, k, j) = 0.0
909        sij(i, k, j) = 0.0
910      END DO
911    END DO
912  END DO
913
914  ! -------------------------------------------------------------------
915  ! --- Calculate rates of mixing,  m(i)
916  ! -------------------------------------------------------------------
917
918  CALL zilch(work, ncum)
919
920  DO j = minorig + 1, nl
921    DO i = 1, ncum
922      IF ((j>=(icb(i) + 1)) .AND. (j<=inb(i))) THEN
923        k = min(j, inb1(i))
924        dbo = abs(tv(i, k + 1) - tvp(i, k + 1) - tv(i, k - 1) + tvp(i, k - 1)) + &
925                entp * 0.04 * (ph(i, k) - ph(i, k + 1))
926        work(i) = work(i) + dbo
927        m(i, j) = cbmf(i) * dbo
928      END IF
929    END DO
930  END DO
931  DO k = minorig + 1, nl
932    DO i = 1, ncum
933      IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i))) THEN
934        m(i, k) = m(i, k) / work(i)
935      END IF
936    END DO
937  END DO
938
939
940  ! =====================================================================
941  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
942  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
943  ! --- FRACTION (sij)
944  ! =====================================================================
945
946  DO i = minorig + 1, nl
947    DO j = minorig + 1, nl
948      DO ij = 1, ncum
949        IF ((i>=(icb(ij) + 1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
950                inb(ij))) THEN
951          qti = qnk(ij) - ep(ij, i) * clw(ij, i)
952          bf2 = 1. + lv(ij, j) * lv(ij, j) * qs(ij, j) / (rrv * t(ij, j) * t(ij, j) * cpd)
953          anum = h(ij, j) - hp(ij, i) + (cpv - cpd) * t(ij, j) * (qti - q(ij, j))
954          denom = h(ij, i) - hp(ij, i) + (cpd - cpv) * (q(ij, i) - qti) * t(ij, j)
955          dei = denom
956          IF (abs(dei)<0.01) dei = 0.01
957          sij(ij, i, j) = anum / dei
958          sij(ij, i, i) = 1.0
959          altem = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti - qs(ij, j)
960          altem = altem / bf2
961          cwat = clw(ij, j) * (1. - ep(ij, j))
962          stemp = sij(ij, i, j)
963          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
964            anum = anum - lv(ij, j) * (qti - qs(ij, j) - cwat * bf2)
965            denom = denom + lv(ij, j) * (q(ij, i) - qti)
966            IF (abs(denom)<0.01) denom = 0.01
967            sij(ij, i, j) = anum / denom
968            altem = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti - qs(ij, j)
969            altem = altem - (bf2 - 1.) * cwat
970          END IF
971          IF (sij(ij, i, j)>0.0 .AND. sij(ij, i, j)<0.9) THEN
972            qent(ij, i, j) = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti
973            uent(ij, i, j) = sij(ij, i, j) * u(ij, i) + &
974                    (1. - sij(ij, i, j)) * u(ij, nk(ij))
975            vent(ij, i, j) = sij(ij, i, j) * v(ij, i) + &
976                    (1. - sij(ij, i, j)) * v(ij, nk(ij))
977            elij(ij, i, j) = altem
978            elij(ij, i, j) = max(0.0, elij(ij, i, j))
979            ment(ij, i, j) = m(ij, i) / (1. - sij(ij, i, j))
980            nent(ij, i) = nent(ij, i) + 1
981          END IF
982          sij(ij, i, j) = max(0.0, sij(ij, i, j))
983          sij(ij, i, j) = min(1.0, sij(ij, i, j))
984        END IF
985      END DO
986    END DO
987
988    ! ***   If no air can entrain at level i assume that updraft detrains
989    ! ***
990    ! ***   at that level and calculate detrained air flux and properties
991    ! ***
992
993    DO ij = 1, ncum
994      IF ((i>=(icb(ij) + 1)) .AND. (i<=inb(ij)) .AND. (nent(ij, i)==0)) THEN
995        ment(ij, i, i) = m(ij, i)
996        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
997        uent(ij, i, i) = u(ij, nk(ij))
998        vent(ij, i, i) = v(ij, nk(ij))
999        elij(ij, i, i) = clw(ij, i)
1000        sij(ij, i, i) = 1.0
1001      END IF
1002    END DO
1003  END DO
1004
1005  DO i = 1, ncum
1006    sij(i, inb(i), inb(i)) = 1.0
1007  END DO
1008
1009  ! =====================================================================
1010  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
1011  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
1012  ! =====================================================================
1013
1014  CALL zilch(bsum, ncum * nlp)
1015  DO ij = 1, ncum
1016    lwork(ij) = .FALSE.
1017  END DO
1018  DO i = minorig + 1, nl
1019
1020    num1 = 0
1021    DO ij = 1, ncum
1022      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) num1 = num1 + 1
1023    END DO
1024    IF (num1<=0) GO TO 789
1025
1026    DO ij = 1, ncum
1027      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) THEN
1028        lwork(ij) = (nent(ij, i)/=0)
1029        qp1 = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
1030        anum = h(ij, i) - hp(ij, i) - lv(ij, i) * (qp1 - qs(ij, i))
1031        denom = h(ij, i) - hp(ij, i) + lv(ij, i) * (q(ij, i) - qp1)
1032        IF (abs(denom)<0.01) denom = 0.01
1033        scrit(ij) = anum / denom
1034        alt = qp1 - qs(ij, i) + scrit(ij) * (q(ij, i) - qp1)
1035        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
1036        asij(ij) = 0.0
1037        smin(ij) = 1.0
1038      END IF
1039    END DO
1040    DO j = minorig, nl
1041
1042      num2 = 0
1043      DO ij = 1, ncum
1044        IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
1045                ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
1046      END DO
1047      IF (num2<=0) GO TO 783
1048
1049      DO ij = 1, ncum
1050        IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
1051                ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1052          IF (sij(ij, i, j)>0.0 .AND. sij(ij, i, j)<0.9) THEN
1053            IF (j>i) THEN
1054              smid = min(sij(ij, i, j), scrit(ij))
1055              sjmax = smid
1056              sjmin = smid
1057              IF (smid<smin(ij) .AND. sij(ij, i, j + 1)<smid) THEN
1058                smin(ij) = smid
1059                sjmax = min(sij(ij, i, j + 1), sij(ij, i, j), scrit(ij))
1060                sjmin = max(sij(ij, i, j - 1), sij(ij, i, j))
1061                sjmin = min(sjmin, scrit(ij))
1062              END IF
1063            ELSE
1064              sjmax = max(sij(ij, i, j + 1), scrit(ij))
1065              smid = max(sij(ij, i, j), scrit(ij))
1066              sjmin = 0.0
1067              IF (j>1) sjmin = sij(ij, i, j - 1)
1068              sjmin = max(sjmin, scrit(ij))
1069            END IF
1070            delp = abs(sjmax - smid)
1071            delm = abs(sjmin - smid)
1072            asij(ij) = asij(ij) + (delp + delm) * (ph(ij, j) - ph(ij, j + 1))
1073            ment(ij, i, j) = ment(ij, i, j) * (delp + delm) * (ph(ij, j) - ph(ij, j + 1))
1074          END IF
1075        END IF
1076      END DO
1077    783 END DO
1078    DO ij = 1, ncum
1079      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
1080        asij(ij) = max(1.0E-21, asij(ij))
1081        asij(ij) = 1.0 / asij(ij)
1082        bsum(ij, i) = 0.0
1083      END IF
1084    END DO
1085    DO j = minorig, nl + 1
1086      DO ij = 1, ncum
1087        IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
1088                ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1089          ment(ij, i, j) = ment(ij, i, j) * asij(ij)
1090          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
1091        END IF
1092      END DO
1093    END DO
1094    DO ij = 1, ncum
1095      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
1096              i)<1.0E-18) .AND. lwork(ij)) THEN
1097        nent(ij, i) = 0
1098        ment(ij, i, i) = m(ij, i)
1099        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
1100        uent(ij, i, i) = u(ij, nk(ij))
1101        vent(ij, i, i) = v(ij, nk(ij))
1102        elij(ij, i, i) = clw(ij, i)
1103        sij(ij, i, i) = 1.0
1104      END IF
1105    END DO
1106  789 END DO
1107
1108END SUBROUTINE cv_mixing
1109
1110SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
1111        ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
1112  USE lmdz_cvthermo
1113
1114  IMPLICIT NONE
1115
1116 
1117  ! inputs:
1118  INTEGER ncum, nd, nloc
1119  INTEGER inb(nloc)
1120  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
1121  REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
1122  REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
1123  REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
1124  REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
1125
1126  ! outputs:
1127  INTEGER iflag(nloc) ! also an input
1128  REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd)
1129  REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd)
1130
1131  ! local variables:
1132  INTEGER i, j, k, ij, num1
1133  INTEGER jtt(nloc)
1134  REAL awat, coeff, qsm, afac, sigt, b6, c6, revap
1135  REAL dhdp, fac, qstm, rat
1136  REAL wdtrain(nloc)
1137  LOGICAL lwork(nloc)
1138
1139  ! =====================================================================
1140  ! --- PRECIPITATING DOWNDRAFT CALCULATION
1141  ! =====================================================================
1142
1143  ! Initializations:
1144
1145  DO i = 1, ncum
1146    DO k = 1, nl + 1
1147      wt(i, k) = omtsnow
1148      mp(i, k) = 0.0
1149      evap(i, k) = 0.0
1150      water(i, k) = 0.0
1151    END DO
1152  END DO
1153
1154  DO i = 1, ncum
1155    qp(i, 1) = q(i, 1)
1156    up(i, 1) = u(i, 1)
1157    vp(i, 1) = v(i, 1)
1158  END DO
1159
1160  DO k = 2, nl + 1
1161    DO i = 1, ncum
1162      qp(i, k) = q(i, k - 1)
1163      up(i, k) = u(i, k - 1)
1164      vp(i, k) = v(i, k - 1)
1165    END DO
1166  END DO
1167
1168
1169  ! ***  Check whether ep(inb)=0, if so, skip precipitating    ***
1170  ! ***             downdraft calculation                      ***
1171
1172
1173  ! ***  Integrate liquid water equation to find condensed water   ***
1174  ! ***                and condensed water flux                    ***
1175
1176  DO i = 1, ncum
1177    jtt(i) = 2
1178    IF (ep(i, inb(i))<=0.0001) iflag(i) = 2
1179    IF (iflag(i)==0) THEN
1180      lwork(i) = .TRUE.
1181    ELSE
1182      lwork(i) = .FALSE.
1183    END IF
1184  END DO
1185
1186  ! ***                    Begin downdraft loop                    ***
1187
1188  CALL zilch(wdtrain, ncum)
1189  DO i = nl + 1, 1, -1
1190
1191    num1 = 0
1192    DO ij = 1, ncum
1193      IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
1194    END DO
1195    IF (num1<=0) GO TO 899
1196
1197
1198    ! ***        Calculate detrained precipitation             ***
1199
1200    DO ij = 1, ncum
1201      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1202        wdtrain(ij) = g * ep(ij, i) * m(ij, i) * clw(ij, i)
1203      END IF
1204    END DO
1205
1206    IF (i>1) THEN
1207      DO j = 1, i - 1
1208        DO ij = 1, ncum
1209          IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1210            awat = elij(ij, j, i) - (1. - ep(ij, i)) * clw(ij, i)
1211            awat = max(0.0, awat)
1212            wdtrain(ij) = wdtrain(ij) + g * awat * ment(ij, j, i)
1213          END IF
1214        END DO
1215      END DO
1216    END IF
1217
1218    ! ***    Find rain water and evaporation using provisional   ***
1219    ! ***              estimates of qp(i)and qp(i-1)             ***
1220
1221
1222    ! ***  Value of terminal velocity and coeffecient of evaporation for snow
1223    ! ***
1224
1225    DO ij = 1, ncum
1226      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1227        coeff = coeffs
1228        wt(ij, i) = omtsnow
1229
1230        ! ***  Value of terminal velocity and coeffecient of evaporation for
1231        ! rain   ***
1232
1233        IF (t(ij, i)>273.0) THEN
1234          coeff = coeffr
1235          wt(ij, i) = omtrain
1236        END IF
1237        qsm = 0.5 * (q(ij, i) + qp(ij, i + 1))
1238        afac = coeff * ph(ij, i) * (qs(ij, i) - qsm) / (1.0E4 + 2.0E3 * ph(ij, i) * qs(ij, i))
1239        afac = max(afac, 0.0)
1240        sigt = sigp(ij, i)
1241        sigt = max(0.0, sigt)
1242        sigt = min(1.0, sigt)
1243        b6 = 100. * (ph(ij, i) - ph(ij, i + 1)) * sigt * afac / wt(ij, i)
1244        c6 = (water(ij, i + 1) * wt(ij, i + 1) + wdtrain(ij) / sigd) / wt(ij, i)
1245        revap = 0.5 * (-b6 + sqrt(b6 * b6 + 4. * c6))
1246        evap(ij, i) = sigt * afac * revap
1247        water(ij, i) = revap * revap
1248
1249        ! ***  Calculate precipitating downdraft mass flux under     ***
1250        ! ***              hydrostatic approximation                 ***
1251
1252        IF (i>1) THEN
1253          dhdp = (h(ij, i) - h(ij, i - 1)) / (p(ij, i - 1) - p(ij, i))
1254          dhdp = max(dhdp, 10.0)
1255          mp(ij, i) = 100. * ginv * lv(ij, i) * sigd * evap(ij, i) / dhdp
1256          mp(ij, i) = max(mp(ij, i), 0.0)
1257
1258          ! ***   Add small amount of inertia to downdraft              ***
1259
1260          fac = 20.0 / (ph(ij, i - 1) - ph(ij, i))
1261          mp(ij, i) = (fac * mp(ij, i + 1) + mp(ij, i)) / (1. + fac)
1262
1263          ! ***      Force mp to decrease linearly to zero
1264          ! ***
1265          ! ***      between about 950 mb and the surface
1266          ! ***
1267
1268          IF (p(ij, i)>(0.949 * p(ij, 1))) THEN
1269            jtt(ij) = max(jtt(ij), i)
1270            mp(ij, i) = mp(ij, jtt(ij)) * (p(ij, 1) - p(ij, i)) / &
1271                    (p(ij, 1) - p(ij, jtt(ij)))
1272          END IF
1273        END IF
1274
1275        ! ***       Find mixing ratio of precipitating downdraft     ***
1276
1277        IF (i/=inb(ij)) THEN
1278          IF (i==1) THEN
1279            qstm = qs(ij, 1)
1280          ELSE
1281            qstm = qs(ij, i - 1)
1282          END IF
1283          IF (mp(ij, i)>mp(ij, i + 1)) THEN
1284            rat = mp(ij, i + 1) / mp(ij, i)
1285            qp(ij, i) = qp(ij, i + 1) * rat + q(ij, i) * (1.0 - rat) + &
1286                    100. * ginv * sigd * (ph(ij, i) - ph(ij, i + 1)) * (evap(ij, i) / mp(ij, i))
1287            up(ij, i) = up(ij, i + 1) * rat + u(ij, i) * (1. - rat)
1288            vp(ij, i) = vp(ij, i + 1) * rat + v(ij, i) * (1. - rat)
1289          ELSE
1290            IF (mp(ij, i + 1)>0.0) THEN
1291              qp(ij, i) = (gz(ij, i + 1) - gz(ij, i) + qp(ij, i + 1) * (lv(ij, i + 1) + t(ij, &
1292                      i + 1) * (cl - cpd)) + cpd * (t(ij, i + 1) - t(ij, &
1293                      i))) / (lv(ij, i) + t(ij, i) * (cl - cpd))
1294              up(ij, i) = up(ij, i + 1)
1295              vp(ij, i) = vp(ij, i + 1)
1296            END IF
1297          END IF
1298          qp(ij, i) = min(qp(ij, i), qstm)
1299          qp(ij, i) = max(qp(ij, i), 0.0)
1300        END IF
1301      END IF
1302    END DO
1303  899 END DO
1304
1305END SUBROUTINE cv_unsat
1306
1307SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
1308        ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
1309        ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
1310        precip, cbmf, ft, fq, fu, fv, ma, qcondc)
1311  USE lmdz_cvthermo
1312
1313  IMPLICIT NONE
1314
1315 
1316  ! inputs
1317  INTEGER ncum, nd, nloc
1318  INTEGER nk(nloc), icb(nloc), inb(nloc)
1319  INTEGER nent(nloc, nd)
1320  REAL delt
1321  REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
1322  REAL gz(nloc, nd)
1323  REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
1324  REAL hp(nloc, nd), lv(nloc, nd)
1325  REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
1326  REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd)
1327  REAL up(nloc, nd), vp(nloc, nd)
1328  REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd)
1329  REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd)
1330  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
1331  REAL tv(nloc, nd), tvp(nloc, nd)
1332
1333  ! outputs
1334  INTEGER iflag(nloc) ! also an input
1335  REAL cbmf(nloc) ! also an input
1336  REAL wd(nloc), tprime(nloc), qprime(nloc)
1337  REAL precip(nloc)
1338  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1339  REAL ma(nloc, nd)
1340  REAL qcondc(nloc, nd)
1341
1342  ! local variables
1343  INTEGER i, j, ij, k, num1
1344  REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti
1345  REAL work(nloc), am(nloc), amp1(nloc), ad(nloc)
1346  REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd)
1347  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
1348  REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld
1349
1350
1351  ! -- initializations:
1352
1353  delti = 1.0 / delt
1354
1355  DO i = 1, ncum
1356    precip(i) = 0.0
1357    wd(i) = 0.0
1358    tprime(i) = 0.0
1359    qprime(i) = 0.0
1360    DO k = 1, nl + 1
1361      ft(i, k) = 0.0
1362      fu(i, k) = 0.0
1363      fv(i, k) = 0.0
1364      fq(i, k) = 0.0
1365      lvcp(i, k) = lv(i, k) / cpn(i, k)
1366      qcondc(i, k) = 0.0 ! cld
1367      qcond(i, k) = 0.0 ! cld
1368      nqcond(i, k) = 0.0 ! cld
1369    END DO
1370  END DO
1371
1372
1373  ! ***  Calculate surface precipitation in mm/day     ***
1374
1375  DO i = 1, ncum
1376    IF (iflag(i)<=1) THEN
1377      ! c            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1378      ! c     &                /(rowl*g)
1379      ! c            precip(i)=precip(i)*delt/86400.
1380      precip(i) = wt(i, 1) * sigd * water(i, 1) * 86400 / g
1381    END IF
1382  END DO
1383
1384
1385  ! ***  Calculate downdraft velocity scale and surface temperature and  ***
1386  ! ***                    water vapor fluctuations                      ***
1387
1388  DO i = 1, ncum
1389    wd(i) = betad * abs(mp(i, icb(i))) * 0.01 * rrd * t(i, icb(i)) / (sigd * p(i, icb(i)))
1390    qprime(i) = 0.5 * (qp(i, 1) - q(i, 1))
1391    tprime(i) = lv0 * qprime(i) / cpd
1392  END DO
1393
1394  ! ***  Calculate tendencies of lowest level potential temperature  ***
1395  ! ***                      and mixing ratio                        ***
1396
1397  DO i = 1, ncum
1398    work(i) = 0.01 / (ph(i, 1) - ph(i, 2))
1399    am(i) = 0.0
1400  END DO
1401  DO k = 2, nl
1402    DO i = 1, ncum
1403      IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
1404        am(i) = am(i) + m(i, k)
1405      END IF
1406    END DO
1407  END DO
1408  DO i = 1, ncum
1409    IF ((g * work(i) * am(i))>=delti) iflag(i) = 1
1410    ft(i, 1) = ft(i, 1) + g * work(i) * am(i) * (t(i, 2) - t(i, 1) + (gz(i, 2) - gz(i, &
1411            1)) / cpn(i, 1))
1412    ft(i, 1) = ft(i, 1) - lvcp(i, 1) * sigd * evap(i, 1)
1413    ft(i, 1) = ft(i, 1) + sigd * wt(i, 2) * (cl - cpd) * water(i, 2) * (t(i, 2) - t(i, 1)) * &
1414            work(i) / cpn(i, 1)
1415    fq(i, 1) = fq(i, 1) + g * mp(i, 2) * (qp(i, 2) - q(i, 1)) * work(i) + &
1416            sigd * evap(i, 1)
1417    fq(i, 1) = fq(i, 1) + g * am(i) * (q(i, 2) - q(i, 1)) * work(i)
1418    fu(i, 1) = fu(i, 1) + g * work(i) * (mp(i, 2) * (up(i, 2) - u(i, 1)) + am(i) * (u(i, &
1419            2) - u(i, 1)))
1420    fv(i, 1) = fv(i, 1) + g * work(i) * (mp(i, 2) * (vp(i, 2) - v(i, 1)) + am(i) * (v(i, &
1421            2) - v(i, 1)))
1422  END DO
1423  DO j = 2, nl
1424    DO i = 1, ncum
1425      IF (j<=inb(i)) THEN
1426        fq(i, 1) = fq(i, 1) + g * work(i) * ment(i, j, 1) * (qent(i, j, 1) - q(i, 1))
1427        fu(i, 1) = fu(i, 1) + g * work(i) * ment(i, j, 1) * (uent(i, j, 1) - u(i, 1))
1428        fv(i, 1) = fv(i, 1) + g * work(i) * ment(i, j, 1) * (vent(i, j, 1) - v(i, 1))
1429      END IF
1430    END DO
1431  END DO
1432
1433  ! ***  Calculate tendencies of potential temperature and mixing ratio  ***
1434  ! ***               at levels above the lowest level                   ***
1435
1436  ! ***  First find the net saturated updraft and downdraft mass fluxes  ***
1437  ! ***                      through each level                          ***
1438
1439  DO i = 2, nl + 1
1440
1441    num1 = 0
1442    DO ij = 1, ncum
1443      IF (i<=inb(ij)) num1 = num1 + 1
1444    END DO
1445    IF (num1<=0) GO TO 1500
1446
1447    CALL zilch(amp1, ncum)
1448    CALL zilch(ad, ncum)
1449
1450    DO k = i + 1, nl + 1
1451      DO ij = 1, ncum
1452        IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij) + 1))) THEN
1453          amp1(ij) = amp1(ij) + m(ij, k)
1454        END IF
1455      END DO
1456    END DO
1457
1458    DO k = 1, i
1459      DO j = i + 1, nl + 1
1460        DO ij = 1, ncum
1461          IF ((j<=(inb(ij) + 1)) .AND. (i<=inb(ij))) THEN
1462            amp1(ij) = amp1(ij) + ment(ij, k, j)
1463          END IF
1464        END DO
1465      END DO
1466    END DO
1467    DO k = 1, i - 1
1468      DO j = i, nl + 1
1469        DO ij = 1, ncum
1470          IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
1471            ad(ij) = ad(ij) + ment(ij, j, k)
1472          END IF
1473        END DO
1474      END DO
1475    END DO
1476
1477    DO ij = 1, ncum
1478      IF (i<=inb(ij)) THEN
1479        dpinv = 0.01 / (ph(ij, i) - ph(ij, i + 1))
1480        cpinv = 1.0 / cpn(ij, i)
1481
1482        ft(ij, i) = ft(ij, i) + g * dpinv * (amp1(ij) * (t(ij, i + 1) - t(ij, &
1483                i) + (gz(ij, i + 1) - gz(ij, i)) * cpinv) - ad(ij) * (t(ij, i) - t(ij, &
1484                i - 1) + (gz(ij, i) - gz(ij, i - 1)) * cpinv)) - sigd * lvcp(ij, i) * evap(ij, i)
1485        ft(ij, i) = ft(ij, i) + g * dpinv * ment(ij, i, i) * (hp(ij, i) - h(ij, i) + t(ij &
1486                , i) * (cpv - cpd) * (q(ij, i) - qent(ij, i, i))) * cpinv
1487        ft(ij, i) = ft(ij, i) + sigd * wt(ij, i + 1) * (cl - cpd) * water(ij, i + 1) * (t(&
1488                ij, i + 1) - t(ij, i)) * dpinv * cpinv
1489        fq(ij, i) = fq(ij, i) + g * dpinv * (amp1(ij) * (q(ij, i + 1) - q(ij, &
1490                i)) - ad(ij) * (q(ij, i) - q(ij, i - 1)))
1491        fu(ij, i) = fu(ij, i) + g * dpinv * (amp1(ij) * (u(ij, i + 1) - u(ij, &
1492                i)) - ad(ij) * (u(ij, i) - u(ij, i - 1)))
1493        fv(ij, i) = fv(ij, i) + g * dpinv * (amp1(ij) * (v(ij, i + 1) - v(ij, &
1494                i)) - ad(ij) * (v(ij, i) - v(ij, i - 1)))
1495      END IF
1496    END DO
1497    DO k = 1, i - 1
1498      DO ij = 1, ncum
1499        IF (i<=inb(ij)) THEN
1500          awat = elij(ij, k, i) - (1. - ep(ij, i)) * clw(ij, i)
1501          awat = max(awat, 0.0)
1502          fq(ij, i) = fq(ij, i) + g * dpinv * ment(ij, k, i) * (qent(ij, k, i) - awat - q &
1503                  (ij, i))
1504          fu(ij, i) = fu(ij, i) + g * dpinv * ment(ij, k, i) * (uent(ij, k, i) - u(ij, i &
1505                  ))
1506          fv(ij, i) = fv(ij, i) + g * dpinv * ment(ij, k, i) * (vent(ij, k, i) - v(ij, i &
1507                  ))
1508          ! (saturated updrafts resulting from mixing)               ! cld
1509          qcond(ij, i) = qcond(ij, i) + (elij(ij, k, i) - awat) ! cld
1510          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1511        END IF
1512      END DO
1513    END DO
1514    DO k = i, nl + 1
1515      DO ij = 1, ncum
1516        IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
1517          fq(ij, i) = fq(ij, i) + g * dpinv * ment(ij, k, i) * (qent(ij, k, i) - q(ij, i &
1518                  ))
1519          fu(ij, i) = fu(ij, i) + g * dpinv * ment(ij, k, i) * (uent(ij, k, i) - u(ij, i &
1520                  ))
1521          fv(ij, i) = fv(ij, i) + g * dpinv * ment(ij, k, i) * (vent(ij, k, i) - v(ij, i &
1522                  ))
1523        END IF
1524      END DO
1525    END DO
1526    DO ij = 1, ncum
1527      IF (i<=inb(ij)) THEN
1528        fq(ij, i) = fq(ij, i) + sigd * evap(ij, i) + g * (mp(ij, i + 1) * (qp(ij, &
1529                i + 1) - q(ij, i)) - mp(ij, i) * (qp(ij, i) - q(ij, i - 1))) * dpinv
1530        fu(ij, i) = fu(ij, i) + g * (mp(ij, i + 1) * (up(ij, i + 1) - u(ij, &
1531                i)) - mp(ij, i) * (up(ij, i) - u(ij, i - 1))) * dpinv
1532        fv(ij, i) = fv(ij, i) + g * (mp(ij, i + 1) * (vp(ij, i + 1) - v(ij, &
1533                i)) - mp(ij, i) * (vp(ij, i) - v(ij, i - 1))) * dpinv
1534        ! (saturated downdrafts resulting from mixing)               ! cld
1535        DO k = i + 1, inb(ij) ! cld
1536          qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld
1537          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1538        END DO ! cld
1539        ! (particular case: no detraining level is found)            ! cld
1540        IF (nent(ij, i)==0) THEN ! cld
1541          qcond(ij, i) = qcond(ij, i) + (1. - ep(ij, i)) * clw(ij, i) ! cld
1542          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1543        END IF ! cld
1544        IF (nqcond(ij, i)/=0.) THEN ! cld
1545          qcond(ij, i) = qcond(ij, i) / nqcond(ij, i) ! cld
1546        END IF ! cld
1547      END IF
1548    END DO
1549  1500 END DO
1550
1551  ! *** Adjust tendencies at top of convection layer to reflect  ***
1552  ! ***       actual position of the level zero cape             ***
1553
1554  DO ij = 1, ncum
1555    fqold = fq(ij, inb(ij))
1556    fq(ij, inb(ij)) = fq(ij, inb(ij)) * (1. - frac(ij))
1557    fq(ij, inb(ij) - 1) = fq(ij, inb(ij) - 1) + frac(ij) * fqold * ((ph(ij, &
1558            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, &
1559            inb(ij)))) * lv(ij, inb(ij)) / lv(ij, inb(ij) - 1)
1560    ftold = ft(ij, inb(ij))
1561    ft(ij, inb(ij)) = ft(ij, inb(ij)) * (1. - frac(ij))
1562    ft(ij, inb(ij) - 1) = ft(ij, inb(ij) - 1) + frac(ij) * ftold * ((ph(ij, &
1563            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, &
1564            inb(ij)))) * cpn(ij, inb(ij)) / cpn(ij, inb(ij) - 1)
1565    fuold = fu(ij, inb(ij))
1566    fu(ij, inb(ij)) = fu(ij, inb(ij)) * (1. - frac(ij))
1567    fu(ij, inb(ij) - 1) = fu(ij, inb(ij) - 1) + frac(ij) * fuold * ((ph(ij, &
1568            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, inb(ij))))
1569    fvold = fv(ij, inb(ij))
1570    fv(ij, inb(ij)) = fv(ij, inb(ij)) * (1. - frac(ij))
1571    fv(ij, inb(ij) - 1) = fv(ij, inb(ij) - 1) + frac(ij) * fvold * ((ph(ij, &
1572            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, inb(ij))))
1573  END DO
1574
1575  ! ***   Very slightly adjust tendencies to force exact   ***
1576  ! ***     enthalpy, momentum and tracer conservation     ***
1577
1578  DO ij = 1, ncum
1579    ents(ij) = 0.0
1580    uav(ij) = 0.0
1581    vav(ij) = 0.0
1582    DO i = 1, inb(ij)
1583      ents(ij) = ents(ij) + (cpn(ij, i) * ft(ij, i) + lv(ij, i) * fq(ij, i)) * (ph(ij, i) - &
1584              ph(ij, i + 1))
1585      uav(ij) = uav(ij) + fu(ij, i) * (ph(ij, i) - ph(ij, i + 1))
1586      vav(ij) = vav(ij) + fv(ij, i) * (ph(ij, i) - ph(ij, i + 1))
1587    END DO
1588  END DO
1589  DO ij = 1, ncum
1590    ents(ij) = ents(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
1591    uav(ij) = uav(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
1592    vav(ij) = vav(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
1593  END DO
1594  DO ij = 1, ncum
1595    DO i = 1, inb(ij)
1596      ft(ij, i) = ft(ij, i) - ents(ij) / cpn(ij, i)
1597      fu(ij, i) = (1. - cu) * (fu(ij, i) - uav(ij))
1598      fv(ij, i) = (1. - cu) * (fv(ij, i) - vav(ij))
1599    END DO
1600  END DO
1601
1602  DO k = 1, nl + 1
1603    DO i = 1, ncum
1604      IF ((q(i, k) + delt * fq(i, k))<0.0) iflag(i) = 10
1605    END DO
1606  END DO
1607
1608  DO i = 1, ncum
1609    IF (iflag(i)>2) THEN
1610      precip(i) = 0.0
1611      cbmf(i) = 0.0
1612    END IF
1613  END DO
1614  DO k = 1, nl
1615    DO i = 1, ncum
1616      IF (iflag(i)>2) THEN
1617        ft(i, k) = 0.0
1618        fq(i, k) = 0.0
1619        fu(i, k) = 0.0
1620        fv(i, k) = 0.0
1621        qcondc(i, k) = 0.0 ! cld
1622      END IF
1623    END DO
1624  END DO
1625
1626  DO k = 1, nl + 1
1627    DO i = 1, ncum
1628      ma(i, k) = 0.
1629    END DO
1630  END DO
1631  DO k = nl, 1, -1
1632    DO i = 1, ncum
1633      ma(i, k) = ma(i, k + 1) + m(i, k)
1634    END DO
1635  END DO
1636
1637
1638  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
1639  ! ***           of condensed water         ***            ! cld
1640  ! cld
1641  DO ij = 1, ncum ! cld
1642    DO i = 1, nd ! cld
1643      mac(ij, i) = 0.0 ! cld
1644      wa(ij, i) = 0.0 ! cld
1645      siga(ij, i) = 0.0 ! cld
1646    END DO ! cld
1647    DO i = nk(ij), inb(ij) ! cld
1648      DO k = i + 1, inb(ij) + 1 ! cld
1649        mac(ij, i) = mac(ij, i) + m(ij, k) ! cld
1650      END DO ! cld
1651    END DO ! cld
1652    DO i = icb(ij), inb(ij) - 1 ! cld
1653      ax(ij, i) = 0. ! cld
1654      DO j = icb(ij), i ! cld
1655        ax(ij, i) = ax(ij, i) + rrd * (tvp(ij, j) - tv(ij, j)) & ! cld
1656                * (ph(ij, j) - ph(ij, j + 1)) / p(ij, j) ! cld
1657      END DO ! cld
1658      IF (ax(ij, i)>0.0) THEN ! cld
1659        wa(ij, i) = sqrt(2. * ax(ij, i)) ! cld
1660      END IF ! cld
1661    END DO ! cld
1662    DO i = 1, nl ! cld
1663      IF (wa(ij, i)>0.0) &          ! cld
1664              siga(ij, i) = mac(ij, i) / wa(ij, i) & ! cld
1665                      * rrd * tvp(ij, i) / p(ij, i) / 100. / delta ! cld
1666      siga(ij, i) = min(siga(ij, i), 1.0) ! cld
1667      qcondc(ij, i) = siga(ij, i) * clw(ij, i) * (1. - ep(ij, i)) & ! cld
1668              + (1. - siga(ij, i)) * qcond(ij, i) ! cld
1669    END DO ! cld
1670  END DO ! cld
1671
1672END SUBROUTINE cv_yield
1673
1674SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
1675        fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
1676        qcondc1)
1677  IMPLICIT NONE
1678
1679 
1680  ! inputs:
1681  INTEGER len, ncum, nd, nloc
1682  INTEGER idcum(nloc)
1683  INTEGER iflag(nloc)
1684  REAL precip(nloc), cbmf(nloc)
1685  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1686  REAL ma(nloc, nd)
1687  REAL qcondc(nloc, nd) !cld
1688
1689  ! outputs:
1690  INTEGER iflag1(len)
1691  REAL precip1(len), cbmf1(len)
1692  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
1693  REAL ma1(len, nd)
1694  REAL qcondc1(len, nd) !cld
1695
1696  ! local variables:
1697  INTEGER i, k
1698
1699  DO i = 1, ncum
1700    precip1(idcum(i)) = precip(i)
1701    cbmf1(idcum(i)) = cbmf(i)
1702    iflag1(idcum(i)) = iflag(i)
1703  END DO
1704
1705  DO k = 1, nl
1706    DO i = 1, ncum
1707      ft1(idcum(i), k) = ft(i, k)
1708      fq1(idcum(i), k) = fq(i, k)
1709      fu1(idcum(i), k) = fu(i, k)
1710      fv1(idcum(i), k) = fv(i, k)
1711      ma1(idcum(i), k) = ma(i, k)
1712      qcondc1(idcum(i), k) = qcondc(i, k)
1713    END DO
1714  END DO
1715
1716END SUBROUTINE cv_uncompress
1717
1718
1719END MODULE lmdz_cv
Note: See TracBrowser for help on using the repository browser.