source: LMDZ5/trunk/libf/phylmd/cv_routines.F90 @ 2679

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

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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