source: LMDZ6/trunk/libf/phylmdiso/cv_routines.f90 @ 5452

Last change on this file since 5452 was 5303, checked in by abarral, 2 months ago

Turn compar1d.h date_cas.h into module
Move fcg_racmo.h to obsolete
  • Property svn:keywords set to Id
File size: 50.1 KB
Line 
1
2! $Id: cv_routines.f90 5303 2024-10-30 17:34:05Z aborella $
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  ! inputs:
260  INTEGER len, nd
261  INTEGER nk(len), icb(len), icbmax
262  REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
263  REAL p(len, nd)
264
265  ! outputs:
266  REAL tp(len, nd), tvp(len, nd), clw(len, nd)
267
268  ! local variables:
269  INTEGER i, k
270  REAL tg, qg, alv, s, ahg, tc, denom, es, rg
271  REAL ah0(len), cpp(len)
272  REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
273
274  ! -------------------------------------------------------------------
275  ! --- Calculates the lifted parcel virtual temperature at nk,
276  ! --- the actual temperature, and the adiabatic
277  ! --- liquid water content. The procedure is to solve the equation.
278  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
279  ! -------------------------------------------------------------------
280
281  DO i = 1, len
282    tnk(i) = t(i, nk(i))
283    qnk(i) = q(i, nk(i))
284    gznk(i) = gz(i, nk(i))
285    ticb(i) = t(i, icb(i))
286    gzicb(i) = gz(i, icb(i))
287  END DO
288
289  ! ***  Calculate certain parcel quantities, including static energy   ***
290
291  DO i = 1, len
292    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
293      273.15)) + gznk(i)
294    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
295  END DO
296
297  ! ***   Calculate lifted parcel quantities below cloud base   ***
298
299  DO k = minorig, icbmax - 1
300    DO i = 1, len
301      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
302      tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
303    END DO
304  END DO
305
306  ! ***  Find lifted parcel quantities above cloud base    ***
307
308  DO i = 1, len
309    tg = ticb(i)
310    qg = qs(i, icb(i))
311    alv = lv0 - clmcpv*(ticb(i)-t0)
312
313    ! First iteration.
314
315    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
316    s = 1./s
317    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
318    tg = tg + s*(ah0(i)-ahg)
319    tg = max(tg, 35.0)
320    tc = tg - t0
321    denom = 243.5 + tc
322    IF (tc>=0.0) THEN
323      es = 6.112*exp(17.67*tc/denom)
324    ELSE
325      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
326    END IF
327    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
328
329    ! Second iteration.
330
331    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
332    s = 1./s
333    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
334    tg = tg + s*(ah0(i)-ahg)
335    tg = max(tg, 35.0)
336    tc = tg - t0
337    denom = 243.5 + tc
338    IF (tc>=0.0) THEN
339      es = 6.112*exp(17.67*tc/denom)
340    ELSE
341      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
342    END IF
343    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
344
345    alv = lv0 - clmcpv*(ticb(i)-273.15)
346    tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
347    clw(i, icb(i)) = qnk(i) - qg
348    clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
349    rg = qg/(1.-qnk(i))
350    tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
351  END DO
352
353  DO k = minorig, icbmax
354    DO i = 1, len
355      tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
356    END DO
357  END DO
358
359  RETURN
360END SUBROUTINE cv_undilute1
361
362SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
363USE cvparam_mod_h
364    IMPLICIT NONE
365
366  ! -------------------------------------------------------------------
367  ! --- Test for instability.
368  ! --- If there was no convection at last time step and parcel
369  ! --- is stable at icb, then set iflag to 4.
370  ! -------------------------------------------------------------------
371
372
373  ! inputs:
374  INTEGER len, nd, icb(len)
375  REAL cbmf(len), tv(len, nd), tvp(len, nd)
376
377  ! outputs:
378  INTEGER iflag(len) ! also an input
379
380  ! local variables:
381  INTEGER i
382
383
384  DO i = 1, len
385    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
386      icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
387  END DO
388
389  RETURN
390END SUBROUTINE cv_trigger
391
392SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
393    tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
394    tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
395    v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
396USE cvparam_mod_h
397    USE print_control_mod, ONLY: lunout
398  IMPLICIT NONE
399
400
401  ! inputs:
402  INTEGER len, ncum, nd, nloc
403  INTEGER iflag1(len), nk1(len), icb1(len)
404  REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
405  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
406  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
407  REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
408  REAL tvp1(len, nd), clw1(len, nd)
409
410  ! outputs:
411  INTEGER iflag(nloc), nk(nloc), icb(nloc)
412  REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
413  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
414  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
415  REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
416  REAL tvp(nloc, nd), clw(nloc, nd)
417  REAL dph(nloc, nd)
418
419  ! local variables:
420  INTEGER i, k, nn
421  CHARACTER (LEN=20) :: modname = 'cv_compress'
422  CHARACTER (LEN=80) :: abort_message
423
424
425  DO k = 1, nl + 1
426    nn = 0
427    DO i = 1, len
428      IF (iflag1(i)==0) THEN
429        nn = nn + 1
430        t(nn, k) = t1(i, k)
431        q(nn, k) = q1(i, k)
432        qs(nn, k) = qs1(i, k)
433        u(nn, k) = u1(i, k)
434        v(nn, k) = v1(i, k)
435        gz(nn, k) = gz1(i, k)
436        h(nn, k) = h1(i, k)
437        lv(nn, k) = lv1(i, k)
438        cpn(nn, k) = cpn1(i, k)
439        p(nn, k) = p1(i, k)
440        ph(nn, k) = ph1(i, k)
441        tv(nn, k) = tv1(i, k)
442        tp(nn, k) = tp1(i, k)
443        tvp(nn, k) = tvp1(i, k)
444        clw(nn, k) = clw1(i, k)
445      END IF
446    END DO
447  END DO
448
449  IF (nn/=ncum) THEN
450    WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
451    abort_message = ''
452    CALL abort_physic(modname, abort_message, 1)
453  END IF
454
455  nn = 0
456  DO i = 1, len
457    IF (iflag1(i)==0) THEN
458      nn = nn + 1
459      cbmf(nn) = cbmf1(i)
460      plcl(nn) = plcl1(i)
461      tnk(nn) = tnk1(i)
462      qnk(nn) = qnk1(i)
463      gznk(nn) = gznk1(i)
464      nk(nn) = nk1(i)
465      icb(nn) = icb1(i)
466      iflag(nn) = iflag1(i)
467    END IF
468  END DO
469
470  DO k = 1, nl
471    DO i = 1, ncum
472      dph(i, k) = ph(i, k) - ph(i, k+1)
473    END DO
474  END DO
475
476  RETURN
477END SUBROUTINE cv_compress
478
479SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
480    gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
481USE cvparam_mod_h
482    USE cvthermo_mod_h
483  IMPLICIT NONE
484
485  ! ---------------------------------------------------------------------
486  ! Purpose:
487  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
488  ! &
489  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
490  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
491  ! &
492  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
493  ! ---------------------------------------------------------------------
494
495  ! inputs:
496  INTEGER ncum, nd, nloc
497  INTEGER icb(nloc), nk(nloc)
498  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
499  REAL p(nloc, nd), dph(nloc, nd)
500  REAL tnk(nloc), qnk(nloc), gznk(nloc)
501  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
502
503  ! outputs:
504  INTEGER inb(nloc), inb1(nloc)
505  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
506  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
507  REAL frac(nloc)
508
509  ! local variables:
510  INTEGER i, k
511  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
512  REAL by, defrac
513  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
514  LOGICAL lcape(nloc)
515
516  ! =====================================================================
517  ! --- SOME INITIALIZATIONS
518  ! =====================================================================
519
520  DO k = 1, nl
521    DO i = 1, ncum
522      ep(i, k) = 0.0
523      sigp(i, k) = sigs
524    END DO
525  END DO
526
527  ! =====================================================================
528  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
529  ! =====================================================================
530
531  ! ---       The procedure is to solve the equation.
532  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
533
534  ! ***  Calculate certain parcel quantities, including static energy   ***
535
536
537  DO i = 1, ncum
538    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
539      t0)) + gznk(i)
540  END DO
541
542
543  ! ***  Find lifted parcel quantities above cloud base    ***
544
545
546  DO k = minorig + 1, nl
547    DO i = 1, ncum
548      IF (k>=(icb(i)+1)) THEN
549        tg = t(i, k)
550        qg = qs(i, k)
551        alv = lv0 - clmcpv*(t(i,k)-t0)
552
553        ! First iteration.
554
555        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
556        s = 1./s
557        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
558        tg = tg + s*(ah0(i)-ahg)
559        tg = max(tg, 35.0)
560        tc = tg - t0
561        denom = 243.5 + tc
562        IF (tc>=0.0) THEN
563          es = 6.112*exp(17.67*tc/denom)
564        ELSE
565          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
566        END IF
567        qg = eps*es/(p(i,k)-es*(1.-eps))
568
569        ! Second iteration.
570
571        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
572        s = 1./s
573        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
574        tg = tg + s*(ah0(i)-ahg)
575        tg = max(tg, 35.0)
576        tc = tg - t0
577        denom = 243.5 + tc
578        IF (tc>=0.0) THEN
579          es = 6.112*exp(17.67*tc/denom)
580        ELSE
581          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
582        END IF
583        qg = eps*es/(p(i,k)-es*(1.-eps))
584
585        alv = lv0 - clmcpv*(t(i,k)-t0)
586        ! print*,'cpd dans convect2 ',cpd
587        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
588        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
589        tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
590        ! if (.not.cpd.gt.1000.) then
591        ! print*,'CPD=',cpd
592        ! stop
593        ! endif
594        clw(i, k) = qnk(i) - qg
595        clw(i, k) = max(0.0, clw(i,k))
596        rg = qg/(1.-qnk(i))
597        tvp(i, k) = tp(i, k)*(1.+rg*epsi)
598      END IF
599    END DO
600  END DO
601
602  ! =====================================================================
603  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
604  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
605  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
606  ! =====================================================================
607
608  DO k = minorig + 1, nl
609    DO i = 1, ncum
610      IF (k>=(nk(i)+1)) THEN
611        tca = tp(i, k) - t0
612        IF (tca>=0.0) THEN
613          elacrit = elcrit
614        ELSE
615          elacrit = elcrit*(1.0-tca/tlcrit)
616        END IF
617        elacrit = max(elacrit, 0.0)
618        ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
619        ep(i, k) = max(ep(i,k), 0.0)
620        ep(i, k) = min(ep(i,k), 1.0)
621        sigp(i, k) = sigs
622      END IF
623    END DO
624  END DO
625
626  ! =====================================================================
627  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
628  ! --- VIRTUAL TEMPERATURE
629  ! =====================================================================
630
631  DO k = minorig + 1, nl
632    DO i = 1, ncum
633      IF (k>=(icb(i)+1)) THEN
634        tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
635        ! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
636        ! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
637      END IF
638    END DO
639  END DO
640  DO i = 1, ncum
641    tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd
642  END DO
643
644  ! =====================================================================
645  ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
646  ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
647  ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
648  ! =====================================================================
649
650  DO i = 1, ncum
651    cape(i) = 0.0
652    capem(i) = 0.0
653    inb(i) = icb(i) + 1
654    inb1(i) = inb(i)
655  END DO
656
657  ! Originial Code
658
659  ! do 530 k=minorig+1,nl-1
660  ! do 520 i=1,ncum
661  ! if(k.ge.(icb(i)+1))then
662  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
663  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
664  ! cape(i)=cape(i)+by
665  ! if(by.ge.0.0)inb1(i)=k+1
666  ! if(cape(i).gt.0.0)then
667  ! inb(i)=k+1
668  ! capem(i)=cape(i)
669  ! endif
670  ! endif
671  ! 520    continue
672  ! 530  continue
673  ! do 540 i=1,ncum
674  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
675  ! cape(i)=capem(i)+byp
676  ! defrac=capem(i)-cape(i)
677  ! defrac=max(defrac,0.001)
678  ! frac(i)=-cape(i)/defrac
679  ! frac(i)=min(frac(i),1.0)
680  ! frac(i)=max(frac(i),0.0)
681  ! 540   continue
682
683  ! K Emanuel fix
684
685  ! call zilch(byp,ncum)
686  ! do 530 k=minorig+1,nl-1
687  ! do 520 i=1,ncum
688  ! if(k.ge.(icb(i)+1))then
689  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
690  ! cape(i)=cape(i)+by
691  ! if(by.ge.0.0)inb1(i)=k+1
692  ! if(cape(i).gt.0.0)then
693  ! inb(i)=k+1
694  ! capem(i)=cape(i)
695  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
696  ! endif
697  ! endif
698  ! 520    continue
699  ! 530  continue
700  ! do 540 i=1,ncum
701  ! inb(i)=max(inb(i),inb1(i))
702  ! cape(i)=capem(i)+byp(i)
703  ! defrac=capem(i)-cape(i)
704  ! defrac=max(defrac,0.001)
705  ! frac(i)=-cape(i)/defrac
706  ! frac(i)=min(frac(i),1.0)
707  ! frac(i)=max(frac(i),0.0)
708  ! 540   continue
709
710  ! J Teixeira fix
711
712  CALL zilch(byp, ncum)
713  DO i = 1, ncum
714    lcape(i) = .TRUE.
715  END DO
716  DO k = minorig + 1, nl - 1
717    DO i = 1, ncum
718      IF (cape(i)<0.0) lcape(i) = .FALSE.
719      IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN
720        by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k)
721        byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1)
722        cape(i) = cape(i) + by
723        IF (by>=0.0) inb1(i) = k + 1
724        IF (cape(i)>0.0) THEN
725          inb(i) = k + 1
726          capem(i) = cape(i)
727        END IF
728      END IF
729    END DO
730  END DO
731  DO i = 1, ncum
732    cape(i) = capem(i) + byp(i)
733    defrac = capem(i) - cape(i)
734    defrac = max(defrac, 0.001)
735    frac(i) = -cape(i)/defrac
736    frac(i) = min(frac(i), 1.0)
737    frac(i) = max(frac(i), 0.0)
738  END DO
739
740  ! =====================================================================
741  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
742  ! =====================================================================
743
744  ! initialization:
745  DO i = 1, ncum*nlp
746    hp(i, 1) = h(i, 1)
747  END DO
748
749  DO k = minorig + 1, nl
750    DO i = 1, ncum
751      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
752        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
753          )
754      END IF
755    END DO
756  END DO
757
758  RETURN
759END SUBROUTINE cv_undilute2
760
761SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
762    cpn, iflag, cbmf)
763  USE cvthermo_mod_h
764  USE cvparam_mod_h
765  IMPLICIT NONE
766
767  ! inputs:
768  INTEGER ncum, nd, nloc
769  INTEGER nk(nloc), icb(nloc)
770  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
771  REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent...
772  REAL plcl(nloc), cpn(nloc, nd)
773
774  ! outputs:
775  INTEGER iflag(nloc)
776  REAL cbmf(nloc) ! also an input
777
778  ! local variables:
779  INTEGER i, k, icbmax
780  REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
781  REAL work(nloc)
782
783  ! -------------------------------------------------------------------
784  ! Compute icbmax.
785  ! -------------------------------------------------------------------
786
787  icbmax = 2
788  DO i = 1, ncum
789    icbmax = max(icbmax, icb(i))
790  END DO
791
792  ! =====================================================================
793  ! ---  CALCULATE CLOUD BASE MASS FLUX
794  ! =====================================================================
795
796  ! tvpplcl = parcel temperature lifted adiabatically from level
797  ! icb-1 to the LCL.
798  ! tvaplcl = virtual temperature at the LCL.
799
800  DO i = 1, ncum
801    dtpbl(i) = 0.0
802    tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( &
803      i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
804    tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
805      ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
806  END DO
807
808  ! -------------------------------------------------------------------
809  ! --- Interpolate difference between lifted parcel and
810  ! --- environmental temperatures to lifted condensation level
811  ! -------------------------------------------------------------------
812
813  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
814
815  DO k = minorig, icbmax
816    DO i = 1, ncum
817      IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
818        dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
819      END IF
820    END DO
821  END DO
822  DO i = 1, ncum
823    dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
824    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
825  END DO
826
827  ! -------------------------------------------------------------------
828  ! --- Adjust cloud base mass flux
829  ! -------------------------------------------------------------------
830
831  DO i = 1, ncum
832    work(i) = cbmf(i)
833    cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
834    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
835      iflag(i) = 3
836    END IF
837  END DO
838
839  RETURN
840END SUBROUTINE cv_closure
841
842SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
843    h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
844    sij, elij)
845USE cvparam_mod_h
846    USE cvthermo_mod_h
847  IMPLICIT NONE
848
849
850  ! inputs:
851  INTEGER ncum, nd, nloc
852  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
853  REAL cbmf(nloc), qnk(nloc)
854  REAL ph(nloc, nd+1)
855  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
856  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
857  REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
858
859  ! outputs:
860  INTEGER nent(nloc, nd)
861  REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
862  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
863  REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
864
865  ! local variables:
866  INTEGER i, j, k, ij
867  INTEGER num1, num2
868  REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
869  REAL alt, qp1, smid, sjmin, sjmax, delp, delm
870  REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
871  REAL bsum(nloc, nd)
872  LOGICAL lwork(nloc)
873
874  ! =====================================================================
875  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
876  ! =====================================================================
877
878  DO i = 1, ncum*nlp
879    nent(i, 1) = 0
880    m(i, 1) = 0.0
881  END DO
882
883  DO k = 1, nlp
884    DO j = 1, nlp
885      DO i = 1, ncum
886        qent(i, k, j) = q(i, j)
887        uent(i, k, j) = u(i, j)
888        vent(i, k, j) = v(i, j)
889        elij(i, k, j) = 0.0
890        ment(i, k, j) = 0.0
891        sij(i, k, j) = 0.0
892      END DO
893    END DO
894  END DO
895
896  ! -------------------------------------------------------------------
897  ! --- Calculate rates of mixing,  m(i)
898  ! -------------------------------------------------------------------
899
900  CALL zilch(work, ncum)
901
902  DO j = minorig + 1, nl
903    DO i = 1, ncum
904      IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
905        k = min(j, inb1(i))
906        dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
907          entp*0.04*(ph(i,k)-ph(i,k+1))
908        work(i) = work(i) + dbo
909        m(i, j) = cbmf(i)*dbo
910      END IF
911    END DO
912  END DO
913  DO k = minorig + 1, nl
914    DO i = 1, ncum
915      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
916        m(i, k) = m(i, k)/work(i)
917      END IF
918    END DO
919  END DO
920
921
922  ! =====================================================================
923  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
924  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
925  ! --- FRACTION (sij)
926  ! =====================================================================
927
928
929  DO i = minorig + 1, nl
930    DO j = minorig + 1, nl
931      DO ij = 1, ncum
932        IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
933            inb(ij))) THEN
934          qti = qnk(ij) - ep(ij, i)*clw(ij, i)
935          bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd)
936          anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
937          denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
938          dei = denom
939          IF (abs(dei)<0.01) dei = 0.01
940          sij(ij, i, j) = anum/dei
941          sij(ij, i, i) = 1.0
942          altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
943          altem = altem/bf2
944          cwat = clw(ij, j)*(1.-ep(ij,j))
945          stemp = sij(ij, i, j)
946          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
947            anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
948            denom = denom + lv(ij, j)*(q(ij,i)-qti)
949            IF (abs(denom)<0.01) denom = 0.01
950            sij(ij, i, j) = anum/denom
951            altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
952            altem = altem - (bf2-1.)*cwat
953          END IF
954          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
955            qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
956            uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
957              (1.-sij(ij,i,j))*u(ij, nk(ij))
958            vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
959              (1.-sij(ij,i,j))*v(ij, nk(ij))
960            elij(ij, i, j) = altem
961            elij(ij, i, j) = max(0.0, elij(ij,i,j))
962            ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
963            nent(ij, i) = nent(ij, i) + 1
964          END IF
965          sij(ij, i, j) = max(0.0, sij(ij,i,j))
966          sij(ij, i, j) = min(1.0, sij(ij,i,j))
967        END IF
968      END DO
969    END DO
970
971    ! ***   If no air can entrain at level i assume that updraft detrains
972    ! ***
973    ! ***   at that level and calculate detrained air flux and properties
974    ! ***
975
976    DO ij = 1, ncum
977      IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
978        ment(ij, i, i) = m(ij, i)
979        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
980        uent(ij, i, i) = u(ij, nk(ij))
981        vent(ij, i, i) = v(ij, nk(ij))
982        elij(ij, i, i) = clw(ij, i)
983        sij(ij, i, i) = 1.0
984      END IF
985    END DO
986  END DO
987
988  DO i = 1, ncum
989    sij(i, inb(i), inb(i)) = 1.0
990  END DO
991
992  ! =====================================================================
993  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
994  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
995  ! =====================================================================
996
997  CALL zilch(bsum, ncum*nlp)
998  DO ij = 1, ncum
999    lwork(ij) = .FALSE.
1000  END DO
1001  DO i = minorig + 1, nl
1002
1003    num1 = 0
1004    DO ij = 1, ncum
1005      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
1006    END DO
1007    IF (num1<=0) GO TO 789
1008
1009    DO ij = 1, ncum
1010      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
1011        lwork(ij) = (nent(ij,i)/=0)
1012        qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1013        anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
1014        denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
1015        IF (abs(denom)<0.01) denom = 0.01
1016        scrit(ij) = anum/denom
1017        alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
1018        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
1019        asij(ij) = 0.0
1020        smin(ij) = 1.0
1021      END IF
1022    END DO
1023    DO j = minorig, nl
1024
1025      num2 = 0
1026      DO ij = 1, ncum
1027        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1028          ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
1029      END DO
1030      IF (num2<=0) GO TO 783
1031
1032      DO ij = 1, ncum
1033        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1034            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1035          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
1036            IF (j>i) THEN
1037              smid = min(sij(ij,i,j), scrit(ij))
1038              sjmax = smid
1039              sjmin = smid
1040              IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
1041                smin(ij) = smid
1042                sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
1043                sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
1044                sjmin = min(sjmin, scrit(ij))
1045              END IF
1046            ELSE
1047              sjmax = max(sij(ij,i,j+1), scrit(ij))
1048              smid = max(sij(ij,i,j), scrit(ij))
1049              sjmin = 0.0
1050              IF (j>1) sjmin = sij(ij, i, j-1)
1051              sjmin = max(sjmin, scrit(ij))
1052            END IF
1053            delp = abs(sjmax-smid)
1054            delm = abs(sjmin-smid)
1055            asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
1056            ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
1057          END IF
1058        END IF
1059      END DO
1060783 END DO
1061    DO ij = 1, ncum
1062      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
1063        asij(ij) = max(1.0E-21, asij(ij))
1064        asij(ij) = 1.0/asij(ij)
1065        bsum(ij, i) = 0.0
1066      END IF
1067    END DO
1068    DO j = minorig, nl + 1
1069      DO ij = 1, ncum
1070        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1071            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1072          ment(ij, i, j) = ment(ij, i, j)*asij(ij)
1073          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
1074        END IF
1075      END DO
1076    END DO
1077    DO ij = 1, ncum
1078      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
1079          i)<1.0E-18) .AND. lwork(ij)) THEN
1080        nent(ij, i) = 0
1081        ment(ij, i, i) = m(ij, i)
1082        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1083        uent(ij, i, i) = u(ij, nk(ij))
1084        vent(ij, i, i) = v(ij, nk(ij))
1085        elij(ij, i, i) = clw(ij, i)
1086        sij(ij, i, i) = 1.0
1087      END IF
1088    END DO
1089789 END DO
1090
1091  RETURN
1092END SUBROUTINE cv_mixing
1093
1094SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
1095    ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
1096USE cvparam_mod_h
1097    USE cvthermo_mod_h
1098  IMPLICIT NONE
1099
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.