source: LMDZ6/trunk/libf/phylmd/cv_routines.f90 @ 5308

Last change on this file since 5308 was 5303, checked in by abarral, 35 hours ago

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