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

Last change on this file since 5278 was 5276, checked in by abarral, 24 hours ago

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