LMDZ
cv_routines.F90
Go to the documentation of this file.
1 
2 ! $Id: cv_routines.F90 2311 2015-06-25 07:45:24Z emillour $
3 
4 SUBROUTINE cv_param(nd)
5  IMPLICIT NONE
6 
7  ! ------------------------------------------------------------
8  ! Set parameters for convectL
9  ! (includes microphysical parameters and parameters that
10  ! control the rate of approach to quasi-equilibrium)
11  ! ------------------------------------------------------------
12 
13  ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
14  ! *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- ***
15  ! *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO ***
16  ! *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY ***
17  ! *** BETWEEN 0 C AND TLCRIT) ***
18  ! *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT ***
19  ! *** FORMULATION ***
20  ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT ***
21  ! *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
22  ! *** OF CLOUD ***
23  ! *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN ***
24  ! *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW ***
25  ! *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
26  ! *** OF RAIN ***
27  ! *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
28  ! *** OF SNOW ***
29  ! *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM ***
30  ! *** TRANSPORT ***
31  ! *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION ***
32  ! *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC ***
33  ! *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF ***
34  ! *** APPROACH TO QUASI-EQUILIBRIUM ***
35  ! *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) ***
36  ! *** (DAMP MUST BE LESS THAN 1) ***
37 
38  include "cvparam.h"
39  INTEGER nd
40  CHARACTER (LEN=20) :: modname = 'cv_routines'
41  CHARACTER (LEN=80) :: abort_message
42 
43  ! noff: integer limit for convection (nd-noff)
44  ! minorig: First level of convection
45 
46  noff = 2
47  minorig = 2
48 
49  nl = nd - noff
50  nlp = nl + 1
51  nlm = nl - 1
52 
53  elcrit = 0.0011
54  tlcrit = -55.0
55  entp = 1.5
56  sigs = 0.12
57  sigd = 0.05
58  omtrain = 50.0
59  omtsnow = 5.5
60  coeffr = 1.0
61  coeffs = 0.8
62  dtmax = 0.9
63 
64  cu = 0.70
65 
66  betad = 10.0
67 
68  damp = 0.1
69  alpha = 0.2
70 
71  delta = 0.01 ! cld
72 
73  RETURN
74 END SUBROUTINE cv_param
75 
76 SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
77  IMPLICIT NONE
78 
79  ! =====================================================================
80  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
81  ! =====================================================================
82 
83  ! inputs:
84  INTEGER len, nd, ndp1
85  REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
86 
87  ! outputs:
88  REAL lv(len, nd), cpn(len, nd), tv(len, nd)
89  REAL gz(len, nd), h(len, nd), hm(len, nd)
90 
91  ! local variables:
92  INTEGER k, i
93  REAL cpx(len, nd)
94 
95  include "cvthermo.h"
96  include "cvparam.h"
97 
98 
99  DO k = 1, nlp
100  DO i = 1, len
101  lv(i, k) = lv0 - clmcpv*(t(i,k)-t0)
102  cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
103  cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
104  tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
105  END DO
106  END DO
107 
108  ! gz = phi at the full levels (same as p).
109 
110  DO i = 1, len
111  gz(i, 1) = 0.0
112  END DO
113  DO k = 2, nlp
114  DO i = 1, len
115  gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
116  k)
117  END DO
118  END DO
119 
120  ! h = phi + cpT (dry static energy).
121  ! hm = phi + cp(T-Tbase)+Lq
122 
123  DO k = 1, nlp
124  DO i = 1, len
125  h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
126  hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
127  END DO
128  END DO
129 
130  RETURN
131 END SUBROUTINE cv_prelim
132 
133 SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, &
134  qnk, gznk, plcl)
135  IMPLICIT NONE
136 
137  ! ================================================================
138  ! Purpose: CONVECTIVE FEED
139  ! ================================================================
140 
141  include "cvparam.h"
142 
143  ! inputs:
144  INTEGER len, nd
145  REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
146  REAL hm(len, nd), gz(len, nd)
147 
148  ! outputs:
149  INTEGER iflag(len), nk(len), icb(len), icbmax
150  REAL tnk(len), qnk(len), gznk(len), plcl(len)
151 
152  ! local variables:
153  INTEGER i, k
154  INTEGER ihmin(len)
155  REAL work(len)
156  REAL pnk(len), qsnk(len), rh(len), chi(len)
157 
158  ! -------------------------------------------------------------------
159  ! --- Find level of minimum moist static energy
160  ! --- If level of minimum moist static energy coincides with
161  ! --- or is lower than minimum allowable parcel origin level,
162  ! --- set iflag to 6.
163  ! -------------------------------------------------------------------
164 
165  DO i = 1, len
166  work(i) = 1.0e12
167  ihmin(i) = nl
168  END DO
169  DO k = 2, nlp
170  DO i = 1, len
171  IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
172  work(i) = hm(i, k)
173  ihmin(i) = k
174  END IF
175  END DO
176  END DO
177  DO i = 1, len
178  ihmin(i) = min(ihmin(i), nlm)
179  IF (ihmin(i)<=minorig) THEN
180  iflag(i) = 6
181  END IF
182  END DO
183 
184  ! -------------------------------------------------------------------
185  ! --- Find that model level below the level of minimum moist static
186  ! --- energy that has the maximum value of moist static energy
187  ! -------------------------------------------------------------------
188 
189  DO i = 1, len
190  work(i) = hm(i, minorig)
191  nk(i) = minorig
192  END DO
193  DO k = minorig + 1, nl
194  DO i = 1, len
195  IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
196  work(i) = hm(i, k)
197  nk(i) = k
198  END IF
199  END DO
200  END DO
201  ! -------------------------------------------------------------------
202  ! --- Check whether parcel level temperature and specific humidity
203  ! --- are reasonable
204  ! -------------------------------------------------------------------
205  DO i = 1, len
206  IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
207  400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
208  END DO
209  ! -------------------------------------------------------------------
210  ! --- Calculate lifted condensation level of air at parcel origin level
211  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
212  ! -------------------------------------------------------------------
213  DO i = 1, len
214  tnk(i) = t(i, nk(i))
215  qnk(i) = q(i, nk(i))
216  gznk(i) = gz(i, nk(i))
217  pnk(i) = p(i, nk(i))
218  qsnk(i) = qs(i, nk(i))
219 
220  rh(i) = qnk(i)/qsnk(i)
221  rh(i) = min(1.0, rh(i))
222  chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
223  plcl(i) = pnk(i)*(rh(i)**chi(i))
224  IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
225  ) = 8
226  END DO
227  ! -------------------------------------------------------------------
228  ! --- Calculate first level above lcl (=icb)
229  ! -------------------------------------------------------------------
230  DO i = 1, len
231  icb(i) = nlm
232  END DO
233 
234  DO k = minorig, nl
235  DO i = 1, len
236  IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
237  END DO
238  END DO
239 
240  DO i = 1, len
241  IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
242  END DO
243 
244  ! Compute icbmax.
245 
246  icbmax = 2
247  DO i = 1, len
248  icbmax = max(icbmax, icb(i))
249  END DO
250 
251  RETURN
252 END SUBROUTINE cv_feed
253 
254 SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, &
255  clw)
256  IMPLICIT NONE
257 
258  include "cvthermo.h"
259  include "cvparam.h"
260 
261  ! inputs:
262  INTEGER len, nd
263  INTEGER nk(len), icb(len), icbmax
264  REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
265  REAL p(len, nd)
266 
267  ! outputs:
268  REAL tp(len, nd), tvp(len, nd), clw(len, nd)
269 
270  ! local variables:
271  INTEGER i, k
272  REAL tg, qg, alv, s, ahg, tc, denom, es, rg
273  REAL ah0(len), cpp(len)
274  REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
275 
276  ! -------------------------------------------------------------------
277  ! --- Calculates the lifted parcel virtual temperature at nk,
278  ! --- the actual temperature, and the adiabatic
279  ! --- liquid water content. The procedure is to solve the equation.
280  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
281  ! -------------------------------------------------------------------
282 
283  DO i = 1, len
284  tnk(i) = t(i, nk(i))
285  qnk(i) = q(i, nk(i))
286  gznk(i) = gz(i, nk(i))
287  ticb(i) = t(i, icb(i))
288  gzicb(i) = gz(i, icb(i))
289  END DO
290 
291  ! *** Calculate certain parcel quantities, including static energy ***
292 
293  DO i = 1, len
294  ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
295  273.15)) + gznk(i)
296  cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
297  END DO
298 
299  ! *** Calculate lifted parcel quantities below cloud base ***
300 
301  DO k = minorig, icbmax - 1
302  DO i = 1, len
303  tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
304  tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
305  END DO
306  END DO
307 
308  ! *** Find lifted parcel quantities above cloud base ***
309 
310  DO i = 1, len
311  tg = ticb(i)
312  qg = qs(i, icb(i))
313  alv = lv0 - clmcpv*(ticb(i)-t0)
314 
315  ! First iteration.
316 
317  s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
318  s = 1./s
319  ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
320  tg = tg + s*(ah0(i)-ahg)
321  tg = max(tg, 35.0)
322  tc = tg - t0
323  denom = 243.5 + tc
324  IF (tc>=0.0) THEN
325  es = 6.112*exp(17.67*tc/denom)
326  ELSE
327  es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
328  END IF
329  qg = eps*es/(p(i,icb(i))-es*(1.-eps))
330 
331  ! Second iteration.
332 
333  s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
334  s = 1./s
335  ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
336  tg = tg + s*(ah0(i)-ahg)
337  tg = max(tg, 35.0)
338  tc = tg - t0
339  denom = 243.5 + tc
340  IF (tc>=0.0) THEN
341  es = 6.112*exp(17.67*tc/denom)
342  ELSE
343  es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
344  END IF
345  qg = eps*es/(p(i,icb(i))-es*(1.-eps))
346 
347  alv = lv0 - clmcpv*(ticb(i)-273.15)
348  tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
349  clw(i, icb(i)) = qnk(i) - qg
350  clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
351  rg = qg/(1.-qnk(i))
352  tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
353  END DO
354 
355  DO k = minorig, icbmax
356  DO i = 1, len
357  tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
358  END DO
359  END DO
360 
361  RETURN
362 END SUBROUTINE cv_undilute1
363 
364 SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
365  IMPLICIT NONE
366 
367  ! -------------------------------------------------------------------
368  ! --- Test for instability.
369  ! --- If there was no convection at last time step and parcel
370  ! --- is stable at icb, then set iflag to 4.
371  ! -------------------------------------------------------------------
372 
373  include "cvparam.h"
374 
375  ! inputs:
376  INTEGER len, nd, icb(len)
377  REAL cbmf(len), tv(len, nd), tvp(len, nd)
378 
379  ! outputs:
380  INTEGER iflag(len) ! also an input
381 
382  ! local variables:
383  INTEGER i
384 
385 
386  DO i = 1, len
387  IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
388  icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
389  END DO
390 
391  RETURN
392 END SUBROUTINE cv_trigger
393 
394 SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
395  tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
396  tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
397  v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
398  USE print_control_mod, ONLY: lunout
399  IMPLICIT NONE
400 
401  include "cvparam.h"
402 
403  ! inputs:
404  INTEGER len, ncum, nd, nloc
405  INTEGER iflag1(len), nk1(len), icb1(len)
406  REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
407  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
408  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
409  REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
410  REAL tvp1(len, nd), clw1(len, nd)
411 
412  ! outputs:
413  INTEGER iflag(nloc), nk(nloc), icb(nloc)
414  REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
415  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
416  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
417  REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
418  REAL tvp(nloc, nd), clw(nloc, nd)
419  REAL dph(nloc, nd)
420 
421  ! local variables:
422  INTEGER i, k, nn
423  CHARACTER (LEN=20) :: modname = 'cv_compress'
424  CHARACTER (LEN=80) :: abort_message
425 
426 
427  DO k = 1, nl + 1
428  nn = 0
429  DO i = 1, len
430  IF (iflag1(i)==0) THEN
431  nn = nn + 1
432  t(nn, k) = t1(i, k)
433  q(nn, k) = q1(i, k)
434  qs(nn, k) = qs1(i, k)
435  u(nn, k) = u1(i, k)
436  v(nn, k) = v1(i, k)
437  gz(nn, k) = gz1(i, k)
438  h(nn, k) = h1(i, k)
439  lv(nn, k) = lv1(i, k)
440  cpn(nn, k) = cpn1(i, k)
441  p(nn, k) = p1(i, k)
442  ph(nn, k) = ph1(i, k)
443  tv(nn, k) = tv1(i, k)
444  tp(nn, k) = tp1(i, k)
445  tvp(nn, k) = tvp1(i, k)
446  clw(nn, k) = clw1(i, k)
447  END IF
448  END DO
449  END DO
450 
451  IF (nn/=ncum) THEN
452  WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
453  abort_message = ''
454  CALL abort_physic(modname, abort_message, 1)
455  END IF
456 
457  nn = 0
458  DO i = 1, len
459  IF (iflag1(i)==0) THEN
460  nn = nn + 1
461  cbmf(nn) = cbmf1(i)
462  plcl(nn) = plcl1(i)
463  tnk(nn) = tnk1(i)
464  qnk(nn) = qnk1(i)
465  gznk(nn) = gznk1(i)
466  nk(nn) = nk1(i)
467  icb(nn) = icb1(i)
468  iflag(nn) = iflag1(i)
469  END IF
470  END DO
471 
472  DO k = 1, nl
473  DO i = 1, ncum
474  dph(i, k) = ph(i, k) - ph(i, k+1)
475  END DO
476  END DO
477 
478  RETURN
479 END SUBROUTINE cv_compress
480 
481 SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
482  gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
483  IMPLICIT NONE
484 
485  ! ---------------------------------------------------------------------
486  ! Purpose:
487  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
488  ! &
489  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
490  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
491  ! &
492  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
493  ! ---------------------------------------------------------------------
494 
495  include "cvthermo.h"
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
762 END SUBROUTINE cv_undilute2
763 
764 SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
765  cpn, iflag, cbmf)
766  IMPLICIT NONE
767 
768  ! inputs:
769  INTEGER ncum, nd, nloc
770  INTEGER nk(nloc), icb(nloc)
771  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
772  REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent...
773  REAL plcl(nloc), cpn(nloc, nd)
774 
775  ! outputs:
776  INTEGER iflag(nloc)
777  REAL cbmf(nloc) ! also an input
778 
779  ! local variables:
780  INTEGER i, k, icbmax
781  REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
782  REAL work(nloc)
783 
784  include "cvthermo.h"
785  include "cvparam.h"
786 
787  ! -------------------------------------------------------------------
788  ! Compute icbmax.
789  ! -------------------------------------------------------------------
790 
791  icbmax = 2
792  DO i = 1, ncum
793  icbmax = max(icbmax, icb(i))
794  END DO
795 
796  ! =====================================================================
797  ! --- CALCULATE CLOUD BASE MASS FLUX
798  ! =====================================================================
799 
800  ! tvpplcl = parcel temperature lifted adiabatically from level
801  ! icb-1 to the LCL.
802  ! tvaplcl = virtual temperature at the LCL.
803 
804  DO i = 1, ncum
805  dtpbl(i) = 0.0
806  tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( &
807  i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
808  tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
809  ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
810  END DO
811 
812  ! -------------------------------------------------------------------
813  ! --- Interpolate difference between lifted parcel and
814  ! --- environmental temperatures to lifted condensation level
815  ! -------------------------------------------------------------------
816 
817  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
818 
819  DO k = minorig, icbmax
820  DO i = 1, ncum
821  IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
822  dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
823  END IF
824  END DO
825  END DO
826  DO i = 1, ncum
827  dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
828  dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
829  END DO
830 
831  ! -------------------------------------------------------------------
832  ! --- Adjust cloud base mass flux
833  ! -------------------------------------------------------------------
834 
835  DO i = 1, ncum
836  work(i) = cbmf(i)
837  cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
838  IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
839  iflag(i) = 3
840  END IF
841  END DO
842 
843  RETURN
844 END SUBROUTINE cv_closure
845 
846 SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
847  h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
848  sij, elij)
849  IMPLICIT NONE
850 
851  include "cvthermo.h"
852  include "cvparam.h"
853 
854  ! inputs:
855  INTEGER ncum, nd, nloc
856  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
857  REAL cbmf(nloc), qnk(nloc)
858  REAL ph(nloc, nd+1)
859  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
860  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
861  REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
862 
863  ! outputs:
864  INTEGER nent(nloc, nd)
865  REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
866  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
867  REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
868 
869  ! local variables:
870  INTEGER i, j, k, ij
871  INTEGER num1, num2
872  REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
873  REAL alt, qp1, smid, sjmin, sjmax, delp, delm
874  REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
875  REAL bsum(nloc, nd)
876  LOGICAL lwork(nloc)
877 
878  ! =====================================================================
879  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
880  ! =====================================================================
881 
882  DO i = 1, ncum*nlp
883  nent(i, 1) = 0
884  m(i, 1) = 0.0
885  END DO
886 
887  DO k = 1, nlp
888  DO j = 1, nlp
889  DO i = 1, ncum
890  qent(i, k, j) = q(i, j)
891  uent(i, k, j) = u(i, j)
892  vent(i, k, j) = v(i, j)
893  elij(i, k, j) = 0.0
894  ment(i, k, j) = 0.0
895  sij(i, k, j) = 0.0
896  END DO
897  END DO
898  END DO
899 
900  ! -------------------------------------------------------------------
901  ! --- Calculate rates of mixing, m(i)
902  ! -------------------------------------------------------------------
903 
904  CALL zilch(work, ncum)
905 
906  DO j = minorig + 1, nl
907  DO i = 1, ncum
908  IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
909  k = min(j, inb1(i))
910  dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
911  entp*0.04*(ph(i,k)-ph(i,k+1))
912  work(i) = work(i) + dbo
913  m(i, j) = cbmf(i)*dbo
914  END IF
915  END DO
916  END DO
917  DO k = minorig + 1, nl
918  DO i = 1, ncum
919  IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
920  m(i, k) = m(i, k)/work(i)
921  END IF
922  END DO
923  END DO
924 
925 
926  ! =====================================================================
927  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
928  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
929  ! --- FRACTION (sij)
930  ! =====================================================================
931 
932 
933  DO i = minorig + 1, nl
934  DO j = minorig + 1, nl
935  DO ij = 1, ncum
936  IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
937  inb(ij))) THEN
938  qti = qnk(ij) - ep(ij, i)*clw(ij, i)
939  bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd)
940  anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
941  denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
942  dei = denom
943  IF (abs(dei)<0.01) dei = 0.01
944  sij(ij, i, j) = anum/dei
945  sij(ij, i, i) = 1.0
946  altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
947  altem = altem/bf2
948  cwat = clw(ij, j)*(1.-ep(ij,j))
949  stemp = sij(ij, i, j)
950  IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
951  anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
952  denom = denom + lv(ij, j)*(q(ij,i)-qti)
953  IF (abs(denom)<0.01) denom = 0.01
954  sij(ij, i, j) = anum/denom
955  altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
956  altem = altem - (bf2-1.)*cwat
957  END IF
958  IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
959  qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
960  uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
961  (1.-sij(ij,i,j))*u(ij, nk(ij))
962  vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
963  (1.-sij(ij,i,j))*v(ij, nk(ij))
964  elij(ij, i, j) = altem
965  elij(ij, i, j) = max(0.0, elij(ij,i,j))
966  ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
967  nent(ij, i) = nent(ij, i) + 1
968  END IF
969  sij(ij, i, j) = max(0.0, sij(ij,i,j))
970  sij(ij, i, j) = min(1.0, sij(ij,i,j))
971  END IF
972  END DO
973  END DO
974 
975  ! *** If no air can entrain at level i assume that updraft detrains
976  ! ***
977  ! *** at that level and calculate detrained air flux and properties
978  ! ***
979 
980  DO ij = 1, ncum
981  IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
982  ment(ij, i, i) = m(ij, i)
983  qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
984  uent(ij, i, i) = u(ij, nk(ij))
985  vent(ij, i, i) = v(ij, nk(ij))
986  elij(ij, i, i) = clw(ij, i)
987  sij(ij, i, i) = 1.0
988  END IF
989  END DO
990  END DO
991 
992  DO i = 1, ncum
993  sij(i, inb(i), inb(i)) = 1.0
994  END DO
995 
996  ! =====================================================================
997  ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
998  ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
999  ! =====================================================================
1000 
1001  CALL zilch(bsum, ncum*nlp)
1002  DO ij = 1, ncum
1003  lwork(ij) = .false.
1004  END DO
1005  DO i = minorig + 1, nl
1006 
1007  num1 = 0
1008  DO ij = 1, ncum
1009  IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
1010  END DO
1011  IF (num1<=0) GO TO 789
1012 
1013  DO ij = 1, ncum
1014  IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
1015  lwork(ij) = (nent(ij,i)/=0)
1016  qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1017  anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
1018  denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
1019  IF (abs(denom)<0.01) denom = 0.01
1020  scrit(ij) = anum/denom
1021  alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
1022  IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
1023  asij(ij) = 0.0
1024  smin(ij) = 1.0
1025  END IF
1026  END DO
1027  DO j = minorig, nl
1028 
1029  num2 = 0
1030  DO ij = 1, ncum
1031  IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1032  ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
1033  END DO
1034  IF (num2<=0) GO TO 783
1035 
1036  DO ij = 1, ncum
1037  IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1038  ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1039  IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
1040  IF (j>i) THEN
1041  smid = min(sij(ij,i,j), scrit(ij))
1042  sjmax = smid
1043  sjmin = smid
1044  IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
1045  smin(ij) = smid
1046  sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
1047  sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
1048  sjmin = min(sjmin, scrit(ij))
1049  END IF
1050  ELSE
1051  sjmax = max(sij(ij,i,j+1), scrit(ij))
1052  smid = max(sij(ij,i,j), scrit(ij))
1053  sjmin = 0.0
1054  IF (j>1) sjmin = sij(ij, i, j-1)
1055  sjmin = max(sjmin, scrit(ij))
1056  END IF
1057  delp = abs(sjmax-smid)
1058  delm = abs(sjmin-smid)
1059  asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
1060  ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
1061  END IF
1062  END IF
1063  END DO
1064 783 END DO
1065  DO ij = 1, ncum
1066  IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
1067  asij(ij) = max(1.0e-21, asij(ij))
1068  asij(ij) = 1.0/asij(ij)
1069  bsum(ij, i) = 0.0
1070  END IF
1071  END DO
1072  DO j = minorig, nl + 1
1073  DO ij = 1, ncum
1074  IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1075  ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1076  ment(ij, i, j) = ment(ij, i, j)*asij(ij)
1077  bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
1078  END IF
1079  END DO
1080  END DO
1081  DO ij = 1, ncum
1082  IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
1083  i)<1.0e-18) .AND. lwork(ij)) THEN
1084  nent(ij, i) = 0
1085  ment(ij, i, i) = m(ij, i)
1086  qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1087  uent(ij, i, i) = u(ij, nk(ij))
1088  vent(ij, i, i) = v(ij, nk(ij))
1089  elij(ij, i, i) = clw(ij, i)
1090  sij(ij, i, i) = 1.0
1091  END IF
1092  END DO
1093 789 END DO
1094 
1095  RETURN
1096 END SUBROUTINE cv_mixing
1097 
1098 SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
1099  ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
1100  IMPLICIT NONE
1101 
1102 
1103  include "cvthermo.h"
1104  include "cvparam.h"
1105 
1106  ! inputs:
1107  INTEGER ncum, nd, nloc
1108  INTEGER inb(nloc)
1109  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
1110  REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
1111  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1112  REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
1113  REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
1114 
1115  ! outputs:
1116  INTEGER iflag(nloc) ! also an input
1117  REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd)
1118  REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd)
1119 
1120  ! local variables:
1121  INTEGER i, j, k, ij, num1
1122  INTEGER jtt(nloc)
1123  REAL awat, coeff, qsm, afac, sigt, b6, c6, revap
1124  REAL dhdp, fac, qstm, rat
1125  REAL wdtrain(nloc)
1126  LOGICAL lwork(nloc)
1127 
1128  ! =====================================================================
1129  ! --- PRECIPITATING DOWNDRAFT CALCULATION
1130  ! =====================================================================
1131 
1132  ! Initializations:
1133 
1134  DO i = 1, ncum
1135  DO k = 1, nl + 1
1136  wt(i, k) = omtsnow
1137  mp(i, k) = 0.0
1138  evap(i, k) = 0.0
1139  water(i, k) = 0.0
1140  END DO
1141  END DO
1142 
1143  DO i = 1, ncum
1144  qp(i, 1) = q(i, 1)
1145  up(i, 1) = u(i, 1)
1146  vp(i, 1) = v(i, 1)
1147  END DO
1148 
1149  DO k = 2, nl + 1
1150  DO i = 1, ncum
1151  qp(i, k) = q(i, k-1)
1152  up(i, k) = u(i, k-1)
1153  vp(i, k) = v(i, k-1)
1154  END DO
1155  END DO
1156 
1157 
1158  ! *** Check whether ep(inb)=0, if so, skip precipitating ***
1159  ! *** downdraft calculation ***
1160 
1161 
1162  ! *** Integrate liquid water equation to find condensed water ***
1163  ! *** and condensed water flux ***
1164 
1165 
1166  DO i = 1, ncum
1167  jtt(i) = 2
1168  IF (ep(i,inb(i))<=0.0001) iflag(i) = 2
1169  IF (iflag(i)==0) THEN
1170  lwork(i) = .true.
1171  ELSE
1172  lwork(i) = .false.
1173  END IF
1174  END DO
1175 
1176  ! *** Begin downdraft loop ***
1177 
1178 
1179  CALL zilch(wdtrain, ncum)
1180  DO i = nl + 1, 1, -1
1181 
1182  num1 = 0
1183  DO ij = 1, ncum
1184  IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
1185  END DO
1186  IF (num1<=0) GO TO 899
1187 
1188 
1189  ! *** Calculate detrained precipitation ***
1190 
1191  DO ij = 1, ncum
1192  IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1193  wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i)
1194  END IF
1195  END DO
1196 
1197  IF (i>1) THEN
1198  DO j = 1, i - 1
1199  DO ij = 1, ncum
1200  IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1201  awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i)
1202  awat = max(0.0, awat)
1203  wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i)
1204  END IF
1205  END DO
1206  END DO
1207  END IF
1208 
1209  ! *** Find rain water and evaporation using provisional ***
1210  ! *** estimates of qp(i)and qp(i-1) ***
1211 
1212 
1213  ! *** Value of terminal velocity and coeffecient of evaporation for snow
1214  ! ***
1215 
1216  DO ij = 1, ncum
1217  IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1218  coeff = coeffs
1219  wt(ij, i) = omtsnow
1220 
1221  ! *** Value of terminal velocity and coeffecient of evaporation for
1222  ! rain ***
1223 
1224  IF (t(ij,i)>273.0) THEN
1225  coeff = coeffr
1226  wt(ij, i) = omtrain
1227  END IF
1228  qsm = 0.5*(q(ij,i)+qp(ij,i+1))
1229  afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0e4+2.0e3*ph(ij,i)*qs(ij,i))
1230  afac = max(afac, 0.0)
1231  sigt = sigp(ij, i)
1232  sigt = max(0.0, sigt)
1233  sigt = min(1.0, sigt)
1234  b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i)
1235  c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i)
1236  revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
1237  evap(ij, i) = sigt*afac*revap
1238  water(ij, i) = revap*revap
1239 
1240  ! *** Calculate precipitating downdraft mass flux under ***
1241  ! *** hydrostatic approximation ***
1242 
1243  IF (i>1) THEN
1244  dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
1245  dhdp = max(dhdp, 10.0)
1246  mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp
1247  mp(ij, i) = max(mp(ij,i), 0.0)
1248 
1249  ! *** Add small amount of inertia to downdraft ***
1250 
1251  fac = 20.0/(ph(ij,i-1)-ph(ij,i))
1252  mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
1253 
1254  ! *** Force mp to decrease linearly to zero
1255  ! ***
1256  ! *** between about 950 mb and the surface
1257  ! ***
1258 
1259  IF (p(ij,i)>(0.949*p(ij,1))) THEN
1260  jtt(ij) = max(jtt(ij), i)
1261  mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ &
1262  (p(ij,1)-p(ij,jtt(ij)))
1263  END IF
1264  END IF
1265 
1266  ! *** Find mixing ratio of precipitating downdraft ***
1267 
1268  IF (i/=inb(ij)) THEN
1269  IF (i==1) THEN
1270  qstm = qs(ij, 1)
1271  ELSE
1272  qstm = qs(ij, i-1)
1273  END IF
1274  IF (mp(ij,i)>mp(ij,i+1)) THEN
1275  rat = mp(ij, i+1)/mp(ij, i)
1276  qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + &
1277  100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
1278  up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat)
1279  vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat)
1280  ELSE
1281  IF (mp(ij,i+1)>0.0) THEN
1282  qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, &
1283  i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, &
1284  i)))/(lv(ij,i)+t(ij,i)*(cl-cpd))
1285  up(ij, i) = up(ij, i+1)
1286  vp(ij, i) = vp(ij, i+1)
1287  END IF
1288  END IF
1289  qp(ij, i) = min(qp(ij,i), qstm)
1290  qp(ij, i) = max(qp(ij,i), 0.0)
1291  END IF
1292  END IF
1293  END DO
1294 899 END DO
1295 
1296  RETURN
1297 END SUBROUTINE cv_unsat
1298 
1299 SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
1300  ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
1301  ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
1302  precip, cbmf, ft, fq, fu, fv, ma, qcondc)
1303  IMPLICIT NONE
1304 
1305  include "cvthermo.h"
1306  include "cvparam.h"
1307 
1308  ! inputs
1309  INTEGER ncum, nd, nloc
1310  INTEGER nk(nloc), icb(nloc), inb(nloc)
1311  INTEGER nent(nloc, nd)
1312  REAL delt
1313  REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
1314  REAL gz(nloc, nd)
1315  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1316  REAL hp(nloc, nd), lv(nloc, nd)
1317  REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
1318  REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd)
1319  REAL up(nloc, nd), vp(nloc, nd)
1320  REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd)
1321  REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd)
1322  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
1323  REAL tv(nloc, nd), tvp(nloc, nd)
1324 
1325  ! outputs
1326  INTEGER iflag(nloc) ! also an input
1327  REAL cbmf(nloc) ! also an input
1328  REAL wd(nloc), tprime(nloc), qprime(nloc)
1329  REAL precip(nloc)
1330  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1331  REAL ma(nloc, nd)
1332  REAL qcondc(nloc, nd)
1333 
1334  ! local variables
1335  INTEGER i, j, ij, k, num1
1336  REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti
1337  REAL work(nloc), am(nloc), amp1(nloc), ad(nloc)
1338  REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd)
1339  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
1340  REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld
1341 
1342 
1343  ! -- initializations:
1344 
1345  delti = 1.0/delt
1346 
1347  DO i = 1, ncum
1348  precip(i) = 0.0
1349  wd(i) = 0.0
1350  tprime(i) = 0.0
1351  qprime(i) = 0.0
1352  DO k = 1, nl + 1
1353  ft(i, k) = 0.0
1354  fu(i, k) = 0.0
1355  fv(i, k) = 0.0
1356  fq(i, k) = 0.0
1357  lvcp(i, k) = lv(i, k)/cpn(i, k)
1358  qcondc(i, k) = 0.0 ! cld
1359  qcond(i, k) = 0.0 ! cld
1360  nqcond(i, k) = 0.0 ! cld
1361  END DO
1362  END DO
1363 
1364 
1365  ! *** Calculate surface precipitation in mm/day ***
1366 
1367  DO i = 1, ncum
1368  IF (iflag(i)<=1) THEN
1369  ! c precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1370  ! c & /(rowl*g)
1371  ! c precip(i)=precip(i)*delt/86400.
1372  precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g
1373  END IF
1374  END DO
1375 
1376 
1377  ! *** Calculate downdraft velocity scale and surface temperature and ***
1378  ! *** water vapor fluctuations ***
1379 
1380  DO i = 1, ncum
1381  wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i)))
1382  qprime(i) = 0.5*(qp(i,1)-q(i,1))
1383  tprime(i) = lv0*qprime(i)/cpd
1384  END DO
1385 
1386  ! *** Calculate tendencies of lowest level potential temperature ***
1387  ! *** and mixing ratio ***
1388 
1389  DO i = 1, ncum
1390  work(i) = 0.01/(ph(i,1)-ph(i,2))
1391  am(i) = 0.0
1392  END DO
1393  DO k = 2, nl
1394  DO i = 1, ncum
1395  IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
1396  am(i) = am(i) + m(i, k)
1397  END IF
1398  END DO
1399  END DO
1400  DO i = 1, ncum
1401  IF ((g*work(i)*am(i))>=delti) iflag(i) = 1
1402  ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, &
1403  1))/cpn(i,1))
1404  ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1)
1405  ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* &
1406  work(i)/cpn(i, 1)
1407  fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + &
1408  sigd*evap(i, 1)
1409  fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i)
1410  fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, &
1411  2)-u(i,1)))
1412  fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, &
1413  2)-v(i,1)))
1414  END DO
1415  DO j = 2, nl
1416  DO i = 1, ncum
1417  IF (j<=inb(i)) THEN
1418  fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1))
1419  fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1))
1420  fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1))
1421  END IF
1422  END DO
1423  END DO
1424 
1425  ! *** Calculate tendencies of potential temperature and mixing ratio ***
1426  ! *** at levels above the lowest level ***
1427 
1428  ! *** First find the net saturated updraft and downdraft mass fluxes ***
1429  ! *** through each level ***
1430 
1431  DO i = 2, nl + 1
1432 
1433  num1 = 0
1434  DO ij = 1, ncum
1435  IF (i<=inb(ij)) num1 = num1 + 1
1436  END DO
1437  IF (num1<=0) GO TO 1500
1438 
1439  CALL zilch(amp1, ncum)
1440  CALL zilch(ad, ncum)
1441 
1442  DO k = i + 1, nl + 1
1443  DO ij = 1, ncum
1444  IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN
1445  amp1(ij) = amp1(ij) + m(ij, k)
1446  END IF
1447  END DO
1448  END DO
1449 
1450  DO k = 1, i
1451  DO j = i + 1, nl + 1
1452  DO ij = 1, ncum
1453  IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN
1454  amp1(ij) = amp1(ij) + ment(ij, k, j)
1455  END IF
1456  END DO
1457  END DO
1458  END DO
1459  DO k = 1, i - 1
1460  DO j = i, nl + 1
1461  DO ij = 1, ncum
1462  IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
1463  ad(ij) = ad(ij) + ment(ij, j, k)
1464  END IF
1465  END DO
1466  END DO
1467  END DO
1468 
1469  DO ij = 1, ncum
1470  IF (i<=inb(ij)) THEN
1471  dpinv = 0.01/(ph(ij,i)-ph(ij,i+1))
1472  cpinv = 1.0/cpn(ij, i)
1473 
1474  ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, &
1475  i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, &
1476  i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i)
1477  ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij &
1478  ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
1479  ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( &
1480  ij,i+1)-t(ij,i))*dpinv*cpinv
1481  fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, &
1482  i))-ad(ij)*(q(ij,i)-q(ij,i-1)))
1483  fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, &
1484  i))-ad(ij)*(u(ij,i)-u(ij,i-1)))
1485  fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, &
1486  i))-ad(ij)*(v(ij,i)-v(ij,i-1)))
1487  END IF
1488  END DO
1489  DO k = 1, i - 1
1490  DO ij = 1, ncum
1491  IF (i<=inb(ij)) THEN
1492  awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i)
1493  awat = max(awat, 0.0)
1494  fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q &
1495  (ij,i))
1496  fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1497  ))
1498  fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1499  ))
1500  ! (saturated updrafts resulting from mixing) ! cld
1501  qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld
1502  nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1503  END IF
1504  END DO
1505  END DO
1506  DO k = i, nl + 1
1507  DO ij = 1, ncum
1508  IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
1509  fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i &
1510  ))
1511  fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1512  ))
1513  fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1514  ))
1515  END IF
1516  END DO
1517  END DO
1518  DO ij = 1, ncum
1519  IF (i<=inb(ij)) THEN
1520  fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, &
1521  i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
1522  fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, &
1523  i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv
1524  fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, &
1525  i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv
1526  ! (saturated downdrafts resulting from mixing) ! cld
1527  DO k = i + 1, inb(ij) ! cld
1528  qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld
1529  nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1530  END DO ! cld
1531  ! (particular case: no detraining level is found) ! cld
1532  IF (nent(ij,i)==0) THEN ! cld
1533  qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld
1534  nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1535  END IF ! cld
1536  IF (nqcond(ij,i)/=0.) THEN ! cld
1537  qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld
1538  END IF ! cld
1539  END IF
1540  END DO
1541 1500 END DO
1542 
1543  ! *** Adjust tendencies at top of convection layer to reflect ***
1544  ! *** actual position of the level zero cape ***
1545 
1546  DO ij = 1, ncum
1547  fqold = fq(ij, inb(ij))
1548  fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij))
1549  fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, &
1550  inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1551  inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1)
1552  ftold = ft(ij, inb(ij))
1553  ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij))
1554  ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, &
1555  inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1556  inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1)
1557  fuold = fu(ij, inb(ij))
1558  fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij))
1559  fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, &
1560  inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1561  fvold = fv(ij, inb(ij))
1562  fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij))
1563  fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, &
1564  inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1565  END DO
1566 
1567  ! *** Very slightly adjust tendencies to force exact ***
1568  ! *** enthalpy, momentum and tracer conservation ***
1569 
1570  DO ij = 1, ncum
1571  ents(ij) = 0.0
1572  uav(ij) = 0.0
1573  vav(ij) = 0.0
1574  DO i = 1, inb(ij)
1575  ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- &
1576  ph(ij,i+1))
1577  uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1))
1578  vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1))
1579  END DO
1580  END DO
1581  DO ij = 1, ncum
1582  ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1583  uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1584  vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1585  END DO
1586  DO ij = 1, ncum
1587  DO i = 1, inb(ij)
1588  ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i)
1589  fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij))
1590  fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij))
1591  END DO
1592  END DO
1593 
1594  DO k = 1, nl + 1
1595  DO i = 1, ncum
1596  IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10
1597  END DO
1598  END DO
1599 
1600 
1601  DO i = 1, ncum
1602  IF (iflag(i)>2) THEN
1603  precip(i) = 0.0
1604  cbmf(i) = 0.0
1605  END IF
1606  END DO
1607  DO k = 1, nl
1608  DO i = 1, ncum
1609  IF (iflag(i)>2) THEN
1610  ft(i, k) = 0.0
1611  fq(i, k) = 0.0
1612  fu(i, k) = 0.0
1613  fv(i, k) = 0.0
1614  qcondc(i, k) = 0.0 ! cld
1615  END IF
1616  END DO
1617  END DO
1618 
1619  DO k = 1, nl + 1
1620  DO i = 1, ncum
1621  ma(i, k) = 0.
1622  END DO
1623  END DO
1624  DO k = nl, 1, -1
1625  DO i = 1, ncum
1626  ma(i, k) = ma(i, k+1) + m(i, k)
1627  END DO
1628  END DO
1629 
1630 
1631  ! *** diagnose the in-cloud mixing ratio *** ! cld
1632  ! *** of condensed water *** ! cld
1633  ! ! cld
1634  DO ij = 1, ncum ! cld
1635  DO i = 1, nd ! cld
1636  mac(ij, i) = 0.0 ! cld
1637  wa(ij, i) = 0.0 ! cld
1638  siga(ij, i) = 0.0 ! cld
1639  END DO ! cld
1640  DO i = nk(ij), inb(ij) ! cld
1641  DO k = i + 1, inb(ij) + 1 ! cld
1642  mac(ij, i) = mac(ij, i) + m(ij, k) ! cld
1643  END DO ! cld
1644  END DO ! cld
1645  DO i = icb(ij), inb(ij) - 1 ! cld
1646  ax(ij, i) = 0. ! cld
1647  DO j = icb(ij), i ! cld
1648  ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld
1649  *(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld
1650  END DO ! cld
1651  IF (ax(ij,i)>0.0) THEN ! cld
1652  wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld
1653  END IF ! cld
1654  END DO ! cld
1655  DO i = 1, nl ! cld
1656  IF (wa(ij,i)>0.0) & ! cld
1657  siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld
1658  *rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld
1659  siga(ij, i) = min(siga(ij,i), 1.0) ! cld
1660  qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld
1661  +(1.-siga(ij,i))*qcond(ij, i) ! cld
1662  END DO ! cld
1663  END DO ! cld
1664 
1665  RETURN
1666 END SUBROUTINE cv_yield
1667 
1668 SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
1669  fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
1670  qcondc1)
1671  IMPLICIT NONE
1672 
1673  include "cvparam.h"
1674 
1675  ! inputs:
1676  INTEGER len, ncum, nd, nloc
1677  INTEGER idcum(nloc)
1678  INTEGER iflag(nloc)
1679  REAL precip(nloc), cbmf(nloc)
1680  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1681  REAL ma(nloc, nd)
1682  REAL qcondc(nloc, nd) !cld
1683 
1684  ! outputs:
1685  INTEGER iflag1(len)
1686  REAL precip1(len), cbmf1(len)
1687  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
1688  REAL ma1(len, nd)
1689  REAL qcondc1(len, nd) !cld
1690 
1691  ! local variables:
1692  INTEGER i, k
1693 
1694  DO i = 1, ncum
1695  precip1(idcum(i)) = precip(i)
1696  cbmf1(idcum(i)) = cbmf(i)
1697  iflag1(idcum(i)) = iflag(i)
1698  END DO
1699 
1700  DO k = 1, nl
1701  DO i = 1, ncum
1702  ft1(idcum(i), k) = ft(i, k)
1703  fq1(idcum(i), k) = fq(i, k)
1704  fu1(idcum(i), k) = fu(i, k)
1705  fv1(idcum(i), k) = fv(i, k)
1706  ma1(idcum(i), k) = ma(i, k)
1707  qcondc1(idcum(i), k) = qcondc(i, k)
1708  END DO
1709  END DO
1710 
1711  RETURN
1712 END SUBROUTINE cv_uncompress
1713 
clmcpv
!$Id!Thermodynamical constants for t0 real clmcpv
Definition: cvthermo.h:6
t0
!$Id!Thermodynamical constants for t0 real clmci real epsim1 real hrd real grav COMMON cvthermo rowl t0
Definition: cvthermo.h:6
cpv
!$Id!Thermodynamical constants for cpv
Definition: cvthermo.h:6
damp
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real coeffs real dtmax real cu real betad real damp real delta COMMON cvparam nlm tlcrit sigd coeffs damp
Definition: cvparam.h:12
minorig
!$Id!Parameters for minorig
Definition: cv30param.h:5
g
!$Id mode_top_bound COMMON comconstr g
Definition: comconst.h:7
cv_param
subroutine cv_param(nd)
Definition: cv_routines.F90:5
dtmax
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real coeffs real dtmax real cu real betad real damp real delta COMMON cvparam nlm tlcrit sigd coeffs dtmax
Definition: cvparam.h:12
cv_closure
subroutine cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, cpn, iflag, cbmf)
Definition: cv_routines.F90:766
delta
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real alpha real delta real betad COMMON cv30param nlm spfac &!IM cf ptcrit omtrain dttrig delta
Definition: cv30param.h:5
hrd
!$Id!Thermodynamical constants for t0 real clmci real epsim1 real hrd real grav COMMON cvthermo rowl cpvmcl hrd
Definition: cvthermo.h:6
sigs
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigs
Definition: cvparam.h:12
cv_mixing
subroutine cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, sij, elij)
Definition: cv_routines.F90:849
ginv
!$Id!Thermodynamical constants for t0 real clmci real epsim1 real ginv
Definition: cvthermo.h:6
eps
!$Id!Thermodynamical constants for t0 real clmci real eps
Definition: cvthermo.h:6
cv_compress
subroutine cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
Definition: cv_routines.F90:398
cv_prelim
subroutine cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
Definition: cv_routines.F90:77
omtrain
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real omtrain
Definition: cvparam.h:12
nlp
!$Id!Parameters for nlp
Definition: cv30param.h:5
rrd
!$Id!Thermodynamical constants for rrd
Definition: cvthermo.h:6
omtsnow
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real omtsnow
Definition: cvparam.h:12
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
cv_trigger
subroutine cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
Definition: cv_routines.F90:365
cv_undilute2
subroutine cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
Definition: cv_routines.F90:483
cv_unsat
subroutine cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
Definition: cv_routines.F90:1100
cl
!$Id!Thermodynamical constants for cl
Definition: cvthermo.h:6
noff
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real alpha real delta real betad COMMON cv30param noff
Definition: cv30param.h:5
betad
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real coeffs real dtmax real cu real betad real damp real delta COMMON cvparam nlm tlcrit sigd coeffs betad
Definition: cvparam.h:12
nl
!$Id!Parameters for nl
Definition: cv30param.h:5
elcrit
!Parameters for nlm real spfac integer flag_wb real ptcrit real elcrit
Definition: cv3param.h:2
u
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm u(l)
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
cu
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real coeffs real dtmax real cu real betad real damp real delta COMMON cvparam nlm tlcrit sigd coeffs cu
Definition: cvparam.h:12
alpha
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real alpha real delta real betad COMMON cv30param nlm spfac &!IM cf ptcrit omtrain dttrig alpha
Definition: cv30param.h:5
print_control_mod
Definition: print_control_mod.F90:2
cv_undilute1
subroutine cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, clw)
Definition: cv_routines.F90:256
cv_yield
subroutine cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, precip, cbmf, ft, fq, fu, fv, ma, qcondc)
Definition: cv_routines.F90:1303
epsim1
!$Id!Thermodynamical constants for t0 real clmci real epsim1 real hrd real grav COMMON cvthermo rowl cpvmcl epsim1
Definition: cvthermo.h:6
cv_uncompress
subroutine cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, qcondc1)
Definition: cv_routines.F90:1671
cv_feed
subroutine cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, qnk, gznk, plcl)
Definition: cv_routines.F90:135
cpd
!$Id!Thermodynamical constants for t0 real clmci real epsim1 real hrd real grav COMMON cvthermo cpd
Definition: cvthermo.h:6
rrv
!$Id!Thermodynamical constants for rrv
Definition: cvthermo.h:6
abort_physic
subroutine abort_physic(modname, message, ierr)
Definition: abort_physic.F90:3
lv0
!$Id!Thermodynamical constants for lv0
Definition: cvthermo.h:6
entp
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real coeffs real dtmax real cu real betad real damp real delta COMMON cvparam nlm tlcrit entp
Definition: cvparam.h:12
epsi
!$Id!Thermodynamical constants for t0 real clmci real epsi
Definition: cvthermo.h:6
coeffr
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real coeffr
Definition: cvparam.h:12
sigd
!$Id!Parameters for nlm real sigd
Definition: cv30param.h:5
lunout
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7
zilch
subroutine zilch(x, m)
Definition: zilch.F90:5