source: LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv_routines.F90 @ 5141

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

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