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

Last change on this file since 5300 was 5285, checked in by abarral, 5 weeks ago

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