GFS Operational Physics Documentation  Revision: 81451
sascnvn.f
Go to the documentation of this file.
1 
11 
14 
56  subroutine sascnvn(im,ix,km,jcap,delt,delp,prslp,psp,phil,ql,
57  & q1,t1,u1,v1,cldwrk,rn,kbot,ktop,kcnv,islimsk,
58  & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc)
59 ! & q1,t1,u1,v1,rcs,cldwrk,rn,kbot,ktop,kcnv,islimsk,
60 ! & dot,ncloud,ud_mf,dd_mf,dt_mf,me)
61 !
62  use machine , only : kind_phys
63  use funcphys , only : fpvs
64  use physcons, grav => con_g, cp => con_cp, hvap => con_hvap
65  &, rv => con_rv, fv => con_fvirt, t0c => con_t0c
66  &, cvap => con_cvap, cliq => con_cliq
67  &, eps => con_eps, epsm1 => con_epsm1
68  implicit none
69 !
70  integer im, ix, km, jcap, ncloud,
71  & kbot(im), ktop(im), kcnv(im)
72 ! &, me
73  real(kind=kind_phys) delt
74  real(kind=kind_phys) psp(im), delp(ix,km), prslp(ix,km)
75  real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km),
76  & ql(ix,km,2),q1(ix,km), t1(ix,km),
77  & u1(ix,km), v1(ix,km), !rcs(im),
78  & cldwrk(im), rn(im),
79  & dot(ix,km), phil(ix,km),
80  & cnvw(ix,km), cnvc(ix,km),
81  & ud_mf(im,km),dd_mf(im,km),dt_mf(im,km) ! hchuang code change mass flux output
82 !
83  integer i, indx, jmn, k, kk, km1
84  integer, dimension(im), intent(in) :: islimsk
85 ! integer latd,lond
86 !
87  real(kind=kind_phys) clam, cxlamu, xlamde, xlamdd
88 !
89 ! real(kind=kind_phys) detad
90  real(kind=kind_phys) adw, aup, aafac,
91  & beta, betal, betas,
92  & c0, dellat, delta,
93  & desdt, dg,
94  & dh, dhh, dp,
95  & dq, dqsdp, dqsdt, dt,
96  & dt2, dtmax, dtmin, dv1h,
97  & dv1q, dv2h, dv2q, dv1u,
98  & dv1v, dv2u, dv2v, dv3q,
99  & dv3h, dv3u, dv3v,
100  & dz, dz1, e1, edtmax,
101  & edtmaxl, edtmaxs, el2orc, elocp,
102  & es, etah, cthk, dthk,
103  & evef, evfact, evfactl, fact1,
104  & fact2, factor, fjcap, fkm,
105  & g, gamma, pprime,
106  & qlk, qrch, qs, c1,
107  & rain, rfact, shear, tem1,
108  & val, val1,
109  & val2, w1, w1l, w1s,
110  & w2, w2l, w2s, w3,
111  & w3l, w3s, w4, w4l,
112  & w4s, xdby, xpw, xpwd,
113  & xqrch, mbdt, tem,
114  & ptem, ptem1, pgcon
115 !
116  integer kb(im), kbcon(im), kbcon1(im),
117  & ktcon(im), ktcon1(im),
118  & jmin(im), lmin(im), kbmax(im),
119  & kbm(im), kmax(im)
120 !
121  real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im),
122  & delhbar(im), delq(im), delq2(im),
123  & delqbar(im), delqev(im), deltbar(im),
124  & deltv(im), dtconv(im), edt(im),
125  & edto(im), edtx(im), fld(im),
126  & hcdo(im,km), hmax(im), hmin(im),
127  & ucdo(im,km), vcdo(im,km),aa2(im),
128  & pbcdif(im), pdot(im), po(im,km),
129  & pwavo(im), pwevo(im), xlamud(im),
130  & qcdo(im,km), qcond(im), qevap(im),
131  & rntot(im), vshear(im), xaa0(im),
132  & xk(im), xlamd(im),
133  & xmb(im), xmbmax(im), xpwav(im),
134  & xpwev(im), delubar(im),delvbar(im)
135 cj
136  real(kind=kind_phys) cincr, cincrmax, cincrmin
137 cj
138 c physical parameters
139  parameter(g=grav)
140  parameter(elocp=hvap/cp,
141  & el2orc=hvap*hvap/(rv*cp))
142  parameter(c0=.002,c1=.002,delta=fv)
143  parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
144  parameter(cthk=150.,cincrmax=180.,cincrmin=120.,dthk=25.)
145 c local variables and arrays
146  real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
147  & uo(im,km), vo(im,km), qeso(im,km)
148 c cloud water
149 ! real(kind=kind_phys) tvo(im,km)
150  real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
151  & dbyo(im,km), zo(im,km), xlamue(im,km),
152  & fent1(im,km), fent2(im,km), frh(im,km),
153  & heo(im,km), heso(im,km),
154  & qrcd(im,km), dellah(im,km), dellaq(im,km),
155  & dellau(im,km), dellav(im,km), hcko(im,km),
156  & ucko(im,km), vcko(im,km), qcko(im,km),
157  & eta(im,km), etad(im,km), zi(im,km),
158  & qrcko(im,km), qrcdo(im,km),
159  & pwo(im,km), pwdo(im,km),
160  & tx1(im), sumx(im), cnvwt(im,km)
161 ! &, rhbar(im)
162 !
163  logical totflg, cnvflg(im), flg(im)
164 !
165  real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
166 ! save pcrit, acritt
167  data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,
168  & 350.,300.,250.,200.,150./
169  data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,
170  & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
171 c gdas derived acrit
172 c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
173 c & .743,.813,.886,.947,1.138,1.377,1.896/
174  real(kind=kind_phys) tf, tcr, tcrf
175  parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
176 !
177 c-----------------------------------------------------------------------
180 !************************************************************************
181 ! convert input pa terms to cb terms -- moorthi
182  ps = psp * 0.001
183  prsl = prslp * 0.001
184  del = delp * 0.001
185 !************************************************************************
186 !
187 !
188  km1 = km - 1
190 c
191 c initialize arrays
192 c
193  do i=1,im
194  cnvflg(i) = .true.
195  rn(i)=0.
196  kbot(i)=km+1
197  ktop(i)=0
198  kbcon(i)=km
199  ktcon(i)=1
200  dtconv(i) = 3600.
201  cldwrk(i) = 0.
202  pdot(i) = 0.
203  pbcdif(i)= 0.
204  lmin(i) = 1
205  jmin(i) = 1
206  qlko_ktcon(i) = 0.
207  edt(i) = 0.
208  edto(i) = 0.
209  edtx(i) = 0.
210  acrt(i) = 0.
211  acrtfct(i) = 1.
212  aa1(i) = 0.
213  aa2(i) = 0.
214  xaa0(i) = 0.
215  pwavo(i)= 0.
216  pwevo(i)= 0.
217  xpwav(i)= 0.
218  xpwev(i)= 0.
219  vshear(i) = 0.
220  enddo
222  do k = 1, km
223  do i = 1, im
224  cnvw(i,k) = 0.
225  cnvc(i,k) = 0.
226  enddo
227  enddo
229 ! hchuang code change
230  do k = 1, km
231  do i = 1, im
232  ud_mf(i,k) = 0.
233  dd_mf(i,k) = 0.
234  dt_mf(i,k) = 0.
235  enddo
236  enddo
238 c
239  do k = 1, 15
240  acrit(k) = acritt(k) * (975. - pcrit(k))
241  enddo
242  dt2 = delt
243  val = 1200.
244  dtmin = max(dt2, val )
245  val = 3600.
246  dtmax = max(dt2, val )
247 c model tunable parameters are all here
248  mbdt = 10.
249  edtmaxl = .3
250  edtmaxs = .3
251  clam = .1
252  aafac = .1
253 ! betal = .15
254 ! betas = .15
255  betal = .05
256  betas = .05
257 c evef = 0.07
258  evfact = 0.3
259  evfactl = 0.3
260 !
261  cxlamu = 1.0e-4
262  xlamde = 1.0e-4
263  xlamdd = 1.0e-4
264 !
265 ! pgcon = 0.7 ! gregory et al. (1997, qjrms)
266  pgcon = 0.55 ! zhang & wu (2003,jas)
267  fjcap = (float(jcap) / 126.) ** 2
268  val = 1.
269  fjcap = max(fjcap,val)
270  fkm = (float(km) / 28.) ** 2
271  fkm = max(fkm,val)
272  w1l = -8.e-3
273  w2l = -4.e-2
274  w3l = -5.e-3
275  w4l = -5.e-4
276  w1s = -2.e-4
277  w2s = -2.e-3
278  w3s = -1.e-3
279  w4s = -2.e-5
281 c
282 c define top layer for search of the downdraft originating layer
283 c and the maximum thetae for updraft
284 c
285  do i=1,im
286  kbmax(i) = km
287  kbm(i) = km
288  kmax(i) = km
289  tx1(i) = 1.0 / ps(i)
290  enddo
291 !
292  do k = 1, km
293  do i=1,im
294  if (prsl(i,k)*tx1(i) .gt. 0.04) kmax(i) = k + 1
295  if (prsl(i,k)*tx1(i) .gt. 0.45) kbmax(i) = k + 1
296  if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
297  enddo
298  enddo
299  do i=1,im
300  kmax(i) = min(km,kmax(i))
301  kbmax(i) = min(kbmax(i),kmax(i))
302  kbm(i) = min(kbm(i),kmax(i))
303  enddo
304 c
305 c hydrostatic height assume zero terr and initially assume
306 c updraft entrainment rate as an inverse function of height
307 c
309  do k = 1, km
310  do i=1,im
311  zo(i,k) = phil(i,k) / g
312  enddo
313  enddo
315  do k = 1, km1
316  do i=1,im
317  zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
318  xlamue(i,k) = clam / zi(i,k)
319  enddo
320  enddo
322 c
323 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324 c convert surface pressure to mb from cb
325 c
326  do k = 1, km
327  do i = 1, im
328  if (k .le. kmax(i)) then
329  pfld(i,k) = prsl(i,k) * 10.0
330  eta(i,k) = 1.
331  fent1(i,k)= 1.
332  fent2(i,k)= 1.
333  frh(i,k) = 0.
334  hcko(i,k) = 0.
335  qcko(i,k) = 0.
336  qrcko(i,k)= 0.
337  ucko(i,k) = 0.
338  vcko(i,k) = 0.
339  etad(i,k) = 1.
340  hcdo(i,k) = 0.
341  qcdo(i,k) = 0.
342  ucdo(i,k) = 0.
343  vcdo(i,k) = 0.
344  qrcd(i,k) = 0.
345  qrcdo(i,k)= 0.
346  dbyo(i,k) = 0.
347  pwo(i,k) = 0.
348  pwdo(i,k) = 0.
349  dellal(i,k) = 0.
350  to(i,k) = t1(i,k)
351  qo(i,k) = q1(i,k)
352  uo(i,k) = u1(i,k)
353  vo(i,k) = v1(i,k)
354 ! uo(i,k) = u1(i,k) * rcs(i)
355 ! vo(i,k) = v1(i,k) * rcs(i)
356  cnvwt(i,k)= 0.
357  endif
358  enddo
359  enddo
360 c
361 c column variables
362 c p is pressure of the layer (mb)
363 c t is temperature at t-dt (k)..tn
364 c q is mixing ratio at t-dt (kg/kg)..qn
365 c to is temperature at t+dt (k)... this is after advection and turbulan
366 c qo is mixing ratio at t+dt (kg/kg)..q1
367 c
369  do k = 1, km
370  do i=1,im
371  if (k .le. kmax(i)) then
372  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
373  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
374  val1 = 1.e-8
375  qeso(i,k) = max(qeso(i,k), val1)
376  val2 = 1.e-10
377  qo(i,k) = max(qo(i,k), val2 )
378 ! qo(i,k) = min(qo(i,k),qeso(i,k))
379 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
380  endif
381  enddo
382  enddo
383 c
384 c compute moist static energy
385 c
387  do k = 1, km
388  do i=1,im
389  if (k .le. kmax(i)) then
390 ! tem = g * zo(i,k) + cp * to(i,k)
391  tem = phil(i,k) + cp * to(i,k)
392  heo(i,k) = tem + hvap * qo(i,k)
393  heso(i,k) = tem + hvap * qeso(i,k)
394 c heo(i,k) = min(heo(i,k),heso(i,k))
395  endif
396  enddo
397  enddo
398 
400 c
401 c determine level with largest moist static energy
402 c this is the level where updraft starts
403 c
405  do i=1,im
406  hmax(i) = heo(i,1)
407  kb(i) = 1
408  enddo
409  do k = 2, km
410  do i=1,im
411  if (k .le. kbm(i)) then
412  if(heo(i,k).gt.hmax(i)) then
413  kb(i) = k
414  hmax(i) = heo(i,k)
415  endif
416  endif
417  enddo
418  enddo
420 c
421  do k = 1, km1
422  do i=1,im
423  if (k .le. kmax(i)-1) then
424  dz = .5 * (zo(i,k+1) - zo(i,k))
425  dp = .5 * (pfld(i,k+1) - pfld(i,k))
426  es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
427  pprime = pfld(i,k+1) + epsm1 * es
428  qs = eps * es / pprime
429  dqsdp = - qs / pprime
430  desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
431  dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
432  gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
433  dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
434  dq = dqsdt * dt + dqsdp * dp
435  to(i,k) = to(i,k+1) + dt
436  qo(i,k) = qo(i,k+1) + dq
437  po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
438  endif
439  enddo
440  enddo
442 !
443  do k = 1, km1
444  do i=1,im
445  if (k .le. kmax(i)-1) then
446  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
447  qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
448  val1 = 1.e-8
449  qeso(i,k) = max(qeso(i,k), val1)
450  val2 = 1.e-10
451  qo(i,k) = max(qo(i,k), val2 )
452 ! qo(i,k) = min(qo(i,k),qeso(i,k))
453  frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), 1.)
454  heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
455  & cp * to(i,k) + hvap * qo(i,k)
456  heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
457  & cp * to(i,k) + hvap * qeso(i,k)
458  uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
459  vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
460  endif
461  enddo
462  enddo
463 c
464 c look for the level of free convection as cloud base
465 c
467  do i=1,im
468  flg(i) = .true.
469  kbcon(i) = kmax(i)
470  enddo
471  do k = 1, km1
472  do i=1,im
473  if (flg(i).and.k.le.kbmax(i)) then
474  if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
475  kbcon(i) = k
476  flg(i) = .false.
477  endif
478  endif
479  enddo
480  enddo
482 c
483  do i=1,im
484  if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
485  enddo
486 !!
487  totflg = .true.
488  do i=1,im
489  totflg = totflg .and. (.not. cnvflg(i))
490  enddo
491  if(totflg) return
492 !!
493 c
494 c determine critical convective inhibition
495 c as a function of vertical velocity at cloud base.
496 c
498  do i=1,im
499  if(cnvflg(i)) then
500 ! pdot(i) = 10.* dot(i,kbcon(i))
501  pdot(i) = 0.01 * dot(i,kbcon(i)) ! now dot is in pa/s
502  endif
503  enddo
504  do i=1,im
505  if(cnvflg(i)) then
506  if(islimsk(i) == 1) then
507  w1 = w1l
508  w2 = w2l
509  w3 = w3l
510  w4 = w4l
511  else
512  w1 = w1s
513  w2 = w2s
514  w3 = w3s
515  w4 = w4s
516  endif
517  if(pdot(i).le.w4) then
518  tem = (pdot(i) - w4) / (w3 - w4)
519  elseif(pdot(i).ge.-w4) then
520  tem = - (pdot(i) + w4) / (w4 - w3)
521  else
522  tem = 0.
523  endif
524  val1 = -1.
525  tem = max(tem,val1)
526  val2 = 1.
527  tem = min(tem,val2)
528  tem = 1. - tem
529  tem1= .5*(cincrmax-cincrmin)
530  cincr = cincrmax - tem * tem1
531  pbcdif(i) = pfld(i,kb(i)) - pfld(i,kbcon(i))
532  if(pbcdif(i).gt.cincr) then
533  cnvflg(i) = .false.
534  endif
535  endif
536  enddo
537 !!
538  totflg = .true.
539  do i=1,im
540  totflg = totflg .and. (.not. cnvflg(i))
541  enddo
542  if(totflg) return
543 !!
544 c
545 c assume that updraft entrainment rate above cloud base is
546 c same as that at cloud base
547 c
553  do k = 2, km1
554  do i=1,im
555  if(cnvflg(i).and.
556  & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
557  xlamue(i,k) = xlamue(i,kbcon(i))
558  endif
559  enddo
560  enddo
561 c
562 c assume the detrainment rate for the updrafts to be same as
563 c the entrainment rate at cloud base
564 c
566  do i = 1, im
567  if(cnvflg(i)) then
568  xlamud(i) = xlamue(i,kbcon(i))
569  endif
570  enddo
571 c
572 c functions rapidly decreasing with height, mimicking a cloud ensemble
573 c (bechtold et al., 2008)
574 c
575  do k = 2, km1
576  do i=1,im
577  if(cnvflg(i).and.
578  & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
579  tem = qeso(i,k)/qeso(i,kbcon(i))
580  fent1(i,k) = tem**2
581  fent2(i,k) = tem**3
582  endif
583  enddo
584  enddo
585 c
586 c final entrainment rate as the sum of turbulent part and organized entrainment
587 c depending on the environmental relative humidity
588 c (bechtold et al., 2008)
589 c
590  do k = 2, km1
591  do i=1,im
592  if(cnvflg(i).and.
593  & (k.ge.kbcon(i).and.k.lt.kmax(i))) then
594  tem = cxlamu * frh(i,k) * fent2(i,k)
595  xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
596  endif
597  enddo
598  enddo
599 c
600 c determine updraft mass flux for the subcloud layers
601 c
607  do k = km1, 1, -1
608  do i = 1, im
609  if (cnvflg(i)) then
610  if(k.lt.kbcon(i).and.k.ge.kb(i)) then
611  dz = zi(i,k+1) - zi(i,k)
612  ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
613  eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
614  endif
615  endif
616  enddo
617  enddo
618 c
619 c compute mass flux above cloud base
620 c
621  do k = 2, km1
622  do i = 1, im
623  if(cnvflg(i))then
624  if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
625  dz = zi(i,k) - zi(i,k-1)
626  ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
627  eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
628  endif
629  endif
630  enddo
631  enddo
632 c
633 c compute updraft cloud properties
634 c
636  do i = 1, im
637  if(cnvflg(i)) then
638  indx = kb(i)
639  hcko(i,indx) = heo(i,indx)
640  ucko(i,indx) = uo(i,indx)
641  vcko(i,indx) = vo(i,indx)
642  pwavo(i) = 0.
643  endif
644  enddo
645 c
646 c cloud property is modified by the entrainment process
647 c
649  do k = 2, km1
650  do i = 1, im
651  if (cnvflg(i)) then
652  if(k.gt.kb(i).and.k.lt.kmax(i)) then
653  dz = zi(i,k) - zi(i,k-1)
654  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
655  tem1 = 0.5 * xlamud(i) * dz
656  factor = 1. + tem - tem1
657  ptem = 0.5 * tem + pgcon
658  ptem1= 0.5 * tem - pgcon
659  hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
660  & (heo(i,k)+heo(i,k-1)))/factor
661  ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)
662  & +ptem1*uo(i,k-1))/factor
663  vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)
664  & +ptem1*vo(i,k-1))/factor
665  dbyo(i,k) = hcko(i,k) - heso(i,k)
666  endif
667  endif
668  enddo
669  enddo
670 c
671 c taking account into convection inhibition due to existence of
672 c dry layers below cloud base
673 c
675  do i=1,im
676  flg(i) = cnvflg(i)
677  kbcon1(i) = kmax(i)
678  enddo
679  do k = 2, km1
680  do i=1,im
681  if (flg(i).and.k.lt.kmax(i)) then
682  if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
683  kbcon1(i) = k
684  flg(i) = .false.
685  endif
686  endif
687  enddo
688  enddo
689  do i=1,im
690  if(cnvflg(i)) then
691  if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
692  endif
693  enddo
694  do i=1,im
695  if(cnvflg(i)) then
696  tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
697  if(tem.gt.dthk) then
698  cnvflg(i) = .false.
699  endif
700  endif
701  enddo
702 !!
703  totflg = .true.
704  do i = 1, im
705  totflg = totflg .and. (.not. cnvflg(i))
706  enddo
707  if(totflg) return
708 !!
709 c
710 c determine first guess cloud top as the level of zero buoyancy
711 c
713  do i = 1, im
714  flg(i) = cnvflg(i)
715  ktcon(i) = 1
716  enddo
717  do k = 2, km1
718  do i = 1, im
719  if (flg(i).and.k .lt. kmax(i)) then
720  if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
721  ktcon(i) = k
722  flg(i) = .false.
723  endif
724  endif
725  enddo
726  enddo
727 c
728  do i = 1, im
729  if(cnvflg(i)) then
730  tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
731  if(tem.lt.cthk) cnvflg(i) = .false.
732  endif
733  enddo
734 !!
735  totflg = .true.
736  do i = 1, im
737  totflg = totflg .and. (.not. cnvflg(i))
738  enddo
739  if(totflg) return
740 !!
741 c
742 c search for downdraft originating level above theta-e minimum
743 c
745  do i = 1, im
746  if(cnvflg(i)) then
747  hmin(i) = heo(i,kbcon1(i))
748  lmin(i) = kbmax(i)
749  jmin(i) = kbmax(i)
750  endif
751  enddo
752  do k = 2, km1
753  do i = 1, im
754  if (cnvflg(i) .and. k .le. kbmax(i)) then
755  if(k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
756  lmin(i) = k + 1
757  hmin(i) = heo(i,k)
758  endif
759  endif
760  enddo
761  enddo
762 c
763 c make sure that jmin(i) is within the cloud
764 c
765  do i = 1, im
766  if(cnvflg(i)) then
767  jmin(i) = min(lmin(i),ktcon(i)-1)
768  jmin(i) = max(jmin(i),kbcon1(i)+1)
769  if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
770  endif
771  enddo
772 c
773 c specify upper limit of mass flux at cloud base
774 c
776  do i = 1, im
777  if(cnvflg(i)) then
778 ! xmbmax(i) = .1
779 !
780  k = kbcon(i)
781  dp = 1000. * del(i,k)
782  xmbmax(i) = dp / (g * dt2)
783 !
784 ! tem = dp / (g * dt2)
785 ! xmbmax(i) = min(tem, xmbmax(i))
786  endif
787  enddo
788 c
789 c compute cloud moisture property and precipitation
790 c
792  do i = 1, im
793  if (cnvflg(i)) then
794  aa1(i) = 0.
795  qcko(i,kb(i)) = qo(i,kb(i))
796  qrcko(i,kb(i)) = qo(i,kb(i))
797 ! rhbar(i) = 0.
798  endif
799  enddo
801  do k = 2, km1
802  do i = 1, im
803  if (cnvflg(i)) then
804  if(k.gt.kb(i).and.k.lt.ktcon(i)) then
805  dz = zi(i,k) - zi(i,k-1)
806  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
807  qrch = qeso(i,k)
808  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
809 cj
810  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
811  tem1 = 0.5 * xlamud(i) * dz
812  factor = 1. + tem - tem1
813  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
814  & (qo(i,k)+qo(i,k-1)))/factor
815  qrcko(i,k) = qcko(i,k)
816 cj
817  dq = eta(i,k) * (qcko(i,k) - qrch)
818 c
819 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
820 c
821 c check if there is excess moisture to release latent heat
822 c
823  if(k.ge.kbcon(i).and.dq.gt.0.) then
824  etah = .5 * (eta(i,k) + eta(i,k-1))
825  if(ncloud.gt.0..and.k.gt.jmin(i)) then
826  dp = 1000. * del(i,k)
827  qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
828  dellal(i,k) = etah * c1 * dz * qlk * g / dp
829  else
830  qlk = dq / (eta(i,k) + etah * c0 * dz)
831  endif
832  aa1(i) = aa1(i) - dz * g * qlk
833  qcko(i,k) = qlk + qrch
834  pwo(i,k) = etah * c0 * dz * qlk
835  pwavo(i) = pwavo(i) + pwo(i,k)
836 ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
837  cnvwt(i,k) = etah * qlk * g / dp
838  endif
839  endif
840  endif
841  enddo
842  enddo
843 c
844 ! do i = 1, im
845 ! if(cnvflg(i)) then
846 ! indx = ktcon(i) - kb(i) - 1
847 ! rhbar(i) = rhbar(i) / float(indx)
848 ! endif
849 ! enddo
850 c
851 c calculate cloud work function
852 c
858  do k = 2, km1
859  do i = 1, im
860  if (cnvflg(i)) then
861  if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
862  dz1 = zo(i,k+1) - zo(i,k)
863  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
864  rfact = 1. + delta * cp * gamma
865  & * to(i,k) / hvap
866  aa1(i) = aa1(i) +
867  & dz1 * (g / (cp * to(i,k)))
868  & * dbyo(i,k) / (1. + gamma)
869  & * rfact
870  val = 0.
871  aa1(i)=aa1(i)+
872  & dz1 * g * delta *
873  & max(val,(qeso(i,k) - qo(i,k)))
874  endif
875  endif
876  enddo
877  enddo
879  do i = 1, im
880  if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
881  enddo
882 !!
883  totflg = .true.
884  do i=1,im
885  totflg = totflg .and. (.not. cnvflg(i))
886  enddo
887  if(totflg) return
888 !!
889 c
890 c estimate the onvective overshooting as the level
891 c where the [aafac * cloud work function] becomes zero,
892 c which is the final cloud top
893 c
895  do i = 1, im
896  if (cnvflg(i)) then
897  aa2(i) = aafac * aa1(i)
898  endif
899  enddo
900 c
901  do i = 1, im
902  flg(i) = cnvflg(i)
903  ktcon1(i) = kmax(i) - 1
904  enddo
905  do k = 2, km1
906  do i = 1, im
907  if (flg(i)) then
908  if(k.ge.ktcon(i).and.k.lt.kmax(i)) then
909  dz1 = zo(i,k+1) - zo(i,k)
910  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
911  rfact = 1. + delta * cp * gamma
912  & * to(i,k) / hvap
913  aa2(i) = aa2(i) +
914  & dz1 * (g / (cp * to(i,k)))
915  & * dbyo(i,k) / (1. + gamma)
916  & * rfact
917  if(aa2(i).lt.0.) then
918  ktcon1(i) = k
919  flg(i) = .false.
920  endif
921  endif
922  endif
923  enddo
924  enddo
925 c
926 c compute cloud moisture property, detraining cloud water
927 c and precipitation in overshooting layers
928 c
930  do k = 2, km1
931  do i = 1, im
932  if (cnvflg(i)) then
933  if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
934  dz = zi(i,k) - zi(i,k-1)
935  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
936  qrch = qeso(i,k)
937  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
938 cj
939  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
940  tem1 = 0.5 * xlamud(i) * dz
941  factor = 1. + tem - tem1
942  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
943  & (qo(i,k)+qo(i,k-1)))/factor
944  qrcko(i,k) = qcko(i,k)
945 cj
946  dq = eta(i,k) * (qcko(i,k) - qrch)
947 c
948 c check if there is excess moisture to release latent heat
949 c
950  if(dq.gt.0.) then
951  etah = .5 * (eta(i,k) + eta(i,k-1))
952  if(ncloud.gt.0.) then
953  dp = 1000. * del(i,k)
954  qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
955  dellal(i,k) = etah * c1 * dz * qlk * g / dp
956  else
957  qlk = dq / (eta(i,k) + etah * c0 * dz)
958  endif
959  qcko(i,k) = qlk + qrch
960  pwo(i,k) = etah * c0 * dz * qlk
961  pwavo(i) = pwavo(i) + pwo(i,k)
962 ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
963  cnvwt(i,k) = etah * qlk * g / dp
964  endif
965  endif
966  endif
967  enddo
968  enddo
969 c
970 c exchange ktcon with ktcon1
971 c
973  do i = 1, im
974  if(cnvflg(i)) then
975  kk = ktcon(i)
976  ktcon(i) = ktcon1(i)
977  ktcon1(i) = kk
978  endif
979  enddo
980 c
981 c this section is ready for cloud water
982 c
984  if(ncloud.gt.0) then
985 c
986 c compute liquid and vapor separation at cloud top
987 c
988  do i = 1, im
989  if(cnvflg(i)) then
990  k = ktcon(i) - 1
991  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
992  qrch = qeso(i,k)
993  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
994  dq = qcko(i,k) - qrch
995 c
996 c check if there is excess moisture to release latent heat
997 c
998  if(dq.gt.0.) then
999  qlko_ktcon(i) = dq
1000  qcko(i,k) = qrch
1001  endif
1002  endif
1003  enddo
1004  endif
1005 c
1006 ccccc if(lat.eq.latd.and.lon.eq.lond.and.cnvflg(i)) then
1007 ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
1008 ccccc endif
1009 c
1010 c------- downdraft calculations
1011 c
1012 c--- compute precipitation efficiency in terms of windshear
1013 c
1020  do i = 1, im
1021  if(cnvflg(i)) then
1022  vshear(i) = 0.
1023  endif
1024  enddo
1025  do k = 2, km
1026  do i = 1, im
1027  if (cnvflg(i)) then
1028  if(k.gt.kb(i).and.k.le.ktcon(i)) then
1029  shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
1030  & + (vo(i,k)-vo(i,k-1)) ** 2)
1031  vshear(i) = vshear(i) + shear
1032  endif
1033  endif
1034  enddo
1035  enddo
1036  do i = 1, im
1037  if(cnvflg(i)) then
1038  vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
1039  e1=1.591-.639*vshear(i)
1040  & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1041  edt(i)=1.-e1
1042  val = .9
1043  edt(i) = min(edt(i),val)
1044  val = .0
1045  edt(i) = max(edt(i),val)
1046  edto(i)=edt(i)
1047  edtx(i)=edt(i)
1048  endif
1049  enddo
1050 c
1051 c determine detrainment rate between 1 and kbcon
1052 c
1058  do i = 1, im
1059  if(cnvflg(i)) then
1060  sumx(i) = 0.
1061  endif
1062  enddo
1063  do k = 1, km1
1064  do i = 1, im
1065  if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
1066  dz = zi(i,k+1) - zi(i,k)
1067  sumx(i) = sumx(i) + dz
1068  endif
1069  enddo
1070  enddo
1071  do i = 1, im
1072  beta = betas
1073  if(islimsk(i) == 1) beta = betal
1074  if(cnvflg(i)) then
1075  dz = (sumx(i)+zi(i,1))/float(kbcon(i))
1076  tem = 1./float(kbcon(i))
1077  xlamd(i) = (1.-beta**tem)/dz
1078  endif
1079  enddo
1080 c
1081 c determine downdraft mass flux
1082 c
1084  do k = km1, 1, -1
1085  do i = 1, im
1086  if (cnvflg(i) .and. k .le. kmax(i)-1) then
1087  if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
1088  dz = zi(i,k+1) - zi(i,k)
1089  ptem = xlamdd - xlamde
1090  etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1091  else if(k.lt.kbcon(i)) then
1092  dz = zi(i,k+1) - zi(i,k)
1093  ptem = xlamd(i) + xlamdd - xlamde
1094  etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1095  endif
1096  endif
1097  enddo
1098  enddo
1099 c
1100 c--- downdraft moisture properties
1101 c
1103  do i = 1, im
1104  if(cnvflg(i)) then
1105  jmn = jmin(i)
1106  hcdo(i,jmn) = heo(i,jmn)
1107  qcdo(i,jmn) = qo(i,jmn)
1108  qrcdo(i,jmn)= qo(i,jmn)
1109  ucdo(i,jmn) = uo(i,jmn)
1110  vcdo(i,jmn) = vo(i,jmn)
1111  pwevo(i) = 0.
1112  endif
1113  enddo
1114 cj
1116  do k = km1, 1, -1
1117  do i = 1, im
1118  if (cnvflg(i) .and. k.lt.jmin(i)) then
1119  dz = zi(i,k+1) - zi(i,k)
1120  if(k.ge.kbcon(i)) then
1121  tem = xlamde * dz
1122  tem1 = 0.5 * xlamdd * dz
1123  else
1124  tem = xlamde * dz
1125  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1126  endif
1127  factor = 1. + tem - tem1
1128  ptem = 0.5 * tem - pgcon
1129  ptem1= 0.5 * tem + pgcon
1130  hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
1131  & (heo(i,k)+heo(i,k+1)))/factor
1132  ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1)
1133  & +ptem1*uo(i,k))/factor
1134  vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1)
1135  & +ptem1*vo(i,k))/factor
1136  dbyo(i,k) = hcdo(i,k) - heso(i,k)
1137  endif
1138  enddo
1139  enddo
1140 c
1142  do k = km1, 1, -1
1143  do i = 1, im
1144  if (cnvflg(i).and.k.lt.jmin(i)) then
1145  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1146  qrcdo(i,k) = qeso(i,k)+
1147  & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
1148 ! detad = etad(i,k+1) - etad(i,k)
1149 cj
1150  dz = zi(i,k+1) - zi(i,k)
1151  if(k.ge.kbcon(i)) then
1152  tem = xlamde * dz
1153  tem1 = 0.5 * xlamdd * dz
1154  else
1155  tem = xlamde * dz
1156  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1157  endif
1158  factor = 1. + tem - tem1
1159  qcdo(i,k) = ((1.-tem1)*qrcdo(i,k+1)+tem*0.5*
1160  & (qo(i,k)+qo(i,k+1)))/factor
1161 cj
1162 ! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
1163 ! & etad(i,k) * qrcdo(i,k)
1164 ! pwdo(i,k) = pwdo(i,k) - detad *
1165 ! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
1166 cj
1167  pwdo(i,k) = etad(i,k) * (qcdo(i,k) - qrcdo(i,k))
1168  pwevo(i) = pwevo(i) + pwdo(i,k)
1169  endif
1170  enddo
1171  enddo
1172 c
1173 c--- final downdraft strength dependent on precip
1174 c--- efficiency (edt), normalized condensate (pwav), and
1175 c--- evaporate (pwev)
1176 c
1178  do i = 1, im
1179  edtmax = edtmaxl
1180  if(islimsk(i) == 0) edtmax = edtmaxs
1181  if(cnvflg(i)) then
1182  if(pwevo(i).lt.0.) then
1183  edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1184  edto(i) = min(edto(i),edtmax)
1185  else
1186  edto(i) = 0.
1187  endif
1188  endif
1189  enddo
1190 c
1191 c--- downdraft cloudwork functions
1192 c
1194  do k = km1, 1, -1
1195  do i = 1, im
1196  if (cnvflg(i) .and. k .lt. jmin(i)) then
1197  gamma = el2orc * qeso(i,k) / to(i,k)**2
1198  dhh=hcdo(i,k)
1199  dt=to(i,k)
1200  dg=gamma
1201  dh=heso(i,k)
1202  dz=-1.*(zo(i,k+1)-zo(i,k))
1203  aa1(i)=aa1(i)+edto(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg))
1204  & *(1.+delta*cp*dg*dt/hvap)
1205  val=0.
1206  aa1(i)=aa1(i)+edto(i)*
1207  & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
1208  endif
1209  enddo
1210  enddo
1212  do i = 1, im
1213  if(cnvflg(i).and.aa1(i).le.0.) then
1214  cnvflg(i) = .false.
1215  endif
1216  enddo
1217 !!
1218  totflg = .true.
1219  do i=1,im
1220  totflg = totflg .and. (.not. cnvflg(i))
1221  enddo
1222  if(totflg) return
1223 !!
1224 c
1225 c--- what would the change be, that a cloud with unit mass
1226 c--- will do to the environment?
1227 c
1229  do k = 1, km
1230  do i = 1, im
1231  if(cnvflg(i) .and. k .le. kmax(i)) then
1232  dellah(i,k) = 0.
1233  dellaq(i,k) = 0.
1234  dellau(i,k) = 0.
1235  dellav(i,k) = 0.
1236  endif
1237  enddo
1238  enddo
1239  do i = 1, im
1240  if(cnvflg(i)) then
1241  dp = 1000. * del(i,1)
1242  dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1)
1243  & - heo(i,1)) * g / dp
1244  dellaq(i,1) = edto(i) * etad(i,1) * (qrcdo(i,1)
1245  & - qo(i,1)) * g / dp
1246  dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1)
1247  & - uo(i,1)) * g / dp
1248  dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1)
1249  & - vo(i,1)) * g / dp
1250  endif
1251  enddo
1252 c
1253 c--- changed due to subsidence and entrainment
1254 c
1255  do k = 2, km1
1256  do i = 1, im
1257  if (cnvflg(i).and.k.lt.ktcon(i)) then
1258  aup = 1.
1259  if(k.le.kb(i)) aup = 0.
1260  adw = 1.
1261  if(k.gt.jmin(i)) adw = 0.
1262  dp = 1000. * del(i,k)
1263  dz = zi(i,k) - zi(i,k-1)
1264 c
1265  dv1h = heo(i,k)
1266  dv2h = .5 * (heo(i,k) + heo(i,k-1))
1267  dv3h = heo(i,k-1)
1268  dv1q = qo(i,k)
1269  dv2q = .5 * (qo(i,k) + qo(i,k-1))
1270  dv3q = qo(i,k-1)
1271  dv1u = uo(i,k)
1272  dv2u = .5 * (uo(i,k) + uo(i,k-1))
1273  dv3u = uo(i,k-1)
1274  dv1v = vo(i,k)
1275  dv2v = .5 * (vo(i,k) + vo(i,k-1))
1276  dv3v = vo(i,k-1)
1277 c
1278  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1279  tem1 = xlamud(i)
1280 c
1281  if(k.le.kbcon(i)) then
1282  ptem = xlamde
1283  ptem1 = xlamd(i)+xlamdd
1284  else
1285  ptem = xlamde
1286  ptem1 = xlamdd
1287  endif
1288 cj
1289  dellah(i,k) = dellah(i,k) +
1290  & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h
1291  & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h
1292  & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz
1293  & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1294  & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz
1295  & ) *g/dp
1296 cj
1297  dellaq(i,k) = dellaq(i,k) +
1298  & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q
1299  & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q
1300  & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz
1301  & + aup*tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz
1302  & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qcdo(i,k-1))*dz
1303  & ) *g/dp
1304 cj
1305  dellau(i,k) = dellau(i,k) +
1306  & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1u
1307  & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3u
1308  & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz
1309  & + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz
1310  & + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz
1311  & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u)
1312  & ) *g/dp
1313 cj
1314  dellav(i,k) = dellav(i,k) +
1315  & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1v
1316  & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3v
1317  & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz
1318  & + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz
1319  & + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz
1320  & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v)
1321  & ) *g/dp
1322 cj
1323  endif
1324  enddo
1325  enddo
1326 c
1327 c------- cloud top
1328 c
1329  do i = 1, im
1330  if(cnvflg(i)) then
1331  indx = ktcon(i)
1332  dp = 1000. * del(i,indx)
1333  dv1h = heo(i,indx-1)
1334  dellah(i,indx) = eta(i,indx-1) *
1335  & (hcko(i,indx-1) - dv1h) * g / dp
1336  dv1q = qo(i,indx-1)
1337  dellaq(i,indx) = eta(i,indx-1) *
1338  & (qcko(i,indx-1) - dv1q) * g / dp
1339  dv1u = uo(i,indx-1)
1340  dellau(i,indx) = eta(i,indx-1) *
1341  & (ucko(i,indx-1) - dv1u) * g / dp
1342  dv1v = vo(i,indx-1)
1343  dellav(i,indx) = eta(i,indx-1) *
1344  & (vcko(i,indx-1) - dv1v) * g / dp
1345 c
1346 c cloud water
1347 c
1348  dellal(i,indx) = eta(i,indx-1) *
1349  & qlko_ktcon(i) * g / dp
1350  endif
1351  enddo
1352 c
1353 c------- final changed variable per unit mass flux
1354 c
1356  do k = 1, km
1357  do i = 1, im
1358  if (cnvflg(i).and.k .le. kmax(i)) then
1359  if(k.gt.ktcon(i)) then
1360  qo(i,k) = q1(i,k)
1361  to(i,k) = t1(i,k)
1362  endif
1363  if(k.le.ktcon(i)) then
1364  qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
1365  dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1366  to(i,k) = dellat * mbdt + t1(i,k)
1367  val = 1.e-10
1368  qo(i,k) = max(qo(i,k), val )
1369  endif
1370  endif
1371  enddo
1372  enddo
1373 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1374 c
1375 c--- the above changed environment is now used to calulate the
1376 c--- effect the arbitrary cloud (with unit mass flux)
1377 c--- would have on the stability,
1378 c--- which then is used to calculate the real mass flux,
1379 c--- necessary to keep this change in balance with the large-scale
1380 c--- destabilization.
1381 c
1382 c--- environmental conditions again, first heights
1383 c
1387  do k = 1, km
1388  do i = 1, im
1389  if(cnvflg(i) .and. k .le. kmax(i)) then
1390  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1391  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
1392  val = 1.e-8
1393  qeso(i,k) = max(qeso(i,k), val )
1394 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
1395  endif
1396  enddo
1397  enddo
1398 c
1399 c--- moist static energy
1400 c
1401 !! - Recalculate moist static energy and saturation moist static energy.
1402  do k = 1, km1
1403  do i = 1, im
1404  if(cnvflg(i) .and. k .le. kmax(i)-1) then
1405  dz = .5 * (zo(i,k+1) - zo(i,k))
1406  dp = .5 * (pfld(i,k+1) - pfld(i,k))
1407  es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
1408  pprime = pfld(i,k+1) + epsm1 * es
1409  qs = eps * es / pprime
1410  dqsdp = - qs / pprime
1411  desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1412  dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
1413  gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1414  dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
1415  dq = dqsdt * dt + dqsdp * dp
1416  to(i,k) = to(i,k+1) + dt
1417  qo(i,k) = qo(i,k+1) + dq
1418  po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
1419  endif
1420  enddo
1421  enddo
1422  do k = 1, km1
1423  do i = 1, im
1424  if(cnvflg(i) .and. k .le. kmax(i)-1) then
1425  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1426  qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
1427  val1 = 1.e-8
1428  qeso(i,k) = max(qeso(i,k), val1)
1429  val2 = 1.e-10
1430  qo(i,k) = max(qo(i,k), val2 )
1431 ! qo(i,k) = min(qo(i,k),qeso(i,k))
1432  heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
1433  & cp * to(i,k) + hvap * qo(i,k)
1434  heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
1435  & cp * to(i,k) + hvap * qeso(i,k)
1436  endif
1437  enddo
1438  enddo
1439  do i = 1, im
1440  if(cnvflg(i)) then
1441  k = kmax(i)
1442  heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
1443  heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
1444 c heo(i,k) = min(heo(i,k),heso(i,k))
1445  endif
1446  enddo
1447 c
1448 c**************************** static control
1449 c
1450 c------- moisture and cloud work functions
1451 c
1453  do i = 1, im
1454  if(cnvflg(i)) then
1455  xaa0(i) = 0.
1456  xpwav(i) = 0.
1457  endif
1458  enddo
1459 c
1460  do i = 1, im
1461  if(cnvflg(i)) then
1462  indx = kb(i)
1463  hcko(i,indx) = heo(i,indx)
1464  qcko(i,indx) = qo(i,indx)
1465  endif
1466  enddo
1467  do k = 2, km1
1468  do i = 1, im
1469  if (cnvflg(i)) then
1470  if(k.gt.kb(i).and.k.le.ktcon(i)) then
1471  dz = zi(i,k) - zi(i,k-1)
1472  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1473  tem1 = 0.5 * xlamud(i) * dz
1474  factor = 1. + tem - tem1
1475  hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
1476  & (heo(i,k)+heo(i,k-1)))/factor
1477  endif
1478  endif
1479  enddo
1480  enddo
1481  do k = 2, km1
1482  do i = 1, im
1483  if (cnvflg(i)) then
1484  if(k.gt.kb(i).and.k.lt.ktcon(i)) then
1485  dz = zi(i,k) - zi(i,k-1)
1486  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1487  xdby = hcko(i,k) - heso(i,k)
1488  xqrch = qeso(i,k)
1489  & + gamma * xdby / (hvap * (1. + gamma))
1490 cj
1491  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1492  tem1 = 0.5 * xlamud(i) * dz
1493  factor = 1. + tem - tem1
1494  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
1495  & (qo(i,k)+qo(i,k-1)))/factor
1496 cj
1497  dq = eta(i,k) * (qcko(i,k) - xqrch)
1498 c
1499  if(k.ge.kbcon(i).and.dq.gt.0.) then
1500  etah = .5 * (eta(i,k) + eta(i,k-1))
1501  if(ncloud.gt.0..and.k.gt.jmin(i)) then
1502  qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
1503  else
1504  qlk = dq / (eta(i,k) + etah * c0 * dz)
1505  endif
1506  if(k.lt.ktcon1(i)) then
1507  xaa0(i) = xaa0(i) - dz * g * qlk
1508  endif
1509  qcko(i,k) = qlk + xqrch
1510  xpw = etah * c0 * dz * qlk
1511  xpwav(i) = xpwav(i) + xpw
1512  endif
1513  endif
1514  if(k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
1515  dz1 = zo(i,k+1) - zo(i,k)
1516  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1517  rfact = 1. + delta * cp * gamma
1518  & * to(i,k) / hvap
1519  xaa0(i) = xaa0(i)
1520  & + dz1 * (g / (cp * to(i,k)))
1521  & * xdby / (1. + gamma)
1522  & * rfact
1523  val=0.
1524  xaa0(i)=xaa0(i)+
1525  & dz1 * g * delta *
1526  & max(val,(qeso(i,k) - qo(i,k)))
1527  endif
1528  endif
1529  enddo
1530  enddo
1531 c
1532 c------- downdraft calculations
1533 c
1534 c--- downdraft moisture properties
1535 c
1537  do i = 1, im
1538  if(cnvflg(i)) then
1539  jmn = jmin(i)
1540  hcdo(i,jmn) = heo(i,jmn)
1541  qcdo(i,jmn) = qo(i,jmn)
1542  qrcd(i,jmn) = qo(i,jmn)
1543  xpwev(i) = 0.
1544  endif
1545  enddo
1546 cj
1547  do k = km1, 1, -1
1548  do i = 1, im
1549  if (cnvflg(i) .and. k.lt.jmin(i)) then
1550  dz = zi(i,k+1) - zi(i,k)
1551  if(k.ge.kbcon(i)) then
1552  tem = xlamde * dz
1553  tem1 = 0.5 * xlamdd * dz
1554  else
1555  tem = xlamde * dz
1556  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1557  endif
1558  factor = 1. + tem - tem1
1559  hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
1560  & (heo(i,k)+heo(i,k+1)))/factor
1561  endif
1562  enddo
1563  enddo
1564 cj
1565  do k = km1, 1, -1
1566  do i = 1, im
1567  if (cnvflg(i) .and. k .lt. jmin(i)) then
1568  dq = qeso(i,k)
1569  dt = to(i,k)
1570  gamma = el2orc * dq / dt**2
1571  dh = hcdo(i,k) - heso(i,k)
1572  qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
1573 ! detad = etad(i,k+1) - etad(i,k)
1574 cj
1575  dz = zi(i,k+1) - zi(i,k)
1576  if(k.ge.kbcon(i)) then
1577  tem = xlamde * dz
1578  tem1 = 0.5 * xlamdd * dz
1579  else
1580  tem = xlamde * dz
1581  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1582  endif
1583  factor = 1. + tem - tem1
1584  qcdo(i,k) = ((1.-tem1)*qrcd(i,k+1)+tem*0.5*
1585  & (qo(i,k)+qo(i,k+1)))/factor
1586 cj
1587 ! xpwd = etad(i,k+1) * qcdo(i,k+1) -
1588 ! & etad(i,k) * qrcd(i,k)
1589 ! xpwd = xpwd - detad *
1590 ! & .5 * (qrcd(i,k) + qrcd(i,k+1))
1591 cj
1592  xpwd = etad(i,k) * (qcdo(i,k) - qrcd(i,k))
1593  xpwev(i) = xpwev(i) + xpwd
1594  endif
1595  enddo
1596  enddo
1597 c
1598  do i = 1, im
1599  edtmax = edtmaxl
1600  if(islimsk(i) == 0) edtmax = edtmaxs
1601  if(cnvflg(i)) then
1602  if(xpwev(i).ge.0.) then
1603  edtx(i) = 0.
1604  else
1605  edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1606  edtx(i) = min(edtx(i),edtmax)
1607  endif
1608  endif
1609  enddo
1610 c
1611 c
1612 c--- downdraft cloudwork functions
1613 c
1614 c
1615  do k = km1, 1, -1
1616  do i = 1, im
1617  if (cnvflg(i) .and. k.lt.jmin(i)) then
1618  gamma = el2orc * qeso(i,k) / to(i,k)**2
1619  dhh=hcdo(i,k)
1620  dt= to(i,k)
1621  dg= gamma
1622  dh= heso(i,k)
1623  dz=-1.*(zo(i,k+1)-zo(i,k))
1624  xaa0(i)=xaa0(i)+edtx(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg))
1625  & *(1.+delta*cp*dg*dt/hvap)
1626  val=0.
1627  xaa0(i)=xaa0(i)+edtx(i)*
1628  & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
1629  endif
1630  enddo
1631  enddo
1632 c
1633 c calculate critical cloud work function
1634 c
1637  do i = 1, im
1638  if(cnvflg(i)) then
1639  if(pfld(i,ktcon(i)).lt.pcrit(15))then
1640  acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i)))
1641  & /(975.-pcrit(15))
1642  else if(pfld(i,ktcon(i)).gt.pcrit(1))then
1643  acrt(i)=acrit(1)
1644  else
1645  k = int((850. - pfld(i,ktcon(i)))/50.) + 2
1646  k = min(k,15)
1647  k = max(k,2)
1648  acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*
1649  & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
1650  endif
1651  endif
1652  enddo
1654  do i = 1, im
1655  if(cnvflg(i)) then
1656  if(islimsk(i) == 1) then
1657  w1 = w1l
1658  w2 = w2l
1659  w3 = w3l
1660  w4 = w4l
1661  else
1662  w1 = w1s
1663  w2 = w2s
1664  w3 = w3s
1665  w4 = w4s
1666  endif
1667 c
1668 c modify critical cloud workfunction by cloud base vertical velocity
1669 c
1670  if(pdot(i).le.w4) then
1671  acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
1672  elseif(pdot(i).ge.-w4) then
1673  acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
1674  else
1675  acrtfct(i) = 0.
1676  endif
1677  val1 = -1.
1678  acrtfct(i) = max(acrtfct(i),val1)
1679  val2 = 1.
1680  acrtfct(i) = min(acrtfct(i),val2)
1681  acrtfct(i) = 1. - acrtfct(i)
1682 c
1683 c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
1684 c
1685 c if(rhbar(i).ge..8) then
1686 c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
1687 c endif
1688 c
1689 c modify adjustment time scale by cloud base vertical velocity
1690 c
1692  dtconv(i) = dt2 + max((1800. - dt2),0.) *
1693  & (pdot(i) - w2) / (w1 - w2)
1694 c dtconv(i) = max(dtconv(i), dt2)
1695 c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
1696  dtconv(i) = max(dtconv(i),dtmin)
1697  dtconv(i) = min(dtconv(i),dtmax)
1698 c
1699  endif
1700  enddo
1701 c
1702 c--- large scale forcing
1703 c
1709  do i= 1, im
1710  if(cnvflg(i)) then
1711  fld(i)=(aa1(i)-acrt(i)* acrtfct(i))/dtconv(i)
1712  if(fld(i).le.0.) cnvflg(i) = .false.
1713  endif
1719  if(cnvflg(i)) then
1720 c xaa0(i) = max(xaa0(i),0.)
1721  xk(i) = (xaa0(i) - aa1(i)) / mbdt
1722  if(xk(i).ge.0.) cnvflg(i) = .false.
1723  endif
1724 c
1725 c--- kernel, cloud base mass flux
1726 c
1731  if(cnvflg(i)) then
1732  xmb(i) = -fld(i) / xk(i)
1733  xmb(i) = min(xmb(i),xmbmax(i))
1734  endif
1735  enddo
1736 !!
1738  totflg = .true.
1739  do i=1,im
1740  totflg = totflg .and. (.not. cnvflg(i))
1741  enddo
1742  if(totflg) return
1743 !!
1744 c
1745 c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
1746 c
1747 
1748  do k = 1, km
1749  do i = 1, im
1750  if (cnvflg(i) .and. k .le. kmax(i)) then
1751  to(i,k) = t1(i,k)
1752  qo(i,k) = q1(i,k)
1753  uo(i,k) = u1(i,k)
1754  vo(i,k) = v1(i,k)
1755  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1756  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
1757  val = 1.e-8
1758  qeso(i,k) = max(qeso(i,k), val )
1759  endif
1760  enddo
1761  enddo
1762 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1763 c
1764 c--- feedback: simply the changes from the cloud with unit mass flux
1765 c--- multiplied by the mass flux necessary to keep the
1766 c--- equilibrium with the larger-scale.
1767 c
1772  do i = 1, im
1773  delhbar(i) = 0.
1774  delqbar(i) = 0.
1775  deltbar(i) = 0.
1776  delubar(i) = 0.
1777  delvbar(i) = 0.
1778  qcond(i) = 0.
1779  enddo
1780  do k = 1, km
1781  do i = 1, im
1782  if (cnvflg(i) .and. k .le. kmax(i)) then
1783  if(k.le.ktcon(i)) then
1784  dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1785  t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1786  q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1787 ! tem = 1./rcs(i)
1788 ! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1789 ! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1790  u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
1791  v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
1792  dp = 1000. * del(i,k)
1793  delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
1794  delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
1795  deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
1796  delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
1797  delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
1798  endif
1799  endif
1800  enddo
1801  enddo
1803  do k = 1, km
1804  do i = 1, im
1805  if (cnvflg(i) .and. k .le. kmax(i)) then
1806  if(k.le.ktcon(i)) then
1807  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1808  qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
1809  val = 1.e-8
1810  qeso(i,k) = max(qeso(i,k), val )
1811  endif
1812  endif
1813  enddo
1814  enddo
1815 c
1817  do i = 1, im
1818  rntot(i) = 0.
1819  delqev(i) = 0.
1820  delq2(i) = 0.
1821  flg(i) = cnvflg(i)
1822  enddo
1823  do k = km, 1, -1
1824  do i = 1, im
1825  if (cnvflg(i) .and. k .le. kmax(i)) then
1826  if(k.lt.ktcon(i)) then
1827  aup = 1.
1828  if(k.le.kb(i)) aup = 0.
1829  adw = 1.
1830  if(k.ge.jmin(i)) adw = 0.
1831  rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
1832  rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
1833  endif
1834  endif
1835  enddo
1836  enddo
1840  do k = km, 1, -1
1841  do i = 1, im
1842  if (k .le. kmax(i)) then
1843  deltv(i) = 0.
1844  delq(i) = 0.
1845  qevap(i) = 0.
1846  if(cnvflg(i).and.k.lt.ktcon(i)) then
1847  aup = 1.
1848  if(k.le.kb(i)) aup = 0.
1849  adw = 1.
1850  if(k.ge.jmin(i)) adw = 0.
1851  rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
1852  rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
1853  endif
1854  if(flg(i).and.k.lt.ktcon(i)) then
1855  evef = edt(i) * evfact
1856  if(islimsk(i) == 1) evef=edt(i) * evfactl
1857 ! if(islimsk(i) == 1) evef=.07
1858 c if(islimsk(i) == 1) evef = 0.
1859  qcond(i) = evef * (q1(i,k) - qeso(i,k))
1860  & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
1861  dp = 1000. * del(i,k)
1862  if(rn(i).gt.0..and.qcond(i).lt.0.) then
1863  qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
1864  qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
1865  delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
1866  endif
1867  if(rn(i).gt.0..and.qcond(i).lt.0..and.
1868  & delq2(i).gt.rntot(i)) then
1869  qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
1870  flg(i) = .false.
1871  endif
1872  if(rn(i).gt.0..and.qevap(i).gt.0.) then
1873  q1(i,k) = q1(i,k) + qevap(i)
1874  t1(i,k) = t1(i,k) - elocp * qevap(i)
1875  rn(i) = rn(i) - .001 * qevap(i) * dp / g
1876  deltv(i) = - elocp*qevap(i)/dt2
1877  delq(i) = + qevap(i)/dt2
1878  delqev(i) = delqev(i) + .001*dp*qevap(i)/g
1879  endif
1880  dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
1881  delqbar(i) = delqbar(i) + delq(i)*dp/g
1882  deltbar(i) = deltbar(i) + deltv(i)*dp/g
1883  endif
1884  endif
1885  enddo
1886  enddo
1887 cj
1888 ! do i = 1, im
1889 ! if(me.eq.31.and.cnvflg(i)) then
1890 ! if(cnvflg(i)) then
1891 ! print *, ' deep delhbar, delqbar, deltbar = ',
1892 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
1893 ! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
1894 ! print *, ' precip =', hvap*rn(i)*1000./dt2
1895 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
1896 ! endif
1897 ! enddo
1898 c
1899 c precipitation rate converted to actual precip
1900 c in unit of m instead of kg
1901 c
1902  do i = 1, im
1903  if(cnvflg(i)) then
1904 c
1905 c in the event of upper level rain evaporation and lower level downdraft
1906 c moistening, rn can become negative, in this case, we back out of the
1907 c heating and the moistening
1908 c
1909 
1910  if(rn(i).lt.0..and..not.flg(i)) rn(i) = 0.
1911  if(rn(i).le.0.) then
1912  rn(i) = 0.
1913  else
1914  ktop(i) = ktcon(i)
1915  kbot(i) = kbcon(i)
1916  kcnv(i) = 1
1917  cldwrk(i) = aa1(i)
1918  endif
1919  endif
1920  enddo
1921 c
1922 c convective cloud water
1923 c
1925  do k = 1, km
1926  do i = 1, im
1927  if (cnvflg(i) .and. rn(i).gt.0.) then
1928  if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1929  cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
1930  endif
1931  endif
1932  enddo
1933  enddo
1934 c
1935 c convective cloud cover
1936 c
1938  do k = 1, km
1939  do i = 1, im
1940  if (cnvflg(i) .and. rn(i).gt.0.) then
1941  if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1942  cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
1943  cnvc(i,k) = min(cnvc(i,k), 0.6)
1944  cnvc(i,k) = max(cnvc(i,k), 0.0)
1945  endif
1946  endif
1947  enddo
1948  enddo
1949 
1950 c
1951 c cloud water
1952 c
1954  if (ncloud.gt.0) then
1955 !
1956  do k = 1, km
1957  do i = 1, im
1958  if (cnvflg(i) .and. rn(i).gt.0.) then
1959  if (k.gt.kb(i).and.k.le.ktcon(i)) then
1960  tem = dellal(i,k) * xmb(i) * dt2
1961  tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
1962  if (ql(i,k,2) .gt. -999.0) then
1963  ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
1964  ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
1965  else
1966  ql(i,k,1) = ql(i,k,1) + tem
1967  endif
1968  endif
1969  endif
1970  enddo
1971  enddo
1972 !
1973  endif
1974 c
1976  do k = 1, km
1977  do i = 1, im
1978  if(cnvflg(i).and.rn(i).le.0.) then
1979  if (k .le. kmax(i)) then
1980  t1(i,k) = to(i,k)
1981  q1(i,k) = qo(i,k)
1982  u1(i,k) = uo(i,k)
1983  v1(i,k) = vo(i,k)
1984  endif
1985  endif
1986  enddo
1987  enddo
1988 !
1989 ! hchuang code change
1990 !
1992  do k = 1, km
1993  do i = 1, im
1994  if(cnvflg(i).and.rn(i).gt.0.) then
1995  if(k.ge.kb(i) .and. k.lt.ktop(i)) then
1996  ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
1997  endif
1998  endif
1999  enddo
2000  enddo
2002  do i = 1, im
2003  if(cnvflg(i).and.rn(i).gt.0.) then
2004  k = ktop(i)-1
2005  dt_mf(i,k) = ud_mf(i,k)
2006  endif
2007  enddo
2009  do k = 1, km
2010  do i = 1, im
2011  if(cnvflg(i).and.rn(i).gt.0.) then
2012  if(k.ge.1 .and. k.le.jmin(i)) then
2013  dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
2014  endif
2015  endif
2016  enddo
2017  enddo
2018 !!
2019  return
2022  end
2023 ! \section original Original Documentation
2024 ! Penetrative convection is simulated following Pan and Wu (1994), which is based on Arakawa and Schubert(1974) as simplified by Grell (1993) and with a saturated downdraft. Convection occurs when the cloud work function (CWF) exceeds a certain threshold. Mass flux of the cloud is determined using a quasi-equilibrium assumption based on this threshold CWF. The CWF is a function of temperature and moisture in each air column of the model gridpoint. The temperature and moisture profiles are adjusted towards the equilibrium CWF within a specified time scale using the deduced mass flux. A major simplification of the original Arakawa-Shubert scheme is to consider only the deepest cloud and not the spectrum of clouds. The cloud model incorporates a downdraft mechanism as well as the evaporation of precipitation. Entrainment of the updraft and detrainment of the downdraft in the sub-cloud layers are included. Downdraft strength is based on the vertical wind shear through the cloud. The critical CWF is a function of the cloud base vertical motion. As the large-scale rising motion becomes strong, the CWF [similar to convective available potential energy (CAPE)] is allowed to approach zero (therefore approaching neutral stability).
2025 !
2026 ! Mass fluxes induced in the updraft and the downdraft are allowed to transport momentum. The momentum exchange is calculated through the mass flux formulation in a manner similar to that for heat and moisture. The effect of the convection-induced pressure gradient force on cumulus momentum transport is parameterized in terms of mass flux and vertical wind shear (Han and Pan, 2006). As a result, the cumulus momentum exchange is reduced by about 55 % compared to the full exchange.
2027 !
2028 ! The entrainment rate in cloud layers is dependent upon environmental humidity (Han and Pan, 2010). A drier environment increases the entrainment, suppressing the convection. The entrainment rate in sub-cloud layers is given as inversely proportional to height. The detrainment rate is assumed to be a constant in all layers and equal to the entrainment rate value at cloud base, which is O(10-4). The liquid water in the updraft layer is assumed to be detrained from the layers above the level of the minimum moist static energy into the grid-scale cloud water with conversion parameter of 0.002 m-1, which is same as the rain conversion parameter.
2029 !
2030 ! Following Han and Pan (2010), the trigger condition is that a parcel lifted from the convection starting level without entrainment must reach its level of free convection within 120-180 hPa of ascent, proportional to the large-scale vertical velocity. This is intended to produce more convection in large-scale convergent regions but less convection in large-scale subsidence regions. Another important trigger mechanism is to include the effect of environmental humidity in the sub-cloud layer, taking into account convection inhibition due to existence of dry layers below cloud base. On the other hand, the cloud parcel might overshoot beyond the level of neutral buoyancy due to its inertia, eventually stopping its overshoot at cloud top. The CWF is used to model the overshoot. The overshoot of the cloud top is stopped at the height where a parcel lifted from the neutral buoyancy level with energy equal to 10% of the CWF would first have zero energy.
2031 !
2032 ! Deep convection parameterization (SAS) modifications include:
2033 ! - Detraining cloud water from every updraft layer
2034 ! - Starting convection from the level of maximum moist static energy within PBL
2035 ! - Random cloud top is eliminated and only deepest cloud is considered
2036 ! - Cloud water is detrained from every cloud layer
2037 ! - Finite entrainment and detrainment rates for heat, moisture, and momentum are specified
2038 ! - Similar to shallow convection scheme,
2039 ! - entrainment rate is given to be inversely proportional to height in sub-cloud layers
2040 ! - detrainment rate is set to be a constant as entrainment rate at the cloud base.
2041 ! -Above cloud base, an organized entrainment is added, which is a function of environmental relative humidity
physcons::con_cliq
real(kind=kind_phys), parameter con_cliq
spec heat H2O liq ( )
Definition: physcons.f:86
physcons::con_g
real(kind=kind_phys), parameter con_g
gravity ( )
Definition: physcons.f:59
physcons::con_rv
real(kind=kind_phys), parameter con_rv
gas constant H2O ( )
Definition: physcons.f:78
physcons::con_t0c
real(kind=kind_phys), parameter con_t0c
temp at 0C (K)
Definition: physcons.f:96
sascnvn
subroutine sascnvn(im, ix, km, jcap, delt, delp, prslp, psp, phil, ql, q1, t1, u1, v1, cldwrk, rn, kbot, ktop, kcnv, islimsk, dot, ncloud, ud_mf, dd_mf, dt_mf, cnvw, cnvc)
This subroutine contains the entirety of the SAS deep convection scheme.
Definition: sascnvn.f:59
physcons::con_cvap
real(kind=kind_phys), parameter con_cvap
spec heat H2O gas ( )
Definition: physcons.f:84
physcons::con_cp
real(kind=kind_phys), parameter con_cp
spec heat air at p ( )
Definition: physcons.f:80
physcons::con_hvap
real(kind=kind_phys), parameter con_hvap
lat heat H2O cond ( )
Definition: physcons.f:90