source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/cv_routines.F90 @ 3809

Last change on this file since 3809 was 3809, checked in by ymipsl, 10 years ago

Add LMDZ in aquaplanet configuration
YM
File size: 50.2 KB
Line 
1
2! $Id: cv_routines.F90 1992 2014-03-05 13:19:12Z lguez $
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  IMPLICIT NONE
399
400  include "cvparam.h"
401
402  ! inputs:
403  INTEGER len, ncum, nd, nloc
404  INTEGER iflag1(len), nk1(len), icb1(len)
405  REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
406  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
407  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
408  REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
409  REAL tvp1(len, nd), clw1(len, nd)
410
411  ! outputs:
412  INTEGER iflag(nloc), nk(nloc), icb(nloc)
413  REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
414  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
415  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
416  REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
417  REAL tvp(nloc, nd), clw(nloc, nd)
418  REAL dph(nloc, nd)
419
420  ! local variables:
421  INTEGER i, k, nn
422  CHARACTER (LEN=20) :: modname = 'cv_compress'
423  CHARACTER (LEN=80) :: abort_message
424
425  include 'iniprint.h'
426
427
428  DO k = 1, nl + 1
429    nn = 0
430    DO i = 1, len
431      IF (iflag1(i)==0) THEN
432        nn = nn + 1
433        t(nn, k) = t1(i, k)
434        q(nn, k) = q1(i, k)
435        qs(nn, k) = qs1(i, k)
436        u(nn, k) = u1(i, k)
437        v(nn, k) = v1(i, k)
438        gz(nn, k) = gz1(i, k)
439        h(nn, k) = h1(i, k)
440        lv(nn, k) = lv1(i, k)
441        cpn(nn, k) = cpn1(i, k)
442        p(nn, k) = p1(i, k)
443        ph(nn, k) = ph1(i, k)
444        tv(nn, k) = tv1(i, k)
445        tp(nn, k) = tp1(i, k)
446        tvp(nn, k) = tvp1(i, k)
447        clw(nn, k) = clw1(i, k)
448      END IF
449    END DO
450  END DO
451
452  IF (nn/=ncum) THEN
453    WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
454    abort_message = ''
455    CALL abort_gcm(modname, abort_message, 1)
456  END IF
457
458  nn = 0
459  DO i = 1, len
460    IF (iflag1(i)==0) THEN
461      nn = nn + 1
462      cbmf(nn) = cbmf1(i)
463      plcl(nn) = plcl1(i)
464      tnk(nn) = tnk1(i)
465      qnk(nn) = qnk1(i)
466      gznk(nn) = gznk1(i)
467      nk(nn) = nk1(i)
468      icb(nn) = icb1(i)
469      iflag(nn) = iflag1(i)
470    END IF
471  END DO
472
473  DO k = 1, nl
474    DO i = 1, ncum
475      dph(i, k) = ph(i, k) - ph(i, k+1)
476    END DO
477  END DO
478
479  RETURN
480END SUBROUTINE cv_compress
481
482SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
483    gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
484  IMPLICIT NONE
485
486  ! ---------------------------------------------------------------------
487  ! Purpose:
488  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
489  ! &
490  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
491  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
492  ! &
493  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
494  ! ---------------------------------------------------------------------
495
496  include "cvthermo.h"
497  include "cvparam.h"
498
499  ! inputs:
500  INTEGER ncum, nd, nloc
501  INTEGER icb(nloc), nk(nloc)
502  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
503  REAL p(nloc, nd), dph(nloc, nd)
504  REAL tnk(nloc), qnk(nloc), gznk(nloc)
505  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
506
507  ! outputs:
508  INTEGER inb(nloc), inb1(nloc)
509  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
510  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
511  REAL frac(nloc)
512
513  ! local variables:
514  INTEGER i, k
515  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
516  REAL by, defrac
517  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
518  LOGICAL lcape(nloc)
519
520  ! =====================================================================
521  ! --- SOME INITIALIZATIONS
522  ! =====================================================================
523
524  DO k = 1, nl
525    DO i = 1, ncum
526      ep(i, k) = 0.0
527      sigp(i, k) = sigs
528    END DO
529  END DO
530
531  ! =====================================================================
532  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
533  ! =====================================================================
534
535  ! ---       The procedure is to solve the equation.
536  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
537
538  ! ***  Calculate certain parcel quantities, including static energy   ***
539
540
541  DO i = 1, ncum
542    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
543      t0)) + gznk(i)
544  END DO
545
546
547  ! ***  Find lifted parcel quantities above cloud base    ***
548
549
550  DO k = minorig + 1, nl
551    DO i = 1, ncum
552      IF (k>=(icb(i)+1)) THEN
553        tg = t(i, k)
554        qg = qs(i, k)
555        alv = lv0 - clmcpv*(t(i,k)-t0)
556
557        ! First iteration.
558
559        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
560        s = 1./s
561        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
562        tg = tg + s*(ah0(i)-ahg)
563        tg = max(tg, 35.0)
564        tc = tg - t0
565        denom = 243.5 + tc
566        IF (tc>=0.0) THEN
567          es = 6.112*exp(17.67*tc/denom)
568        ELSE
569          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
570        END IF
571        qg = eps*es/(p(i,k)-es*(1.-eps))
572
573        ! Second iteration.
574
575        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
576        s = 1./s
577        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
578        tg = tg + s*(ah0(i)-ahg)
579        tg = max(tg, 35.0)
580        tc = tg - t0
581        denom = 243.5 + tc
582        IF (tc>=0.0) THEN
583          es = 6.112*exp(17.67*tc/denom)
584        ELSE
585          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
586        END IF
587        qg = eps*es/(p(i,k)-es*(1.-eps))
588
589        alv = lv0 - clmcpv*(t(i,k)-t0)
590        ! print*,'cpd dans convect2 ',cpd
591        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
592        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
593        tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
594        ! if (.not.cpd.gt.1000.) then
595        ! print*,'CPD=',cpd
596        ! stop
597        ! endif
598        clw(i, k) = qnk(i) - qg
599        clw(i, k) = max(0.0, clw(i,k))
600        rg = qg/(1.-qnk(i))
601        tvp(i, k) = tp(i, k)*(1.+rg*epsi)
602      END IF
603    END DO
604  END DO
605
606  ! =====================================================================
607  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
608  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
609  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
610  ! =====================================================================
611
612  DO k = minorig + 1, nl
613    DO i = 1, ncum
614      IF (k>=(nk(i)+1)) THEN
615        tca = tp(i, k) - t0
616        IF (tca>=0.0) THEN
617          elacrit = elcrit
618        ELSE
619          elacrit = elcrit*(1.0-tca/tlcrit)
620        END IF
621        elacrit = max(elacrit, 0.0)
622        ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
623        ep(i, k) = max(ep(i,k), 0.0)
624        ep(i, k) = min(ep(i,k), 1.0)
625        sigp(i, k) = sigs
626      END IF
627    END DO
628  END DO
629
630  ! =====================================================================
631  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
632  ! --- VIRTUAL TEMPERATURE
633  ! =====================================================================
634
635  DO k = minorig + 1, nl
636    DO i = 1, ncum
637      IF (k>=(icb(i)+1)) THEN
638        tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
639        ! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
640        ! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
641      END IF
642    END DO
643  END DO
644  DO i = 1, ncum
645    tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd
646  END DO
647
648  ! =====================================================================
649  ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
650  ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
651  ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
652  ! =====================================================================
653
654  DO i = 1, ncum
655    cape(i) = 0.0
656    capem(i) = 0.0
657    inb(i) = icb(i) + 1
658    inb1(i) = inb(i)
659  END DO
660
661  ! Originial Code
662
663  ! do 530 k=minorig+1,nl-1
664  ! do 520 i=1,ncum
665  ! if(k.ge.(icb(i)+1))then
666  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
667  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
668  ! cape(i)=cape(i)+by
669  ! if(by.ge.0.0)inb1(i)=k+1
670  ! if(cape(i).gt.0.0)then
671  ! inb(i)=k+1
672  ! capem(i)=cape(i)
673  ! endif
674  ! endif
675  ! 520    continue
676  ! 530  continue
677  ! do 540 i=1,ncum
678  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
679  ! cape(i)=capem(i)+byp
680  ! defrac=capem(i)-cape(i)
681  ! defrac=max(defrac,0.001)
682  ! frac(i)=-cape(i)/defrac
683  ! frac(i)=min(frac(i),1.0)
684  ! frac(i)=max(frac(i),0.0)
685  ! 540   continue
686
687  ! K Emanuel fix
688
689  ! call zilch(byp,ncum)
690  ! do 530 k=minorig+1,nl-1
691  ! do 520 i=1,ncum
692  ! if(k.ge.(icb(i)+1))then
693  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
694  ! cape(i)=cape(i)+by
695  ! if(by.ge.0.0)inb1(i)=k+1
696  ! if(cape(i).gt.0.0)then
697  ! inb(i)=k+1
698  ! capem(i)=cape(i)
699  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
700  ! endif
701  ! endif
702  ! 520    continue
703  ! 530  continue
704  ! do 540 i=1,ncum
705  ! inb(i)=max(inb(i),inb1(i))
706  ! cape(i)=capem(i)+byp(i)
707  ! defrac=capem(i)-cape(i)
708  ! defrac=max(defrac,0.001)
709  ! frac(i)=-cape(i)/defrac
710  ! frac(i)=min(frac(i),1.0)
711  ! frac(i)=max(frac(i),0.0)
712  ! 540   continue
713
714  ! J Teixeira fix
715
716  CALL zilch(byp, ncum)
717  DO i = 1, ncum
718    lcape(i) = .TRUE.
719  END DO
720  DO k = minorig + 1, nl - 1
721    DO i = 1, ncum
722      IF (cape(i)<0.0) lcape(i) = .FALSE.
723      IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN
724        by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k)
725        byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1)
726        cape(i) = cape(i) + by
727        IF (by>=0.0) inb1(i) = k + 1
728        IF (cape(i)>0.0) THEN
729          inb(i) = k + 1
730          capem(i) = cape(i)
731        END IF
732      END IF
733    END DO
734  END DO
735  DO i = 1, ncum
736    cape(i) = capem(i) + byp(i)
737    defrac = capem(i) - cape(i)
738    defrac = max(defrac, 0.001)
739    frac(i) = -cape(i)/defrac
740    frac(i) = min(frac(i), 1.0)
741    frac(i) = max(frac(i), 0.0)
742  END DO
743
744  ! =====================================================================
745  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
746  ! =====================================================================
747
748  ! initialization:
749  DO i = 1, ncum*nlp
750    hp(i, 1) = h(i, 1)
751  END DO
752
753  DO k = minorig + 1, nl
754    DO i = 1, ncum
755      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
756        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
757          )
758      END IF
759    END DO
760  END DO
761
762  RETURN
763END SUBROUTINE cv_undilute2
764
765SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
766    cpn, iflag, cbmf)
767  IMPLICIT NONE
768
769  ! inputs:
770  INTEGER ncum, nd, nloc
771  INTEGER nk(nloc), icb(nloc)
772  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
773  REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent...
774  REAL plcl(nloc), cpn(nloc, nd)
775
776  ! outputs:
777  INTEGER iflag(nloc)
778  REAL cbmf(nloc) ! also an input
779
780  ! local variables:
781  INTEGER i, k, icbmax
782  REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
783  REAL work(nloc)
784
785  include "cvthermo.h"
786  include "cvparam.h"
787
788  ! -------------------------------------------------------------------
789  ! Compute icbmax.
790  ! -------------------------------------------------------------------
791
792  icbmax = 2
793  DO i = 1, ncum
794    icbmax = max(icbmax, icb(i))
795  END DO
796
797  ! =====================================================================
798  ! ---  CALCULATE CLOUD BASE MASS FLUX
799  ! =====================================================================
800
801  ! tvpplcl = parcel temperature lifted adiabatically from level
802  ! icb-1 to the LCL.
803  ! tvaplcl = virtual temperature at the LCL.
804
805  DO i = 1, ncum
806    dtpbl(i) = 0.0
807    tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( &
808      i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
809    tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
810      ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
811  END DO
812
813  ! -------------------------------------------------------------------
814  ! --- Interpolate difference between lifted parcel and
815  ! --- environmental temperatures to lifted condensation level
816  ! -------------------------------------------------------------------
817
818  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
819
820  DO k = minorig, icbmax
821    DO i = 1, ncum
822      IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
823        dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
824      END IF
825    END DO
826  END DO
827  DO i = 1, ncum
828    dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
829    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
830  END DO
831
832  ! -------------------------------------------------------------------
833  ! --- Adjust cloud base mass flux
834  ! -------------------------------------------------------------------
835
836  DO i = 1, ncum
837    work(i) = cbmf(i)
838    cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
839    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
840      iflag(i) = 3
841    END IF
842  END DO
843
844  RETURN
845END SUBROUTINE cv_closure
846
847SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
848    h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
849    sij, elij)
850  IMPLICIT NONE
851
852  include "cvthermo.h"
853  include "cvparam.h"
854
855  ! inputs:
856  INTEGER ncum, nd, nloc
857  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
858  REAL cbmf(nloc), qnk(nloc)
859  REAL ph(nloc, nd+1)
860  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
861  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
862  REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
863
864  ! outputs:
865  INTEGER nent(nloc, nd)
866  REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
867  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
868  REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
869
870  ! local variables:
871  INTEGER i, j, k, ij
872  INTEGER num1, num2
873  REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
874  REAL alt, qp1, smid, sjmin, sjmax, delp, delm
875  REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
876  REAL bsum(nloc, nd)
877  LOGICAL lwork(nloc)
878
879  ! =====================================================================
880  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
881  ! =====================================================================
882
883  DO i = 1, ncum*nlp
884    nent(i, 1) = 0
885    m(i, 1) = 0.0
886  END DO
887
888  DO k = 1, nlp
889    DO j = 1, nlp
890      DO i = 1, ncum
891        qent(i, k, j) = q(i, j)
892        uent(i, k, j) = u(i, j)
893        vent(i, k, j) = v(i, j)
894        elij(i, k, j) = 0.0
895        ment(i, k, j) = 0.0
896        sij(i, k, j) = 0.0
897      END DO
898    END DO
899  END DO
900
901  ! -------------------------------------------------------------------
902  ! --- Calculate rates of mixing,  m(i)
903  ! -------------------------------------------------------------------
904
905  CALL zilch(work, ncum)
906
907  DO j = minorig + 1, nl
908    DO i = 1, ncum
909      IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
910        k = min(j, inb1(i))
911        dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
912          entp*0.04*(ph(i,k)-ph(i,k+1))
913        work(i) = work(i) + dbo
914        m(i, j) = cbmf(i)*dbo
915      END IF
916    END DO
917  END DO
918  DO k = minorig + 1, nl
919    DO i = 1, ncum
920      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
921        m(i, k) = m(i, k)/work(i)
922      END IF
923    END DO
924  END DO
925
926
927  ! =====================================================================
928  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
929  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
930  ! --- FRACTION (sij)
931  ! =====================================================================
932
933
934  DO i = minorig + 1, nl
935    DO j = minorig + 1, nl
936      DO ij = 1, ncum
937        IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
938            inb(ij))) THEN
939          qti = qnk(ij) - ep(ij, i)*clw(ij, i)
940          bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd)
941          anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
942          denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
943          dei = denom
944          IF (abs(dei)<0.01) dei = 0.01
945          sij(ij, i, j) = anum/dei
946          sij(ij, i, i) = 1.0
947          altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
948          altem = altem/bf2
949          cwat = clw(ij, j)*(1.-ep(ij,j))
950          stemp = sij(ij, i, j)
951          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
952            anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
953            denom = denom + lv(ij, j)*(q(ij,i)-qti)
954            IF (abs(denom)<0.01) denom = 0.01
955            sij(ij, i, j) = anum/denom
956            altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
957            altem = altem - (bf2-1.)*cwat
958          END IF
959          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
960            qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
961            uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
962              (1.-sij(ij,i,j))*u(ij, nk(ij))
963            vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
964              (1.-sij(ij,i,j))*v(ij, nk(ij))
965            elij(ij, i, j) = altem
966            elij(ij, i, j) = max(0.0, elij(ij,i,j))
967            ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
968            nent(ij, i) = nent(ij, i) + 1
969          END IF
970          sij(ij, i, j) = max(0.0, sij(ij,i,j))
971          sij(ij, i, j) = min(1.0, sij(ij,i,j))
972        END IF
973      END DO
974    END DO
975
976    ! ***   If no air can entrain at level i assume that updraft detrains
977    ! ***
978    ! ***   at that level and calculate detrained air flux and properties
979    ! ***
980
981    DO ij = 1, ncum
982      IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
983        ment(ij, i, i) = m(ij, i)
984        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
985        uent(ij, i, i) = u(ij, nk(ij))
986        vent(ij, i, i) = v(ij, nk(ij))
987        elij(ij, i, i) = clw(ij, i)
988        sij(ij, i, i) = 1.0
989      END IF
990    END DO
991  END DO
992
993  DO i = 1, ncum
994    sij(i, inb(i), inb(i)) = 1.0
995  END DO
996
997  ! =====================================================================
998  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
999  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
1000  ! =====================================================================
1001
1002  CALL zilch(bsum, ncum*nlp)
1003  DO ij = 1, ncum
1004    lwork(ij) = .FALSE.
1005  END DO
1006  DO i = minorig + 1, nl
1007
1008    num1 = 0
1009    DO ij = 1, ncum
1010      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
1011    END DO
1012    IF (num1<=0) GO TO 789
1013
1014    DO ij = 1, ncum
1015      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
1016        lwork(ij) = (nent(ij,i)/=0)
1017        qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1018        anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
1019        denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
1020        IF (abs(denom)<0.01) denom = 0.01
1021        scrit(ij) = anum/denom
1022        alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
1023        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
1024        asij(ij) = 0.0
1025        smin(ij) = 1.0
1026      END IF
1027    END DO
1028    DO j = minorig, nl
1029
1030      num2 = 0
1031      DO ij = 1, ncum
1032        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1033          ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
1034      END DO
1035      IF (num2<=0) GO TO 783
1036
1037      DO ij = 1, ncum
1038        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1039            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1040          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
1041            IF (j>i) THEN
1042              smid = min(sij(ij,i,j), scrit(ij))
1043              sjmax = smid
1044              sjmin = smid
1045              IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
1046                smin(ij) = smid
1047                sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
1048                sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
1049                sjmin = min(sjmin, scrit(ij))
1050              END IF
1051            ELSE
1052              sjmax = max(sij(ij,i,j+1), scrit(ij))
1053              smid = max(sij(ij,i,j), scrit(ij))
1054              sjmin = 0.0
1055              IF (j>1) sjmin = sij(ij, i, j-1)
1056              sjmin = max(sjmin, scrit(ij))
1057            END IF
1058            delp = abs(sjmax-smid)
1059            delm = abs(sjmin-smid)
1060            asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
1061            ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
1062          END IF
1063        END IF
1064      END DO
1065783 END DO
1066    DO ij = 1, ncum
1067      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
1068        asij(ij) = max(1.0E-21, asij(ij))
1069        asij(ij) = 1.0/asij(ij)
1070        bsum(ij, i) = 0.0
1071      END IF
1072    END DO
1073    DO j = minorig, nl + 1
1074      DO ij = 1, ncum
1075        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1076            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1077          ment(ij, i, j) = ment(ij, i, j)*asij(ij)
1078          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
1079        END IF
1080      END DO
1081    END DO
1082    DO ij = 1, ncum
1083      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
1084          i)<1.0E-18) .AND. lwork(ij)) THEN
1085        nent(ij, i) = 0
1086        ment(ij, i, i) = m(ij, i)
1087        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1088        uent(ij, i, i) = u(ij, nk(ij))
1089        vent(ij, i, i) = v(ij, nk(ij))
1090        elij(ij, i, i) = clw(ij, i)
1091        sij(ij, i, i) = 1.0
1092      END IF
1093    END DO
1094789 END DO
1095
1096  RETURN
1097END SUBROUTINE cv_mixing
1098
1099SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
1100    ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
1101  IMPLICIT NONE
1102
1103
1104  include "cvthermo.h"
1105  include "cvparam.h"
1106
1107  ! inputs:
1108  INTEGER ncum, nd, nloc
1109  INTEGER inb(nloc)
1110  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
1111  REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
1112  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1113  REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
1114  REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
1115
1116  ! outputs:
1117  INTEGER iflag(nloc) ! also an input
1118  REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd)
1119  REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd)
1120
1121  ! local variables:
1122  INTEGER i, j, k, ij, num1
1123  INTEGER jtt(nloc)
1124  REAL awat, coeff, qsm, afac, sigt, b6, c6, revap
1125  REAL dhdp, fac, qstm, rat
1126  REAL wdtrain(nloc)
1127  LOGICAL lwork(nloc)
1128
1129  ! =====================================================================
1130  ! --- PRECIPITATING DOWNDRAFT CALCULATION
1131  ! =====================================================================
1132
1133  ! Initializations:
1134
1135  DO i = 1, ncum
1136    DO k = 1, nl + 1
1137      wt(i, k) = omtsnow
1138      mp(i, k) = 0.0
1139      evap(i, k) = 0.0
1140      water(i, k) = 0.0
1141    END DO
1142  END DO
1143
1144  DO i = 1, ncum
1145    qp(i, 1) = q(i, 1)
1146    up(i, 1) = u(i, 1)
1147    vp(i, 1) = v(i, 1)
1148  END DO
1149
1150  DO k = 2, nl + 1
1151    DO i = 1, ncum
1152      qp(i, k) = q(i, k-1)
1153      up(i, k) = u(i, k-1)
1154      vp(i, k) = v(i, k-1)
1155    END DO
1156  END DO
1157
1158
1159  ! ***  Check whether ep(inb)=0, if so, skip precipitating    ***
1160  ! ***             downdraft calculation                      ***
1161
1162
1163  ! ***  Integrate liquid water equation to find condensed water   ***
1164  ! ***                and condensed water flux                    ***
1165
1166
1167  DO i = 1, ncum
1168    jtt(i) = 2
1169    IF (ep(i,inb(i))<=0.0001) iflag(i) = 2
1170    IF (iflag(i)==0) THEN
1171      lwork(i) = .TRUE.
1172    ELSE
1173      lwork(i) = .FALSE.
1174    END IF
1175  END DO
1176
1177  ! ***                    Begin downdraft loop                    ***
1178
1179
1180  CALL zilch(wdtrain, ncum)
1181  DO i = nl + 1, 1, -1
1182
1183    num1 = 0
1184    DO ij = 1, ncum
1185      IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
1186    END DO
1187    IF (num1<=0) GO TO 899
1188
1189
1190    ! ***        Calculate detrained precipitation             ***
1191
1192    DO ij = 1, ncum
1193      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1194        wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i)
1195      END IF
1196    END DO
1197
1198    IF (i>1) THEN
1199      DO j = 1, i - 1
1200        DO ij = 1, ncum
1201          IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1202            awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i)
1203            awat = max(0.0, awat)
1204            wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i)
1205          END IF
1206        END DO
1207      END DO
1208    END IF
1209
1210    ! ***    Find rain water and evaporation using provisional   ***
1211    ! ***              estimates of qp(i)and qp(i-1)             ***
1212
1213
1214    ! ***  Value of terminal velocity and coeffecient of evaporation for snow
1215    ! ***
1216
1217    DO ij = 1, ncum
1218      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1219        coeff = coeffs
1220        wt(ij, i) = omtsnow
1221
1222        ! ***  Value of terminal velocity and coeffecient of evaporation for
1223        ! rain   ***
1224
1225        IF (t(ij,i)>273.0) THEN
1226          coeff = coeffr
1227          wt(ij, i) = omtrain
1228        END IF
1229        qsm = 0.5*(q(ij,i)+qp(ij,i+1))
1230        afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i))
1231        afac = max(afac, 0.0)
1232        sigt = sigp(ij, i)
1233        sigt = max(0.0, sigt)
1234        sigt = min(1.0, sigt)
1235        b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i)
1236        c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i)
1237        revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
1238        evap(ij, i) = sigt*afac*revap
1239        water(ij, i) = revap*revap
1240
1241        ! ***  Calculate precipitating downdraft mass flux under     ***
1242        ! ***              hydrostatic approximation                 ***
1243
1244        IF (i>1) THEN
1245          dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
1246          dhdp = max(dhdp, 10.0)
1247          mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp
1248          mp(ij, i) = max(mp(ij,i), 0.0)
1249
1250          ! ***   Add small amount of inertia to downdraft              ***
1251
1252          fac = 20.0/(ph(ij,i-1)-ph(ij,i))
1253          mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
1254
1255          ! ***      Force mp to decrease linearly to zero
1256          ! ***
1257          ! ***      between about 950 mb and the surface
1258          ! ***
1259
1260          IF (p(ij,i)>(0.949*p(ij,1))) THEN
1261            jtt(ij) = max(jtt(ij), i)
1262            mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ &
1263              (p(ij,1)-p(ij,jtt(ij)))
1264          END IF
1265        END IF
1266
1267        ! ***       Find mixing ratio of precipitating downdraft     ***
1268
1269        IF (i/=inb(ij)) THEN
1270          IF (i==1) THEN
1271            qstm = qs(ij, 1)
1272          ELSE
1273            qstm = qs(ij, i-1)
1274          END IF
1275          IF (mp(ij,i)>mp(ij,i+1)) THEN
1276            rat = mp(ij, i+1)/mp(ij, i)
1277            qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + &
1278              100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
1279            up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat)
1280            vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat)
1281          ELSE
1282            IF (mp(ij,i+1)>0.0) THEN
1283              qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, &
1284                i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, &
1285                i)))/(lv(ij,i)+t(ij,i)*(cl-cpd))
1286              up(ij, i) = up(ij, i+1)
1287              vp(ij, i) = vp(ij, i+1)
1288            END IF
1289          END IF
1290          qp(ij, i) = min(qp(ij,i), qstm)
1291          qp(ij, i) = max(qp(ij,i), 0.0)
1292        END IF
1293      END IF
1294    END DO
1295899 END DO
1296
1297  RETURN
1298END SUBROUTINE cv_unsat
1299
1300SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
1301    ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
1302    ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
1303    precip, cbmf, ft, fq, fu, fv, ma, qcondc)
1304  IMPLICIT NONE
1305
1306  include "cvthermo.h"
1307  include "cvparam.h"
1308
1309  ! inputs
1310  INTEGER ncum, nd, nloc
1311  INTEGER nk(nloc), icb(nloc), inb(nloc)
1312  INTEGER nent(nloc, nd)
1313  REAL delt
1314  REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
1315  REAL gz(nloc, nd)
1316  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1317  REAL hp(nloc, nd), lv(nloc, nd)
1318  REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
1319  REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd)
1320  REAL up(nloc, nd), vp(nloc, nd)
1321  REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd)
1322  REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd)
1323  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
1324  REAL tv(nloc, nd), tvp(nloc, nd)
1325
1326  ! outputs
1327  INTEGER iflag(nloc) ! also an input
1328  REAL cbmf(nloc) ! also an input
1329  REAL wd(nloc), tprime(nloc), qprime(nloc)
1330  REAL precip(nloc)
1331  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1332  REAL ma(nloc, nd)
1333  REAL qcondc(nloc, nd)
1334
1335  ! local variables
1336  INTEGER i, j, ij, k, num1
1337  REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti
1338  REAL work(nloc), am(nloc), amp1(nloc), ad(nloc)
1339  REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd)
1340  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
1341  REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld
1342
1343
1344  ! -- initializations:
1345
1346  delti = 1.0/delt
1347
1348  DO i = 1, ncum
1349    precip(i) = 0.0
1350    wd(i) = 0.0
1351    tprime(i) = 0.0
1352    qprime(i) = 0.0
1353    DO k = 1, nl + 1
1354      ft(i, k) = 0.0
1355      fu(i, k) = 0.0
1356      fv(i, k) = 0.0
1357      fq(i, k) = 0.0
1358      lvcp(i, k) = lv(i, k)/cpn(i, k)
1359      qcondc(i, k) = 0.0 ! cld
1360      qcond(i, k) = 0.0 ! cld
1361      nqcond(i, k) = 0.0 ! cld
1362    END DO
1363  END DO
1364
1365
1366  ! ***  Calculate surface precipitation in mm/day     ***
1367
1368  DO i = 1, ncum
1369    IF (iflag(i)<=1) THEN
1370      ! c            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1371      ! c     &                /(rowl*g)
1372      ! c            precip(i)=precip(i)*delt/86400.
1373      precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g
1374    END IF
1375  END DO
1376
1377
1378  ! ***  Calculate downdraft velocity scale and surface temperature and  ***
1379  ! ***                    water vapor fluctuations                      ***
1380
1381  DO i = 1, ncum
1382    wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i)))
1383    qprime(i) = 0.5*(qp(i,1)-q(i,1))
1384    tprime(i) = lv0*qprime(i)/cpd
1385  END DO
1386
1387  ! ***  Calculate tendencies of lowest level potential temperature  ***
1388  ! ***                      and mixing ratio                        ***
1389
1390  DO i = 1, ncum
1391    work(i) = 0.01/(ph(i,1)-ph(i,2))
1392    am(i) = 0.0
1393  END DO
1394  DO k = 2, nl
1395    DO i = 1, ncum
1396      IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
1397        am(i) = am(i) + m(i, k)
1398      END IF
1399    END DO
1400  END DO
1401  DO i = 1, ncum
1402    IF ((g*work(i)*am(i))>=delti) iflag(i) = 1
1403    ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, &
1404      1))/cpn(i,1))
1405    ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1)
1406    ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* &
1407      work(i)/cpn(i, 1)
1408    fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + &
1409      sigd*evap(i, 1)
1410    fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i)
1411    fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, &
1412      2)-u(i,1)))
1413    fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, &
1414      2)-v(i,1)))
1415  END DO
1416  DO j = 2, nl
1417    DO i = 1, ncum
1418      IF (j<=inb(i)) THEN
1419        fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1))
1420        fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1))
1421        fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1))
1422      END IF
1423    END DO
1424  END DO
1425
1426  ! ***  Calculate tendencies of potential temperature and mixing ratio  ***
1427  ! ***               at levels above the lowest level                   ***
1428
1429  ! ***  First find the net saturated updraft and downdraft mass fluxes  ***
1430  ! ***                      through each level                          ***
1431
1432  DO i = 2, nl + 1
1433
1434    num1 = 0
1435    DO ij = 1, ncum
1436      IF (i<=inb(ij)) num1 = num1 + 1
1437    END DO
1438    IF (num1<=0) GO TO 1500
1439
1440    CALL zilch(amp1, ncum)
1441    CALL zilch(ad, ncum)
1442
1443    DO k = i + 1, nl + 1
1444      DO ij = 1, ncum
1445        IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN
1446          amp1(ij) = amp1(ij) + m(ij, k)
1447        END IF
1448      END DO
1449    END DO
1450
1451    DO k = 1, i
1452      DO j = i + 1, nl + 1
1453        DO ij = 1, ncum
1454          IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN
1455            amp1(ij) = amp1(ij) + ment(ij, k, j)
1456          END IF
1457        END DO
1458      END DO
1459    END DO
1460    DO k = 1, i - 1
1461      DO j = i, nl + 1
1462        DO ij = 1, ncum
1463          IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
1464            ad(ij) = ad(ij) + ment(ij, j, k)
1465          END IF
1466        END DO
1467      END DO
1468    END DO
1469
1470    DO ij = 1, ncum
1471      IF (i<=inb(ij)) THEN
1472        dpinv = 0.01/(ph(ij,i)-ph(ij,i+1))
1473        cpinv = 1.0/cpn(ij, i)
1474
1475        ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, &
1476          i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, &
1477          i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i)
1478        ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij &
1479          ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
1480        ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( &
1481          ij,i+1)-t(ij,i))*dpinv*cpinv
1482        fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, &
1483          i))-ad(ij)*(q(ij,i)-q(ij,i-1)))
1484        fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, &
1485          i))-ad(ij)*(u(ij,i)-u(ij,i-1)))
1486        fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, &
1487          i))-ad(ij)*(v(ij,i)-v(ij,i-1)))
1488      END IF
1489    END DO
1490    DO k = 1, i - 1
1491      DO ij = 1, ncum
1492        IF (i<=inb(ij)) THEN
1493          awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i)
1494          awat = max(awat, 0.0)
1495          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q &
1496            (ij,i))
1497          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1498            ))
1499          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1500            ))
1501          ! (saturated updrafts resulting from mixing)               ! cld
1502          qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld
1503          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1504        END IF
1505      END DO
1506    END DO
1507    DO k = i, nl + 1
1508      DO ij = 1, ncum
1509        IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
1510          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i &
1511            ))
1512          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1513            ))
1514          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1515            ))
1516        END IF
1517      END DO
1518    END DO
1519    DO ij = 1, ncum
1520      IF (i<=inb(ij)) THEN
1521        fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, &
1522          i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
1523        fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, &
1524          i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv
1525        fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, &
1526          i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv
1527        ! (saturated downdrafts resulting from mixing)               ! cld
1528        DO k = i + 1, inb(ij) ! cld
1529          qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld
1530          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1531        END DO ! cld
1532        ! (particular case: no detraining level is found)            ! cld
1533        IF (nent(ij,i)==0) THEN ! cld
1534          qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld
1535          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1536        END IF ! cld
1537        IF (nqcond(ij,i)/=0.) THEN ! cld
1538          qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld
1539        END IF ! cld
1540      END IF
1541    END DO
15421500 END DO
1543
1544  ! *** Adjust tendencies at top of convection layer to reflect  ***
1545  ! ***       actual position of the level zero cape             ***
1546
1547  DO ij = 1, ncum
1548    fqold = fq(ij, inb(ij))
1549    fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij))
1550    fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, &
1551      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1552      inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1)
1553    ftold = ft(ij, inb(ij))
1554    ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij))
1555    ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, &
1556      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1557      inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1)
1558    fuold = fu(ij, inb(ij))
1559    fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij))
1560    fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, &
1561      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1562    fvold = fv(ij, inb(ij))
1563    fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij))
1564    fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, &
1565      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1566  END DO
1567
1568  ! ***   Very slightly adjust tendencies to force exact   ***
1569  ! ***     enthalpy, momentum and tracer conservation     ***
1570
1571  DO ij = 1, ncum
1572    ents(ij) = 0.0
1573    uav(ij) = 0.0
1574    vav(ij) = 0.0
1575    DO i = 1, inb(ij)
1576      ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- &
1577        ph(ij,i+1))
1578      uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1))
1579      vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1))
1580    END DO
1581  END DO
1582  DO ij = 1, ncum
1583    ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1584    uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1585    vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1586  END DO
1587  DO ij = 1, ncum
1588    DO i = 1, inb(ij)
1589      ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i)
1590      fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij))
1591      fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij))
1592    END DO
1593  END DO
1594
1595  DO k = 1, nl + 1
1596    DO i = 1, ncum
1597      IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10
1598    END DO
1599  END DO
1600
1601
1602  DO i = 1, ncum
1603    IF (iflag(i)>2) THEN
1604      precip(i) = 0.0
1605      cbmf(i) = 0.0
1606    END IF
1607  END DO
1608  DO k = 1, nl
1609    DO i = 1, ncum
1610      IF (iflag(i)>2) THEN
1611        ft(i, k) = 0.0
1612        fq(i, k) = 0.0
1613        fu(i, k) = 0.0
1614        fv(i, k) = 0.0
1615        qcondc(i, k) = 0.0 ! cld
1616      END IF
1617    END DO
1618  END DO
1619
1620  DO k = 1, nl + 1
1621    DO i = 1, ncum
1622      ma(i, k) = 0.
1623    END DO
1624  END DO
1625  DO k = nl, 1, -1
1626    DO i = 1, ncum
1627      ma(i, k) = ma(i, k+1) + m(i, k)
1628    END DO
1629  END DO
1630
1631
1632  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
1633  ! ***           of condensed water         ***            ! cld
1634  ! ! cld
1635  DO ij = 1, ncum ! cld
1636    DO i = 1, nd ! cld
1637      mac(ij, i) = 0.0 ! cld
1638      wa(ij, i) = 0.0 ! cld
1639      siga(ij, i) = 0.0 ! cld
1640    END DO ! cld
1641    DO i = nk(ij), inb(ij) ! cld
1642      DO k = i + 1, inb(ij) + 1 ! cld
1643        mac(ij, i) = mac(ij, i) + m(ij, k) ! cld
1644      END DO ! cld
1645    END DO ! cld
1646    DO i = icb(ij), inb(ij) - 1 ! cld
1647      ax(ij, i) = 0. ! cld
1648      DO j = icb(ij), i ! cld
1649        ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld
1650          *(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld
1651      END DO ! cld
1652      IF (ax(ij,i)>0.0) THEN ! cld
1653        wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld
1654      END IF ! cld
1655    END DO ! cld
1656    DO i = 1, nl ! cld
1657      IF (wa(ij,i)>0.0) &          ! cld
1658        siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld
1659        *rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld
1660      siga(ij, i) = min(siga(ij,i), 1.0) ! cld
1661      qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld
1662        +(1.-siga(ij,i))*qcond(ij, i) ! cld
1663    END DO ! cld
1664  END DO ! cld
1665
1666  RETURN
1667END SUBROUTINE cv_yield
1668
1669SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
1670    fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
1671    qcondc1)
1672  IMPLICIT NONE
1673
1674  include "cvparam.h"
1675
1676  ! inputs:
1677  INTEGER len, ncum, nd, nloc
1678  INTEGER idcum(nloc)
1679  INTEGER iflag(nloc)
1680  REAL precip(nloc), cbmf(nloc)
1681  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1682  REAL ma(nloc, nd)
1683  REAL qcondc(nloc, nd) !cld
1684
1685  ! outputs:
1686  INTEGER iflag1(len)
1687  REAL precip1(len), cbmf1(len)
1688  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
1689  REAL ma1(len, nd)
1690  REAL qcondc1(len, nd) !cld
1691
1692  ! local variables:
1693  INTEGER i, k
1694
1695  DO i = 1, ncum
1696    precip1(idcum(i)) = precip(i)
1697    cbmf1(idcum(i)) = cbmf(i)
1698    iflag1(idcum(i)) = iflag(i)
1699  END DO
1700
1701  DO k = 1, nl
1702    DO i = 1, ncum
1703      ft1(idcum(i), k) = ft(i, k)
1704      fq1(idcum(i), k) = fq(i, k)
1705      fu1(idcum(i), k) = fu(i, k)
1706      fv1(idcum(i), k) = fv(i, k)
1707      ma1(idcum(i), k) = ma(i, k)
1708      qcondc1(idcum(i), k) = qcondc(i, k)
1709    END DO
1710  END DO
1711
1712  RETURN
1713END SUBROUTINE cv_uncompress
1714
Note: See TracBrowser for help on using the repository browser.