GFS Physics Documentation
samfdeepcnv.f
Go to the documentation of this file.
1 
13 
16 
75  subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql,
76  & q1,t1,u1,v1,cldwrk,rn,kbot,ktop,kcnv,islimsk,garea,
77  & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc,
78  & clam,c0s,c1,betal,betas,evfact,evfactl,pgcon,asolfac)
79 !
80  use machine , only : kind_phys
81  use funcphys , only : fpvs
82  use physcons, grav => con_g, cp => con_cp, hvap => con_hvap
83  &, rv => con_rv, fv => con_fvirt, t0c => con_t0c
84  &, rd => con_rd, cvap => con_cvap, cliq => con_cliq
85  &, eps => con_eps, epsm1 => con_epsm1
86  implicit none
87 !
88  integer, intent(in) :: im, ix, km, ncloud
89  integer, intent(in) :: islimsk(im)
90  real(kind=kind_phys), intent(in) :: delt
91  real(kind=kind_phys), intent(in) :: psp(im), delp(ix,km),
92  & prslp(ix,km), garea(im), dot(ix,km), phil(ix,km)
93 
94  integer, intent(inout) :: kcnv(im)
95  real(kind=kind_phys), intent(inout) :: ql(ix,km,2),
96  & q1(ix,km), t1(ix,km), u1(ix,km), v1(ix,km)
97 
98  integer, intent(out) :: kbot(im), ktop(im)
99  real(kind=kind_phys), intent(out) :: cldwrk(im),
100  & rn(im), cnvw(ix,km), cnvc(ix,km),
101  & ud_mf(im,km),dd_mf(im,km), dt_mf(im,km)
102 !
103 !------local variables
104  integer i, indx, jmn, k, kk, km1, n
105 ! integer latd,lond
106 !
107  real(kind=kind_phys) clam, cxlamu, cxlamd,
108  & xlamde, xlamdd,
109  & crtlamu, crtlamd
110 !
111 ! real(kind=kind_phys) detad
112  real(kind=kind_phys) adw, aup, aafac,
113  & beta, betal, betas,
114  & c0l, c0s, d0,
115  & c1, asolfac,
116  & dellat, delta, desdt, dg,
117  & dh, dhh, dp,
118  & dq, dqsdp, dqsdt, dt,
119  & dt2, dtmax, dtmin,
120  & dxcrtas, dxcrtuf,
121  & dv1h, dv2h, dv3h,
122  & dv1q, dv2q, dv3q,
123  & dz, dz1, e1, edtmax,
124  & edtmaxl, edtmaxs, el2orc, elocp,
125  & es, etah,
126  & cthk, dthk,
127  & evef, evfact, evfactl, fact1,
128  & fact2, factor,
129  & g, gamma, pprime, cm,
130  & qlk, qrch, qs,
131  & rain, rfact, shear, tfac,
132  & val, val1, val2,
133  & w1, w1l, w1s, w2,
134  & w2l, w2s, w3, w3l,
135  & w3s, w4, w4l, w4s,
136  & rho, betaw,
137  & xdby, xpw, xpwd,
138 ! & xqrch, mbdt, tem,
139  & xqrch, tem, tem1, tem2,
140  & ptem, ptem1, ptem2,
141  & pgcon
142 !
143  integer kb(im), kbcon(im), kbcon1(im),
144  & ktcon(im), ktcon1(im), ktconn(im),
145  & jmin(im), lmin(im), kbmax(im),
146  & kbm(im), kmax(im)
147 !
148 ! real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im),
149  real(kind=kind_phys) aa1(im),
150  & ps(im), del(ix,km), prsl(ix,km),
151  & umean(im), tauadv(im), gdx(im),
152  & delhbar(im), delq(im), delq2(im),
153  & delqbar(im), delqev(im), deltbar(im),
154  & deltv(im), dtconv(im), edt(im),
155  & edto(im), edtx(im), fld(im),
156  & hcdo(im,km), hmax(im), hmin(im),
157  & ucdo(im,km), vcdo(im,km),aa2(im),
158  & pdot(im), po(im,km),
159  & pwavo(im), pwevo(im), mbdt(im),
160  & qcdo(im,km), qcond(im), qevap(im),
161  & rntot(im), vshear(im), xaa0(im),
162  & xk(im), xlamd(im), cina(im),
163  & xmb(im), xmbmax(im), xpwav(im),
164  & xpwev(im), xlamx(im),
165  & delubar(im),delvbar(im)
166 !
167  real(kind=kind_phys) c0(im)
168 cj
169  real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn,
170  & cinacr, cinacrmx, cinacrmn
171 cj
172 !
173 ! parameters for updraft velocity calculation
174  real(kind=kind_phys) bet1, cd1, f1, gam1,
175  & bb1, bb2
176 ! & bb1, bb2, wucb
177 !
178 c physical parameters
179 ! parameter(g=grav,asolfac=0.958)
180  parameter(g=grav)
181  parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
182 ! parameter(c0s=.002,c1=.002,d0=.01)
183  parameter(d0=.01)
184 ! parameter(c0l=c0s*asolfac)
185 !
186 ! asolfac: aerosol-aware parameter based on Lim (2011)
187 ! asolfac= cx / c0s(=.002)
188 ! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s)
189 ! Nccn: CCN number concentration in cm^(-3)
190 ! Until a realistic Nccn is provided, Nccns are assumed
191 ! as Nccn=100 for sea and Nccn=1000 for land
192 !
193  parameter(cm=1.0,delta=fv)
194  parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
195  parameter(cthk=200.,dthk=25.)
196  parameter(cinpcrmx=180.,cinpcrmn=120.)
197 ! parameter(cinacrmx=-120.,cinacrmn=-120.)
198  parameter(cinacrmx=-120.,cinacrmn=-80.)
199  parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
200  parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3)
201 !
202 ! local variables and arrays
203  real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
204  & uo(im,km), vo(im,km), qeso(im,km)
205 ! for updraft velocity calculation
206  real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km)
207  real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im)
208 !
209 c cloud water
210 ! real(kind=kind_phys) tvo(im,km)
211  real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
212  & dbyo(im,km), zo(im,km),
213  & xlamue(im,km), xlamud(im,km),
214  & fent1(im,km), fent2(im,km), frh(im,km),
215  & heo(im,km), heso(im,km),
216  & qrcd(im,km), dellah(im,km), dellaq(im,km),
217  & dellau(im,km), dellav(im,km), hcko(im,km),
218  & ucko(im,km), vcko(im,km), qcko(im,km),
219  & eta(im,km), etad(im,km), zi(im,km),
220  & qrcko(im,km), qrcdo(im,km),
221  & pwo(im,km), pwdo(im,km), c0t(im,km),
222  & tx1(im), sumx(im), cnvwt(im,km)
223 ! &, rhbar(im)
224 !
225  logical totflg, cnvflg(im), asqecflg(im), flg(im)
226 !
227 ! asqecflg: flag for the quasi-equilibrium assumption of Arakawa-Schubert
228 !
229 ! real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
230 !! save pcrit, acritt
231 ! data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,
232 ! & 350.,300.,250.,200.,150./
233 ! data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,
234 ! & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
235 c gdas derived acrit
236 c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
237 c & .743,.813,.886,.947,1.138,1.377,1.896/
238  real(kind=kind_phys) tf, tcr, tcrf
239  parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
240 !
241 c-----------------------------------------------------------------------
244 !************************************************************************
245 ! convert input Pa terms to Cb terms -- Moorthi
246  ps = psp * 0.001
247  prsl = prslp * 0.001
248  del = delp * 0.001
249 !************************************************************************
250 !
251 !
252  km1 = km - 1
254 c
255 c initialize arrays
256 c
257  do i=1,im
258  cnvflg(i) = .true.
259  rn(i)=0.
260  mbdt(i)=10.
261  kbot(i)=km+1
262  ktop(i)=0
263  kbcon(i)=km
264  ktcon(i)=1
265  ktconn(i)=1
266  dtconv(i) = 3600.
267  cldwrk(i) = 0.
268  pdot(i) = 0.
269  lmin(i) = 1
270  jmin(i) = 1
271  qlko_ktcon(i) = 0.
272  edt(i) = 0.
273  edto(i) = 0.
274  edtx(i) = 0.
275 ! acrt(i) = 0.
276 ! acrtfct(i) = 1.
277  aa1(i) = 0.
278  aa2(i) = 0.
279  xaa0(i) = 0.
280  cina(i) = 0.
281  pwavo(i)= 0.
282  pwevo(i)= 0.
283  xpwav(i)= 0.
284  xpwev(i)= 0.
285  vshear(i) = 0.
286  gdx(i) = sqrt(garea(i))
287  enddo
288 !
290  do i=1,im
291  if(islimsk(i) == 1) then
292  c0(i) = c0s*asolfac
293  else
294  c0(i) = c0s
295  endif
296  enddo
298  do k = 1, km
299  do i = 1, im
300  if(t1(i,k) > 273.16) then
301  c0t(i,k) = c0(i)
302  else
303  tem = d0 * (t1(i,k) - 273.16)
304  tem1 = exp(tem)
305  c0t(i,k) = c0(i) * tem1
306  endif
307  enddo
308  enddo
310  do k = 1, km
311  do i = 1, im
312  cnvw(i,k) = 0.
313  cnvc(i,k) = 0.
314  enddo
315  enddo
316 ! hchuang code change
318  do k = 1, km
319  do i = 1, im
320  ud_mf(i,k) = 0.
321  dd_mf(i,k) = 0.
322  dt_mf(i,k) = 0.
323  enddo
324  enddo
325 c
326 ! do k = 1, 15
327 ! acrit(k) = acritt(k) * (975. - pcrit(k))
328 ! enddo
329 !
330  dt2 = delt
331 ! val = 1200.
332  val = 600.
333  dtmin = max(dt2, val )
334 ! val = 5400.
335  val = 10800.
336  dtmax = max(dt2, val )
337 c model tunable parameters are all here
338  edtmaxl = .3
339  edtmaxs = .3
340 ! clam = .1
341  aafac = .1
342 ! betal = .15
343 ! betas = .15
344 ! betal = .05
345 ! betas = .05
346 c evef = 0.07
347 ! evfact = 0.3
348 ! evfactl = 0.3
349 !
350  crtlamu = 1.0e-4
351  crtlamd = 1.0e-4
352 !
353  cxlamu = 1.0e-3
354  cxlamd = 1.0e-4
355  xlamde = 1.0e-4
356  xlamdd = 1.0e-4
357 !
358 ! pgcon = 0.7 ! Gregory et al. (1997, QJRMS)
359 ! pgcon = 0.55 ! Zhang & Wu (2003,JAS)
360 !
361  w1l = -8.e-3
362  w2l = -4.e-2
363  w3l = -5.e-3
364  w4l = -5.e-4
365  w1s = -2.e-4
366  w2s = -2.e-3
367  w3s = -1.e-3
368  w4s = -2.e-5
369 c
370 c define top layer for search of the downdraft originating layer
371 c and the maximum thetae for updraft
372 c
374  do i=1,im
375  kbmax(i) = km
376  kbm(i) = km
377  kmax(i) = km
378  tx1(i) = 1.0 / ps(i)
379  enddo
380 !
381  do k = 1, km
382  do i=1,im
383  if (prsl(i,k)*tx1(i) > 0.04) kmax(i) = k + 1
384  if (prsl(i,k)*tx1(i) > 0.45) kbmax(i) = k + 1
385  if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1
386  enddo
387  enddo
388  do i=1,im
389  kmax(i) = min(km,kmax(i))
390  kbmax(i) = min(kbmax(i),kmax(i))
391  kbm(i) = min(kbm(i),kmax(i))
392  enddo
393 c
394 c hydrostatic height assume zero terr and initially assume
395 c updraft entrainment rate as an inverse function of height
396 c
398  do k = 1, km
399  do i=1,im
400  zo(i,k) = phil(i,k) / g
401  enddo
402  enddo
404  do k = 1, km1
405  do i=1,im
406  zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
407  xlamue(i,k) = clam / zi(i,k)
408 ! xlamue(i,k) = max(xlamue(i,k), crtlamu)
409  enddo
410  enddo
411 c
412 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
413 c convert surface pressure to mb from cb
414 c
416  do k = 1, km
417  do i = 1, im
418  if (k <= kmax(i)) then
419  pfld(i,k) = prsl(i,k) * 10.0
420  eta(i,k) = 1.
421  fent1(i,k)= 1.
422  fent2(i,k)= 1.
423  frh(i,k) = 0.
424  hcko(i,k) = 0.
425  qcko(i,k) = 0.
426  qrcko(i,k)= 0.
427  ucko(i,k) = 0.
428  vcko(i,k) = 0.
429  etad(i,k) = 1.
430  hcdo(i,k) = 0.
431  qcdo(i,k) = 0.
432  ucdo(i,k) = 0.
433  vcdo(i,k) = 0.
434  qrcd(i,k) = 0.
435  qrcdo(i,k)= 0.
436  dbyo(i,k) = 0.
437  pwo(i,k) = 0.
438  pwdo(i,k) = 0.
439  dellal(i,k) = 0.
440  to(i,k) = t1(i,k)
441  qo(i,k) = q1(i,k)
442  uo(i,k) = u1(i,k)
443  vo(i,k) = v1(i,k)
444 ! uo(i,k) = u1(i,k) * rcs(i)
445 ! vo(i,k) = v1(i,k) * rcs(i)
446  wu2(i,k) = 0.
447  buo(i,k) = 0.
448  drag(i,k) = 0.
449  cnvwt(i,k)= 0.
450  endif
451  enddo
452  enddo
454  do k = 1, km
455  do i=1,im
456  if (k <= kmax(i)) then
457  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
458  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
459  val1 = 1.e-8
460  qeso(i,k) = max(qeso(i,k), val1)
461  val2 = 1.e-10
462  qo(i,k) = max(qo(i,k), val2 )
463 ! qo(i,k) = min(qo(i,k),qeso(i,k))
464 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
465  endif
466  enddo
467  enddo
468 c
469 c compute moist static energy
470 c
472  do k = 1, km
473  do i=1,im
474  if (k <= kmax(i)) then
475 ! tem = g * zo(i,k) + cp * to(i,k)
476  tem = phil(i,k) + cp * to(i,k)
477  heo(i,k) = tem + hvap * qo(i,k)
478  heso(i,k) = tem + hvap * qeso(i,k)
479 c heo(i,k) = min(heo(i,k),heso(i,k))
480  endif
481  enddo
482  enddo
483 c
484 c determine level with largest moist static energy
485 c this is the level where updraft starts
486 c
489  do i=1,im
490  hmax(i) = heo(i,1)
491  kb(i) = 1
492  enddo
493  do k = 2, km
494  do i=1,im
495  if (k <= kbm(i)) then
496  if(heo(i,k) > hmax(i)) then
497  kb(i) = k
498  hmax(i) = heo(i,k)
499  endif
500  endif
501  enddo
502  enddo
503 c
505  do k = 1, km1
506  do i=1,im
507  if (k <= kmax(i)-1) then
508  dz = .5 * (zo(i,k+1) - zo(i,k))
509  dp = .5 * (pfld(i,k+1) - pfld(i,k))
510  es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
511  pprime = pfld(i,k+1) + epsm1 * es
512  qs = eps * es / pprime
513  dqsdp = - qs / pprime
514  desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
515  dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
516  gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
517  dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
518  dq = dqsdt * dt + dqsdp * dp
519  to(i,k) = to(i,k+1) + dt
520  qo(i,k) = qo(i,k+1) + dq
521  po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
522  endif
523  enddo
524  enddo
525 !
527  do k = 1, km1
528  do i=1,im
529  if (k <= kmax(i)-1) then
530  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
531  qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
532  val1 = 1.e-8
533  qeso(i,k) = max(qeso(i,k), val1)
534  val2 = 1.e-10
535  qo(i,k) = max(qo(i,k), val2 )
536 ! qo(i,k) = min(qo(i,k),qeso(i,k))
537  frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), 1.)
538  heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
539  & cp * to(i,k) + hvap * qo(i,k)
540  heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
541  & cp * to(i,k) + hvap * qeso(i,k)
542  uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
543  vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
544  endif
545  enddo
546  enddo
547 c
548 c look for the level of free convection as cloud base
549 c
551  do i=1,im
552  flg(i) = .true.
553  kbcon(i) = kmax(i)
554  enddo
555  do k = 1, km1
556  do i=1,im
557  if (flg(i) .and. k <= kbmax(i)) then
558  if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
559  kbcon(i) = k
560  flg(i) = .false.
561  endif
562  endif
563  enddo
564  enddo
565 c
567  do i=1,im
568  if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
569  enddo
570 !!
571  totflg = .true.
572  do i=1,im
573  totflg = totflg .and. (.not. cnvflg(i))
574  enddo
575  if(totflg) return
576 !!
578  do i=1,im
579  if(cnvflg(i)) then
580 ! pdot(i) = 10.* dot(i,kbcon(i))
581  pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
582  endif
583  enddo
584 c
585 c turn off convection if pressure depth between parcel source level
586 c and cloud base is larger than a critical value, cinpcr
587 c
588  do i=1,im
589  if(cnvflg(i)) then
590  if(islimsk(i) == 1) then
591  w1 = w1l
592  w2 = w2l
593  w3 = w3l
594  w4 = w4l
595  else
596  w1 = w1s
597  w2 = w2s
598  w3 = w3s
599  w4 = w4s
600  endif
601  if(pdot(i) <= w4) then
602  tem = (pdot(i) - w4) / (w3 - w4)
603  elseif(pdot(i) >= -w4) then
604  tem = - (pdot(i) + w4) / (w4 - w3)
605  else
606  tem = 0.
607  endif
608  val1 = -1.
609  tem = max(tem,val1)
610  val2 = 1.
611  tem = min(tem,val2)
612  ptem = 1. - tem
613  ptem1= .5*(cinpcrmx-cinpcrmn)
614  cinpcr = cinpcrmx - ptem * ptem1
615  tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
616  if(tem1 > cinpcr) then
617  cnvflg(i) = .false.
618  endif
619  endif
620  enddo
621 !!
622  totflg = .true.
623  do i=1,im
624  totflg = totflg .and. (.not. cnvflg(i))
625  enddo
626  if(totflg) return
627 !!
628 c
629 c assume that updraft entrainment rate above cloud base is
630 c same as that at cloud base
631 c
637  do i=1,im
638  if(cnvflg(i)) then
639  xlamx(i) = xlamue(i,kbcon(i))
640  endif
641  enddo
642  do k = 2, km1
643  do i=1,im
644  if(cnvflg(i).and.
645  & (k > kbcon(i) .and. k < kmax(i))) then
646  xlamue(i,k) = xlamx(i)
647  endif
648  enddo
649  enddo
650 c
651 c specify detrainment rate for the updrafts
652 c
654  do k = 1, km1
655  do i=1,im
656  if(cnvflg(i) .and. k < kmax(i)) then
657  xlamud(i,k) = xlamx(i)
658 ! xlamud(i,k) = crtlamd
659  endif
660  enddo
661  enddo
662 c
663 c functions rapidly decreasing with height, mimicking a cloud ensemble
664 c(bechtold et al., 2008)
665 c
666  do k = 2, km1
667  do i=1,im
668  if(cnvflg(i).and.
669  & (k > kbcon(i) .and. k < kmax(i))) then
670  tem = qeso(i,k)/qeso(i,kbcon(i))
671  fent1(i,k) = tem**2
672  fent2(i,k) = tem**3
673  endif
674  enddo
675  enddo
676 c
677 c final entrainment and detrainment rates as the sum of turbulent part and
678 c organized entrainment depending on the environmental relative humidity
679 c(bechtold et al., 2008)
680 c
681  do k = 2, km1
682  do i=1,im
683  if(cnvflg(i) .and.
684  & (k > kbcon(i) .and. k < kmax(i))) then
685  tem = cxlamu * frh(i,k) * fent2(i,k)
686  xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
687 ! tem1 = cxlamd * frh(i,k)
688 ! xlamud(i,k) = xlamud(i,k) + tem1
689  endif
690  enddo
691  enddo
692 !
693 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
694 c
695 c determine updraft mass flux for the subcloud layers
696 c
702  do k = km1, 1, -1
703  do i = 1, im
704  if (cnvflg(i)) then
705  if(k < kbcon(i) .and. k >= kb(i)) then
706  dz = zi(i,k+1) - zi(i,k)
707  tem = 0.5*(xlamud(i,k)+xlamud(i,k+1))
708  ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-tem
709  eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
710  endif
711  endif
712  enddo
713  enddo
714 c
715 c compute mass flux above cloud base
716 c
717  do i = 1, im
718  flg(i) = cnvflg(i)
719  enddo
720  do k = 2, km1
721  do i = 1, im
722  if(flg(i))then
723  if(k > kbcon(i) .and. k < kmax(i)) then
724  dz = zi(i,k) - zi(i,k-1)
725  tem = 0.5*(xlamud(i,k)+xlamud(i,k-1))
726  ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-tem
727  eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
728  if(eta(i,k) <= 0.) then
729  kmax(i) = k
730  ktconn(i) = k
731  flg(i) = .false.
732  endif
733  endif
734  endif
735  enddo
736  enddo
737 c
738 c compute updraft cloud properties
739 c
741  do i = 1, im
742  if(cnvflg(i)) then
743  indx = kb(i)
744  hcko(i,indx) = heo(i,indx)
745  ucko(i,indx) = uo(i,indx)
746  vcko(i,indx) = vo(i,indx)
747  pwavo(i) = 0.
748  endif
749  enddo
750 c
751 c cloud property is modified by the entrainment process
752 c
753 ! cm is an enhancement factor in entrainment rates for momentum
754 !
756  do k = 2, km1
757  do i = 1, im
758  if (cnvflg(i)) then
759  if(k > kb(i) .and. k < kmax(i)) then
760  dz = zi(i,k) - zi(i,k-1)
761  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
762  tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
763  factor = 1. + tem - tem1
764  hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
765  & (heo(i,k)+heo(i,k-1)))/factor
766  dbyo(i,k) = hcko(i,k) - heso(i,k)
767 !
768  tem = 0.5 * cm * tem
769  factor = 1. + tem
770  ptem = tem + pgcon
771  ptem1= tem - pgcon
772  ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*uo(i,k)
773  & +ptem1*uo(i,k-1))/factor
774  vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*vo(i,k)
775  & +ptem1*vo(i,k-1))/factor
776  endif
777  endif
778  enddo
779  enddo
780 c
781 c taking account into convection inhibition due to existence of
782 c dry layers below cloud base
783 c
785  do i=1,im
786  flg(i) = cnvflg(i)
787  kbcon1(i) = kmax(i)
788  enddo
789  do k = 2, km1
790  do i=1,im
791  if (flg(i) .and. k < kmax(i)) then
792  if(k >= kbcon(i) .and. dbyo(i,k) > 0.) then
793  kbcon1(i) = k
794  flg(i) = .false.
795  endif
796  endif
797  enddo
798  enddo
799  do i=1,im
800  if(cnvflg(i)) then
801  if(kbcon1(i) == kmax(i)) cnvflg(i) = .false.
802  endif
803  enddo
804  do i=1,im
805  if(cnvflg(i)) then
806  tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
807  if(tem > dthk) then
808  cnvflg(i) = .false.
809  endif
810  endif
811  enddo
812 !!
813  totflg = .true.
814  do i = 1, im
815  totflg = totflg .and. (.not. cnvflg(i))
816  enddo
817  if(totflg) return
818 !!
819 c
820 c calculate convective inhibition
821 c
823  do k = 2, km1
824  do i = 1, im
825  if (cnvflg(i)) then
826  if(k > kb(i) .and. k < kbcon1(i)) then
827  dz1 = zo(i,k+1) - zo(i,k)
828  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
829  rfact = 1. + delta * cp * gamma
830  & * to(i,k) / hvap
831  cina(i) = cina(i) +
832 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
833  & dz1 * (g / (cp * to(i,k)))
834  & * dbyo(i,k) / (1. + gamma)
835  & * rfact
836  val = 0.
837  cina(i) = cina(i) +
838 ! & dz1 * eta(i,k) * g * delta *
839  & dz1 * g * delta *
840  & max(val,(qeso(i,k) - qo(i,k)))
841  endif
842  endif
843  enddo
844  enddo
846  do i = 1, im
847  if(cnvflg(i)) then
848 !
849  if(islimsk(i) == 1) then
850  w1 = w1l
851  w2 = w2l
852  w3 = w3l
853  w4 = w4l
854  else
855  w1 = w1s
856  w2 = w2s
857  w3 = w3s
858  w4 = w4s
859  endif
860  if(pdot(i) <= w4) then
861  tem = (pdot(i) - w4) / (w3 - w4)
862  elseif(pdot(i) >= -w4) then
863  tem = - (pdot(i) + w4) / (w4 - w3)
864  else
865  tem = 0.
866  endif
867 
868  val1 = -1.
869  tem = max(tem,val1)
870  val2 = 1.
871  tem = min(tem,val2)
872  tem = 1. - tem
873  tem1= .5*(cinacrmx-cinacrmn)
874  cinacr = cinacrmx - tem * tem1
875 !
876 ! cinacr = cinacrmx
877  if(cina(i) < cinacr) cnvflg(i) = .false.
878  endif
879  enddo
880 !!
881  totflg = .true.
882  do i=1,im
883  totflg = totflg .and. (.not. cnvflg(i))
884  enddo
885  if(totflg) return
886 !!
887 c
888 c determine first guess cloud top as the level of zero buoyancy
889 c
891  do i = 1, im
892  flg(i) = cnvflg(i)
893  ktcon(i) = 1
894  enddo
895  do k = 2, km1
896  do i = 1, im
897  if (flg(i) .and. k < kmax(i)) then
898  if(k > kbcon1(i) .and. dbyo(i,k) < 0.) then
899  ktcon(i) = k
900  flg(i) = .false.
901  endif
902  endif
903  enddo
904  enddo
905 c
906  do i = 1, im
907  if(cnvflg(i)) then
908  if(ktcon(i) == 1 .and. ktconn(i) > 1) then
909  ktcon(i) = ktconn(i)
910  endif
911  tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
912  if(tem < cthk) cnvflg(i) = .false.
913  endif
914  enddo
915 !!
916  totflg = .true.
917  do i = 1, im
918  totflg = totflg .and. (.not. cnvflg(i))
919  enddo
920  if(totflg) return
921 !!
922 c
923 c search for downdraft originating level above theta-e minimum
924 c
926  do i = 1, im
927  if(cnvflg(i)) then
928  hmin(i) = heo(i,kbcon1(i))
929  lmin(i) = kbmax(i)
930  jmin(i) = kbmax(i)
931  endif
932  enddo
933  do k = 2, km1
934  do i = 1, im
935  if (cnvflg(i) .and. k <= kbmax(i)) then
936  if(k > kbcon1(i) .and. heo(i,k) < hmin(i)) then
937  lmin(i) = k + 1
938  hmin(i) = heo(i,k)
939  endif
940  endif
941  enddo
942  enddo
943 c
944 c make sure that jmin is within the cloud
945 c
946  do i = 1, im
947  if(cnvflg(i)) then
948  jmin(i) = min(lmin(i),ktcon(i)-1)
949  jmin(i) = max(jmin(i),kbcon1(i)+1)
950  if(jmin(i) >= ktcon(i)) cnvflg(i) = .false.
951  endif
952  enddo
953 c
954 c specify upper limit of mass flux at cloud base
955 c
957  do i = 1, im
958  if(cnvflg(i)) then
959 ! xmbmax(i) = .1
960 !
961  k = kbcon(i)
962  dp = 1000. * del(i,k)
963  xmbmax(i) = dp / (g * dt2)
964 !
965 ! mbdt(i) = 0.1 * dp / g
966 !
967 ! tem = dp / (g * dt2)
968 ! xmbmax(i) = min(tem, xmbmax(i))
969  endif
970  enddo
971 c
972 c compute cloud moisture property and precipitation
973 c
975  do i = 1, im
976  if (cnvflg(i)) then
977 ! aa1(i) = 0.
978  qcko(i,kb(i)) = qo(i,kb(i))
979  qrcko(i,kb(i)) = qo(i,kb(i))
980 ! rhbar(i) = 0.
981  endif
982  enddo
984  do k = 2, km1
985  do i = 1, im
986  if (cnvflg(i)) then
987  if(k > kb(i) .and. k < ktcon(i)) then
988  dz = zi(i,k) - zi(i,k-1)
989  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
990  qrch = qeso(i,k)
991  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
992 cj
993  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
994  tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
995  factor = 1. + tem - tem1
996  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
997  & (qo(i,k)+qo(i,k-1)))/factor
998  qrcko(i,k) = qcko(i,k)
999 cj
1000  dq = eta(i,k) * (qcko(i,k) - qrch)
1001 c
1002 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
1003 c
1004 c check if there is excess moisture to release latent heat
1005 c
1006  if(k >= kbcon(i) .and. dq > 0.) then
1007  etah = .5 * (eta(i,k) + eta(i,k-1))
1008  dp = 1000. * del(i,k)
1009  if(ncloud > 0 .and. k > jmin(i)) then
1010  ptem = c0t(i,k) + c1
1011  qlk = dq / (eta(i,k) + etah * ptem * dz)
1012  dellal(i,k) = etah * c1 * dz * qlk * g / dp
1013  else
1014  qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1015  endif
1016 ! aa1(i) = aa1(i) - dz * g * qlk * etah
1017 ! aa1(i) = aa1(i) - dz * g * qlk
1018  buo(i,k) = buo(i,k) - g * qlk
1019  qcko(i,k) = qlk + qrch
1020  pwo(i,k) = etah * c0t(i,k) * dz * qlk
1021  pwavo(i) = pwavo(i) + pwo(i,k)
1022 ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
1023  cnvwt(i,k) = etah * qlk * g / dp
1024  endif
1025 !
1026 ! compute buoyancy and drag for updraft velocity
1027 !
1028  if(k >= kbcon(i)) then
1029  rfact = 1. + delta * cp * gamma
1030  & * to(i,k) / hvap
1031  buo(i,k) = buo(i,k) + (g / (cp * to(i,k)))
1032  & * dbyo(i,k) / (1. + gamma)
1033  & * rfact
1034  val = 0.
1035  buo(i,k) = buo(i,k) + g * delta *
1036  & max(val,(qeso(i,k) - qo(i,k)))
1037  drag(i,k) = max(xlamue(i,k),xlamud(i,k))
1038  endif
1039 !
1040  endif
1041  endif
1042  enddo
1043  enddo
1044 c
1045 ! do i = 1, im
1046 ! if(cnvflg(i)) then
1047 ! indx = ktcon(i) - kb(i) - 1
1048 ! rhbar(i) = rhbar(i) / float(indx)
1049 ! endif
1050 ! enddo
1051 c
1052 c calculate cloud work function
1053 c
1054 ! do k = 2, km1
1055 ! do i = 1, im
1056 ! if (cnvflg(i)) then
1057 ! if(k >= kbcon(i) .and. k < ktcon(i)) then
1058 ! dz1 = zo(i,k+1) - zo(i,k)
1059 ! gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1060 ! rfact = 1. + delta * cp * gamma
1061 ! & * to(i,k) / hvap
1062 ! aa1(i) = aa1(i) +
1063 !! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
1064 ! & dz1 * (g / (cp * to(i,k)))
1065 ! & * dbyo(i,k) / (1. + gamma)
1066 ! & * rfact
1067 ! val = 0.
1068 ! aa1(i) = aa1(i) +
1069 !! & dz1 * eta(i,k) * g * delta *
1070 ! & dz1 * g * delta *
1071 ! & max(val,(qeso(i,k) - qo(i,k)))
1072 ! endif
1073 ! endif
1074 ! enddo
1075 ! enddo
1076 !
1077 ! calculate cloud work function
1078 !
1084  do i = 1, im
1085  if (cnvflg(i)) then
1086  aa1(i) = 0.
1087  endif
1088  enddo
1089  do k = 2, km1
1090  do i = 1, im
1091  if (cnvflg(i)) then
1092  if(k >= kbcon(i) .and. k < ktcon(i)) then
1093  dz1 = zo(i,k+1) - zo(i,k)
1094 ! aa1(i) = aa1(i) + buo(i,k) * dz1 * eta(i,k)
1095  aa1(i) = aa1(i) + buo(i,k) * dz1
1096  endif
1097  endif
1098  enddo
1099  enddo
1100 !
1102  do i = 1, im
1103  if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
1104  enddo
1105 !!
1106  totflg = .true.
1107  do i=1,im
1108  totflg = totflg .and. (.not. cnvflg(i))
1109  enddo
1110  if(totflg) return
1111 !!
1112 c
1113 c estimate the onvective overshooting as the level
1114 c where the [aafac * cloud work function] becomes zero,
1115 c which is the final cloud top
1116 c
1118  do i = 1, im
1119  if (cnvflg(i)) then
1120  aa2(i) = aafac * aa1(i)
1121  endif
1122  enddo
1123 c
1124  do i = 1, im
1125  flg(i) = cnvflg(i)
1126  ktcon1(i) = kmax(i)
1127  enddo
1128  do k = 2, km1
1129  do i = 1, im
1130  if (flg(i)) then
1131  if(k >= ktcon(i) .and. k < kmax(i)) then
1132  dz1 = zo(i,k+1) - zo(i,k)
1133  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1134  rfact = 1. + delta * cp * gamma
1135  & * to(i,k) / hvap
1136  aa2(i) = aa2(i) +
1137 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
1138  & dz1 * (g / (cp * to(i,k)))
1139  & * dbyo(i,k) / (1. + gamma)
1140  & * rfact
1141 ! val = 0.
1142 ! aa2(i) = aa2(i) +
1143 !! & dz1 * eta(i,k) * g * delta *
1144 ! & dz1 * g * delta *
1145 ! & max(val,(qeso(i,k) - qo(i,k)))
1146  if(aa2(i) < 0.) then
1147  ktcon1(i) = k
1148  flg(i) = .false.
1149  endif
1150  endif
1151  endif
1152  enddo
1153  enddo
1154 c
1155 c compute cloud moisture property, detraining cloud water
1156 c and precipitation in overshooting layers
1157 c
1159  do k = 2, km1
1160  do i = 1, im
1161  if (cnvflg(i)) then
1162  if(k >= ktcon(i) .and. k < ktcon1(i)) then
1163  dz = zi(i,k) - zi(i,k-1)
1164  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1165  qrch = qeso(i,k)
1166  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1167 cj
1168  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1169  tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
1170  factor = 1. + tem - tem1
1171  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
1172  & (qo(i,k)+qo(i,k-1)))/factor
1173  qrcko(i,k) = qcko(i,k)
1174 cj
1175  dq = eta(i,k) * (qcko(i,k) - qrch)
1176 c
1177 c check if there is excess moisture to release latent heat
1178 c
1179  if(dq > 0.) then
1180  etah = .5 * (eta(i,k) + eta(i,k-1))
1181  dp = 1000. * del(i,k)
1182  if(ncloud > 0) then
1183  ptem = c0t(i,k) + c1
1184  qlk = dq / (eta(i,k) + etah * ptem * dz)
1185  dellal(i,k) = etah * c1 * dz * qlk * g / dp
1186  else
1187  qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1188  endif
1189  qcko(i,k) = qlk + qrch
1190  pwo(i,k) = etah * c0t(i,k) * dz * qlk
1191  pwavo(i) = pwavo(i) + pwo(i,k)
1192 ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
1193  cnvwt(i,k) = etah * qlk * g / dp
1194  endif
1195  endif
1196  endif
1197  enddo
1198  enddo
1199 !
1200 ! compute updraft velocity square(wu2)
1202 !
1203 ! bb1 = 2. * (1.+bet1*cd1)
1204 ! bb2 = 2. / (f1*(1.+gam1))
1205 !
1206 ! bb1 = 3.9
1207 ! bb2 = 0.67
1208 !
1209 ! bb1 = 2.0
1210 ! bb2 = 4.0
1211 !
1212  bb1 = 4.0
1213  bb2 = 0.8
1214 !
1215 ! do i = 1, im
1216 ! if (cnvflg(i)) then
1217 ! k = kbcon1(i)
1218 ! tem = po(i,k) / (rd * to(i,k))
1219 ! wucb = -0.01 * dot(i,k) / (tem * g)
1220 ! if(wucb > 0.) then
1221 ! wu2(i,k) = wucb * wucb
1222 ! else
1223 ! wu2(i,k) = 0.
1224 ! endif
1225 ! endif
1226 ! enddo
1227  do k = 2, km1
1228  do i = 1, im
1229  if (cnvflg(i)) then
1230  if(k > kbcon1(i) .and. k < ktcon(i)) then
1231  dz = zi(i,k) - zi(i,k-1)
1232  tem = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz
1233  tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz
1234  ptem = (1. - tem) * wu2(i,k-1)
1235  ptem1 = 1. + tem
1236  wu2(i,k) = (ptem + tem1) / ptem1
1237  wu2(i,k) = max(wu2(i,k), 0.)
1238  endif
1239  endif
1240  enddo
1241  enddo
1242 !
1243 ! compute updraft velocity average over the whole cumulus
1244 !
1246  do i = 1, im
1247  wc(i) = 0.
1248  sumx(i) = 0.
1249  enddo
1250  do k = 2, km1
1251  do i = 1, im
1252  if (cnvflg(i)) then
1253  if(k > kbcon1(i) .and. k < ktcon(i)) then
1254  dz = zi(i,k) - zi(i,k-1)
1255  tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1256  wc(i) = wc(i) + tem * dz
1257  sumx(i) = sumx(i) + dz
1258  endif
1259  endif
1260  enddo
1261  enddo
1262  do i = 1, im
1263  if(cnvflg(i)) then
1264  if(sumx(i) == 0.) then
1265  cnvflg(i)=.false.
1266  else
1267  wc(i) = wc(i) / sumx(i)
1268  endif
1269  val = 1.e-4
1270  if (wc(i) < val) cnvflg(i)=.false.
1271  endif
1272  enddo
1273 c
1274 c exchange ktcon with ktcon1
1275 c
1277  do i = 1, im
1278  if(cnvflg(i)) then
1279  kk = ktcon(i)
1280  ktcon(i) = ktcon1(i)
1281  ktcon1(i) = kk
1282  endif
1283  enddo
1284 c
1285 c this section is ready for cloud water
1286 c
1288  if(ncloud > 0) then
1289 c
1290 c compute liquid and vapor separation at cloud top
1291 c
1292  do i = 1, im
1293  if(cnvflg(i)) then
1294  k = ktcon(i) - 1
1295  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1296  qrch = qeso(i,k)
1297  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1298  dq = qcko(i,k) - qrch
1299 c
1300 c check if there is excess moisture to release latent heat
1301 c
1302  if(dq > 0.) then
1303  qlko_ktcon(i) = dq
1304  qcko(i,k) = qrch
1305  endif
1306  endif
1307  enddo
1308  endif
1309 c
1310 ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then
1311 ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
1312 ccccc endif
1313 c
1314 c------- downdraft calculations
1315 c
1316 c--- compute precipitation efficiency in terms of windshear
1317 c
1324  do i = 1, im
1325  if(cnvflg(i)) then
1326  vshear(i) = 0.
1327  endif
1328  enddo
1329  do k = 2, km
1330  do i = 1, im
1331  if (cnvflg(i)) then
1332  if(k > kb(i) .and. k <= ktcon(i)) then
1333  shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
1334  & + (vo(i,k)-vo(i,k-1)) ** 2)
1335  vshear(i) = vshear(i) + shear
1336  endif
1337  endif
1338  enddo
1339  enddo
1340  do i = 1, im
1341  if(cnvflg(i)) then
1342  vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
1343  e1=1.591-.639*vshear(i)
1344  & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1345  edt(i)=1.-e1
1346  val = .9
1347  edt(i) = min(edt(i),val)
1348  val = .0
1349  edt(i) = max(edt(i),val)
1350  edto(i)=edt(i)
1351  edtx(i)=edt(i)
1352  endif
1353  enddo
1354 c
1355 c determine detrainment rate between 1 and kbcon
1356 c
1362  do i = 1, im
1363  if(cnvflg(i)) then
1364  sumx(i) = 0.
1365  endif
1366  enddo
1367  do k = 1, km1
1368  do i = 1, im
1369  if(cnvflg(i)) then
1370  if(k >= 1 .and. k < kbcon(i)) then
1371  dz = zi(i,k+1) - zi(i,k)
1372  sumx(i) = sumx(i) + dz
1373  endif
1374  endif
1375  enddo
1376  enddo
1377  do i = 1, im
1378  beta = betas
1379  if(islimsk(i) == 1) beta = betal
1380  if(cnvflg(i)) then
1381  dz = (sumx(i)+zi(i,1))/float(kbcon(i))
1382  tem = 1./float(kbcon(i))
1383  xlamd(i) = (1.-beta**tem)/dz
1384  endif
1385  enddo
1386 c
1387 c determine downdraft mass flux
1388 c
1390  do k = km1, 1, -1
1391  do i = 1, im
1392  if (cnvflg(i) .and. k <= kmax(i)-1) then
1393  if(k < jmin(i) .and. k >= kbcon(i)) then
1394  dz = zi(i,k+1) - zi(i,k)
1395  ptem = xlamdd - xlamde
1396  etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1397  else if(k < kbcon(i)) then
1398  dz = zi(i,k+1) - zi(i,k)
1399  ptem = xlamd(i) + xlamdd - xlamde
1400  etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1401  endif
1402  endif
1403  enddo
1404  enddo
1405 c
1406 c--- downdraft moisture properties
1407 c
1409  do i = 1, im
1410  if(cnvflg(i)) then
1411  jmn = jmin(i)
1412  hcdo(i,jmn) = heo(i,jmn)
1413  qcdo(i,jmn) = qo(i,jmn)
1414  qrcdo(i,jmn)= qo(i,jmn)
1415  ucdo(i,jmn) = uo(i,jmn)
1416  vcdo(i,jmn) = vo(i,jmn)
1417  pwevo(i) = 0.
1418  endif
1419  enddo
1420 cj
1422  do k = km1, 1, -1
1423  do i = 1, im
1424  if (cnvflg(i) .and. k < jmin(i)) then
1425  dz = zi(i,k+1) - zi(i,k)
1426  if(k >= kbcon(i)) then
1427  tem = xlamde * dz
1428  tem1 = 0.5 * xlamdd * dz
1429  else
1430  tem = xlamde * dz
1431  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1432  endif
1433  factor = 1. + tem - tem1
1434  hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
1435  & (heo(i,k)+heo(i,k+1)))/factor
1436  dbyo(i,k) = hcdo(i,k) - heso(i,k)
1437 !
1438  tem = 0.5 * cm * tem
1439  factor = 1. + tem
1440  ptem = tem - pgcon
1441  ptem1= tem + pgcon
1442  ucdo(i,k) = ((1.-tem)*ucdo(i,k+1)+ptem*uo(i,k+1)
1443  & +ptem1*uo(i,k))/factor
1444  vcdo(i,k) = ((1.-tem)*vcdo(i,k+1)+ptem*vo(i,k+1)
1445  & +ptem1*vo(i,k))/factor
1446  endif
1447  enddo
1448  enddo
1449 c
1451  do k = km1, 1, -1
1452  do i = 1, im
1453  if (cnvflg(i) .and. k < jmin(i)) then
1454  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1455  qrcdo(i,k) = qeso(i,k)+
1456  & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
1457 ! detad = etad(i,k+1) - etad(i,k)
1458 cj
1459  dz = zi(i,k+1) - zi(i,k)
1460  if(k >= kbcon(i)) then
1461  tem = xlamde * dz
1462  tem1 = 0.5 * xlamdd * dz
1463  else
1464  tem = xlamde * dz
1465  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1466  endif
1467  factor = 1. + tem - tem1
1468  qcdo(i,k) = ((1.-tem1)*qrcdo(i,k+1)+tem*0.5*
1469  & (qo(i,k)+qo(i,k+1)))/factor
1470 cj
1471 ! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
1472 ! & etad(i,k) * qrcdo(i,k)
1473 ! pwdo(i,k) = pwdo(i,k) - detad *
1474 ! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
1475 cj
1476  pwdo(i,k) = etad(i,k) * (qcdo(i,k) - qrcdo(i,k))
1477  pwevo(i) = pwevo(i) + pwdo(i,k)
1478  endif
1479  enddo
1480  enddo
1481 c
1482 c--- final downdraft strength dependent on precip
1483 c--- efficiency(edt), normalized condensate(pwav), and
1484 c--- evaporate(pwev)
1485 c
1487  do i = 1, im
1488  edtmax = edtmaxl
1489  if(islimsk(i) == 0) edtmax = edtmaxs
1490  if(cnvflg(i)) then
1491  if(pwevo(i) < 0.) then
1492  edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1493  edto(i) = min(edto(i),edtmax)
1494  else
1495  edto(i) = 0.
1496  endif
1497  endif
1498  enddo
1499 c
1500 c--- downdraft cloudwork functions
1501 c
1503  do k = km1, 1, -1
1504  do i = 1, im
1505  if (cnvflg(i) .and. k < jmin(i)) then
1506  gamma = el2orc * qeso(i,k) / to(i,k)**2
1507  dhh=hcdo(i,k)
1508  dt=to(i,k)
1509  dg=gamma
1510  dh=heso(i,k)
1511  dz=-1.*(zo(i,k+1)-zo(i,k))
1512 ! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k)
1513  aa1(i)=aa1(i)+edto(i)*dz
1514  & *(g/(cp*dt))*((dhh-dh)/(1.+dg))
1515  & *(1.+delta*cp*dg*dt/hvap)
1516  val=0.
1517 ! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k)
1518  aa1(i)=aa1(i)+edto(i)*dz
1519  & *g*delta*max(val,(qeso(i,k)-qo(i,k)))
1520  endif
1521  enddo
1522  enddo
1524  do i = 1, im
1525  if(cnvflg(i) .and. aa1(i) <= 0.) then
1526  cnvflg(i) = .false.
1527  endif
1528  enddo
1529 !!
1530  totflg = .true.
1531  do i=1,im
1532  totflg = totflg .and. (.not. cnvflg(i))
1533  enddo
1534  if(totflg) return
1535 !!
1536 c
1537 c--- what would the change be, that a cloud with unit mass
1538 c--- will do to the environment?
1539 c
1541  do k = 1, km
1542  do i = 1, im
1543  if(cnvflg(i) .and. k <= kmax(i)) then
1544  dellah(i,k) = 0.
1545  dellaq(i,k) = 0.
1546  dellau(i,k) = 0.
1547  dellav(i,k) = 0.
1548  endif
1549  enddo
1550  enddo
1551  do i = 1, im
1552  if(cnvflg(i)) then
1553  dp = 1000. * del(i,1)
1554  dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1)
1555  & - heo(i,1)) * g / dp
1556  dellaq(i,1) = edto(i) * etad(i,1) * (qrcdo(i,1)
1557  & - qo(i,1)) * g / dp
1558  dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1)
1559  & - uo(i,1)) * g / dp
1560  dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1)
1561  & - vo(i,1)) * g / dp
1562  endif
1563  enddo
1564 c
1565 c--- changed due to subsidence and entrainment
1566 c
1567  do k = 2, km1
1568  do i = 1, im
1569  if (cnvflg(i) .and. k < ktcon(i)) then
1570  aup = 1.
1571  if(k <= kb(i)) aup = 0.
1572  adw = 1.
1573  if(k > jmin(i)) adw = 0.
1574  dp = 1000. * del(i,k)
1575  dz = zi(i,k) - zi(i,k-1)
1576 c
1577  dv1h = heo(i,k)
1578  dv2h = .5 * (heo(i,k) + heo(i,k-1))
1579  dv3h = heo(i,k-1)
1580  dv1q = qo(i,k)
1581  dv2q = .5 * (qo(i,k) + qo(i,k-1))
1582  dv3q = qo(i,k-1)
1583 c
1584  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1585  tem1 = 0.5 * (xlamud(i,k)+xlamud(i,k-1))
1586 c
1587  if(k <= kbcon(i)) then
1588  ptem = xlamde
1589  ptem1 = xlamd(i)+xlamdd
1590  else
1591  ptem = xlamde
1592  ptem1 = xlamdd
1593  endif
1594 cj
1595  dellah(i,k) = dellah(i,k) +
1596  & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h
1597  & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h
1598  & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz
1599  & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1600  & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz
1601  & ) *g/dp
1602 cj
1603  dellaq(i,k) = dellaq(i,k) +
1604  & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q
1605  & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q
1606  & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz
1607  & + aup*tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz
1608  & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qcdo(i,k-1))*dz
1609  & ) *g/dp
1610 cj
1611  tem1=eta(i,k)*(uo(i,k)-ucko(i,k))
1612  tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1))
1613  ptem1=etad(i,k)*(uo(i,k)-ucdo(i,k))
1614  ptem2=etad(i,k-1)*(uo(i,k-1)-ucdo(i,k-1))
1615  dellau(i,k) = dellau(i,k) +
1616  & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*g/dp
1617 cj
1618  tem1=eta(i,k)*(vo(i,k)-vcko(i,k))
1619  tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1))
1620  ptem1=etad(i,k)*(vo(i,k)-vcdo(i,k))
1621  ptem2=etad(i,k-1)*(vo(i,k-1)-vcdo(i,k-1))
1622  dellav(i,k) = dellav(i,k) +
1623  & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*g/dp
1624 cj
1625  endif
1626  enddo
1627  enddo
1628 c
1629 c------- cloud top
1630 c
1631  do i = 1, im
1632  if(cnvflg(i)) then
1633  indx = ktcon(i)
1634  dp = 1000. * del(i,indx)
1635  dv1h = heo(i,indx-1)
1636  dellah(i,indx) = eta(i,indx-1) *
1637  & (hcko(i,indx-1) - dv1h) * g / dp
1638  dv1q = qo(i,indx-1)
1639  dellaq(i,indx) = eta(i,indx-1) *
1640  & (qcko(i,indx-1) - dv1q) * g / dp
1641  dellau(i,indx) = eta(i,indx-1) *
1642  & (ucko(i,indx-1) - uo(i,indx-1)) * g / dp
1643  dellav(i,indx) = eta(i,indx-1) *
1644  & (vcko(i,indx-1) - vo(i,indx-1)) * g / dp
1645 c
1646 c cloud water
1647 c
1648  dellal(i,indx) = eta(i,indx-1) *
1649  & qlko_ktcon(i) * g / dp
1650  endif
1651  enddo
1652 c
1653 c------- final changed variable per unit mass flux
1654 c
1656 !
1657  do i = 1, im
1658  asqecflg(i) = cnvflg(i)
1659  if(asqecflg(i) .and. gdx(i) < dxcrtas) then
1660  asqecflg(i) = .false.
1661  endif
1662  enddo
1663 !
1665  do k = 1, km
1666  do i = 1, im
1667  if (asqecflg(i) .and. k <= kmax(i)) then
1668  if(k > ktcon(i)) then
1669  qo(i,k) = q1(i,k)
1670  to(i,k) = t1(i,k)
1671  endif
1672  if(k <= ktcon(i)) then
1673  qo(i,k) = dellaq(i,k) * mbdt(i) + q1(i,k)
1674  dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1675  to(i,k) = dellat * mbdt(i) + t1(i,k)
1676  val = 1.e-10
1677  qo(i,k) = max(qo(i,k), val )
1678  endif
1679  endif
1680  enddo
1681  enddo
1682 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1683 c
1684 c--- the above changed environment is now used to calulate the
1685 c--- effect the arbitrary cloud(with unit mass flux)
1686 c--- would have on the stability,
1687 c--- which then is used to calculate the real mass flux,
1688 c--- necessary to keep this change in balance with the large-scale
1689 c--- destabilization.
1690 c
1691 c--- environmental conditions again, first heights
1692 c
1696  do k = 1, km
1697  do i = 1, im
1698  if(asqecflg(i) .and. k <= kmax(i)) then
1699  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1700  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
1701  val = 1.e-8
1702  qeso(i,k) = max(qeso(i,k), val )
1703 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
1704  endif
1705  enddo
1706  enddo
1707 c
1708 c--- moist static energy
1709 c
1710 !! - Recalculate moist static energy and saturation moist static energy.
1711  do k = 1, km1
1712  do i = 1, im
1713  if(asqecflg(i) .and. k <= kmax(i)-1) then
1714  dz = .5 * (zo(i,k+1) - zo(i,k))
1715  dp = .5 * (pfld(i,k+1) - pfld(i,k))
1716  es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
1717  pprime = pfld(i,k+1) + epsm1 * es
1718  qs = eps * es / pprime
1719  dqsdp = - qs / pprime
1720  desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1721  dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
1722  gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1723  dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
1724  dq = dqsdt * dt + dqsdp * dp
1725  to(i,k) = to(i,k+1) + dt
1726  qo(i,k) = qo(i,k+1) + dq
1727  po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
1728  endif
1729  enddo
1730  enddo
1731  do k = 1, km1
1732  do i = 1, im
1733  if(asqecflg(i) .and. k <= kmax(i)-1) then
1734  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1735  qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
1736  val1 = 1.e-8
1737  qeso(i,k) = max(qeso(i,k), val1)
1738  val2 = 1.e-10
1739  qo(i,k) = max(qo(i,k), val2 )
1740 ! qo(i,k) = min(qo(i,k),qeso(i,k))
1741  heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
1742  & cp * to(i,k) + hvap * qo(i,k)
1743  heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
1744  & cp * to(i,k) + hvap * qeso(i,k)
1745  endif
1746  enddo
1747  enddo
1748  do i = 1, im
1749  if(asqecflg(i)) then
1750  k = kmax(i)
1751  heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
1752  heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
1753 c heo(i,k) = min(heo(i,k),heso(i,k))
1754  endif
1755  enddo
1756 c
1757 c**************************** static control
1758 c
1759 c------- moisture and cloud work functions
1760 c
1762  do i = 1, im
1763  if(asqecflg(i)) then
1764  xaa0(i) = 0.
1765  xpwav(i) = 0.
1766  endif
1767  enddo
1768 c
1769  do i = 1, im
1770  if(asqecflg(i)) then
1771  indx = kb(i)
1772  hcko(i,indx) = heo(i,indx)
1773  qcko(i,indx) = qo(i,indx)
1774  endif
1775  enddo
1776  do k = 2, km1
1777  do i = 1, im
1778  if (asqecflg(i)) then
1779  if(k > kb(i) .and. k <= ktcon(i)) then
1780  dz = zi(i,k) - zi(i,k-1)
1781  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1782  tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
1783  factor = 1. + tem - tem1
1784  hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
1785  & (heo(i,k)+heo(i,k-1)))/factor
1786  endif
1787  endif
1788  enddo
1789  enddo
1790  do k = 2, km1
1791  do i = 1, im
1792  if (asqecflg(i)) then
1793  if(k > kb(i) .and. k < ktcon(i)) then
1794  dz = zi(i,k) - zi(i,k-1)
1795  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1796  xdby = hcko(i,k) - heso(i,k)
1797  xqrch = qeso(i,k)
1798  & + gamma * xdby / (hvap * (1. + gamma))
1799 cj
1800  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1801  tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
1802  factor = 1. + tem - tem1
1803  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
1804  & (qo(i,k)+qo(i,k-1)))/factor
1805 cj
1806  dq = eta(i,k) * (qcko(i,k) - xqrch)
1807 c
1808  if(k >= kbcon(i) .and. dq > 0.) then
1809  etah = .5 * (eta(i,k) + eta(i,k-1))
1810  if(ncloud > 0 .and. k > jmin(i)) then
1811  ptem = c0t(i,k) + c1
1812  qlk = dq / (eta(i,k) + etah * ptem * dz)
1813  else
1814  qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1815  endif
1816  if(k < ktcon1(i)) then
1817 ! xaa0(i) = xaa0(i) - dz * g * qlk * etah
1818  xaa0(i) = xaa0(i) - dz * g * qlk
1819  endif
1820  qcko(i,k) = qlk + xqrch
1821  xpw = etah * c0t(i,k) * dz * qlk
1822  xpwav(i) = xpwav(i) + xpw
1823  endif
1824  endif
1825  if(k >= kbcon(i) .and. k < ktcon1(i)) then
1826  dz1 = zo(i,k+1) - zo(i,k)
1827  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1828  rfact = 1. + delta * cp * gamma
1829  & * to(i,k) / hvap
1830  xaa0(i) = xaa0(i)
1831 ! & + dz1 * eta(i,k) * (g / (cp * to(i,k)))
1832  & + dz1 * (g / (cp * to(i,k)))
1833  & * xdby / (1. + gamma)
1834  & * rfact
1835  val=0.
1836  xaa0(i) = xaa0(i) +
1837 ! & dz1 * eta(i,k) * g * delta *
1838  & dz1 * g * delta *
1839  & max(val,(qeso(i,k) - qo(i,k)))
1840  endif
1841  endif
1842  enddo
1843  enddo
1844 c
1845 c------- downdraft calculations
1846 c
1847 c--- downdraft moisture properties
1848 c
1850  do i = 1, im
1851  if(asqecflg(i)) then
1852  jmn = jmin(i)
1853  hcdo(i,jmn) = heo(i,jmn)
1854  qcdo(i,jmn) = qo(i,jmn)
1855  qrcd(i,jmn) = qo(i,jmn)
1856  xpwev(i) = 0.
1857  endif
1858  enddo
1859 cj
1860  do k = km1, 1, -1
1861  do i = 1, im
1862  if (asqecflg(i) .and. k < jmin(i)) then
1863  dz = zi(i,k+1) - zi(i,k)
1864  if(k >= kbcon(i)) then
1865  tem = xlamde * dz
1866  tem1 = 0.5 * xlamdd * dz
1867  else
1868  tem = xlamde * dz
1869  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1870  endif
1871  factor = 1. + tem - tem1
1872  hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
1873  & (heo(i,k)+heo(i,k+1)))/factor
1874  endif
1875  enddo
1876  enddo
1877 cj
1878  do k = km1, 1, -1
1879  do i = 1, im
1880  if (asqecflg(i) .and. k < jmin(i)) then
1881  dq = qeso(i,k)
1882  dt = to(i,k)
1883  gamma = el2orc * dq / dt**2
1884  dh = hcdo(i,k) - heso(i,k)
1885  qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
1886 ! detad = etad(i,k+1) - etad(i,k)
1887 cj
1888  dz = zi(i,k+1) - zi(i,k)
1889  if(k >= kbcon(i)) then
1890  tem = xlamde * dz
1891  tem1 = 0.5 * xlamdd * dz
1892  else
1893  tem = xlamde * dz
1894  tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1895  endif
1896  factor = 1. + tem - tem1
1897  qcdo(i,k) = ((1.-tem1)*qrcd(i,k+1)+tem*0.5*
1898  & (qo(i,k)+qo(i,k+1)))/factor
1899 cj
1900 ! xpwd = etad(i,k+1) * qcdo(i,k+1) -
1901 ! & etad(i,k) * qrcd(i,k)
1902 ! xpwd = xpwd - detad *
1903 ! & .5 * (qrcd(i,k) + qrcd(i,k+1))
1904 cj
1905  xpwd = etad(i,k) * (qcdo(i,k) - qrcd(i,k))
1906  xpwev(i) = xpwev(i) + xpwd
1907  endif
1908  enddo
1909  enddo
1910 c
1911  do i = 1, im
1912  edtmax = edtmaxl
1913  if(islimsk(i) == 0) edtmax = edtmaxs
1914  if(asqecflg(i)) then
1915  if(xpwev(i) >= 0.) then
1916  edtx(i) = 0.
1917  else
1918  edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1919  edtx(i) = min(edtx(i),edtmax)
1920  endif
1921  endif
1922  enddo
1923 c
1924 c
1925 c--- downdraft cloudwork functions
1926 c
1927 c
1928  do k = km1, 1, -1
1929  do i = 1, im
1930  if (asqecflg(i) .and. k < jmin(i)) then
1931  gamma = el2orc * qeso(i,k) / to(i,k)**2
1932  dhh=hcdo(i,k)
1933  dt= to(i,k)
1934  dg= gamma
1935  dh= heso(i,k)
1936  dz=-1.*(zo(i,k+1)-zo(i,k))
1937 ! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k)
1938  xaa0(i)=xaa0(i)+edtx(i)*dz
1939  & *(g/(cp*dt))*((dhh-dh)/(1.+dg))
1940  & *(1.+delta*cp*dg*dt/hvap)
1941  val=0.
1942 ! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k)
1943  xaa0(i)=xaa0(i)+edtx(i)*dz
1944  & *g*delta*max(val,(qeso(i,k)-qo(i,k)))
1945  endif
1946  enddo
1947  enddo
1948 c
1949 c calculate critical cloud work function
1950 c
1951 ! do i = 1, im
1952 ! if(cnvflg(i)) then
1953 ! if(pfld(i,ktcon(i)) < pcrit(15))then
1954 ! acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i)))
1955 ! & /(975.-pcrit(15))
1956 ! else if(pfld(i,ktcon(i)) > pcrit(1))then
1957 ! acrt(i)=acrit(1)
1958 ! else
1959 ! k = int((850. - pfld(i,ktcon(i)))/50.) + 2
1960 ! k = min(k,15)
1961 ! k = max(k,2)
1962 ! acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*
1963 ! & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
1964 ! endif
1965 ! endif
1966 ! enddo
1967 ! do i = 1, im
1968 ! if(cnvflg(i)) then
1969 ! if(islimsk(i) == 1) then
1970 ! w1 = w1l
1971 ! w2 = w2l
1972 ! w3 = w3l
1973 ! w4 = w4l
1974 ! else
1975 ! w1 = w1s
1976 ! w2 = w2s
1977 ! w3 = w3s
1978 ! w4 = w4s
1979 ! endif
1980 c
1981 c modify critical cloud workfunction by cloud base vertical velocity
1982 c
1983 ! if(pdot(i) <= w4) then
1984 ! acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
1985 ! elseif(pdot(i) >= -w4) then
1986 ! acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
1987 ! else
1988 ! acrtfct(i) = 0.
1989 ! endif
1990 ! val1 = -1.
1991 ! acrtfct(i) = max(acrtfct(i),val1)
1992 ! val2 = 1.
1993 ! acrtfct(i) = min(acrtfct(i),val2)
1994 ! acrtfct(i) = 1. - acrtfct(i)
1995 c
1996 c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
1997 c
1998 c if(rhbar(i) >= .8) then
1999 c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
2000 c endif
2001 c
2002 c modify adjustment time scale by cloud base vertical velocity
2003 c
2004 ! dtconv(i) = dt2 + max((1800. - dt2),0.) *
2005 ! & (pdot(i) - w2) / (w1 - w2)
2006 c dtconv(i) = max(dtconv(i), dt2)
2007 c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
2008 !
2009 ! dtconv(i) = max(dtconv(i),dtmin)
2010 ! dtconv(i) = min(dtconv(i),dtmax)
2011 c
2012 ! endif
2013 ! enddo
2014 !
2015 ! compute convective turn-over time
2016 !
2018  do i= 1, im
2019  if(cnvflg(i)) then
2020  tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
2021  dtconv(i) = tem / wc(i)
2022  tfac = 1. + gdx(i) / 75000.
2023  dtconv(i) = tfac * dtconv(i)
2024  dtconv(i) = max(dtconv(i),dtmin)
2025  dtconv(i) = min(dtconv(i),dtmax)
2026  endif
2027  enddo
2028 !
2030  do i= 1, im
2031  if(cnvflg(i)) then
2032  sumx(i) = 0.
2033  umean(i) = 0.
2034  endif
2035  enddo
2036  do k = 2, km1
2037  do i = 1, im
2038  if(cnvflg(i)) then
2039  if(k >= kbcon1(i) .and. k < ktcon1(i)) then
2040  dz = zi(i,k) - zi(i,k-1)
2041  tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
2042  umean(i) = umean(i) + tem * dz
2043  sumx(i) = sumx(i) + dz
2044  endif
2045  endif
2046  enddo
2047  enddo
2048  do i= 1, im
2049  if(cnvflg(i)) then
2050  umean(i) = umean(i) / sumx(i)
2051  umean(i) = max(umean(i), 1.)
2052  tauadv(i) = gdx(i) / umean(i)
2053  endif
2054  enddo
2057  do i= 1, im
2058  if(cnvflg(i) .and. .not.asqecflg(i)) then
2059  k = kbcon(i)
2060  rho = po(i,k)*100. / (rd*to(i,k))
2061  tfac = tauadv(i) / dtconv(i)
2062  tfac = min(tfac, 1.)
2063  xmb(i) = tfac*betaw*rho*wc(i)
2064  endif
2065  enddo
2071  do i= 1, im
2072  if(asqecflg(i)) then
2073 ! fld(i)=(aa1(i)-acrt(i)*acrtfct(i))/dtconv(i)
2074  fld(i)=aa1(i)/dtconv(i)
2075  if(fld(i) <= 0.) then
2076  asqecflg(i) = .false.
2077  cnvflg(i) = .false.
2078  endif
2079  endif
2085  if(asqecflg(i)) then
2086 c xaa0(i) = max(xaa0(i),0.)
2087  xk(i) = (xaa0(i) - aa1(i)) / mbdt(i)
2088  if(xk(i) >= 0.) then
2089  asqecflg(i) = .false.
2090  cnvflg(i) = .false.
2091  endif
2092  endif
2093 c
2094 c--- kernel, cloud base mass flux
2095 c
2102  if(asqecflg(i)) then
2103  tfac = tauadv(i) / dtconv(i)
2104  tfac = min(tfac, 1.)
2105  xmb(i) = -tfac * fld(i) / xk(i)
2106 ! xmb(i) = min(xmb(i),xmbmax(i))
2107  endif
2108  enddo
2109 !!
2111  totflg = .true.
2112  do i=1,im
2113  totflg = totflg .and. (.not. cnvflg(i))
2114  enddo
2115  if(totflg) return
2116 !!
2117 !
2119  do i = 1, im
2120  if(cnvflg(i)) then
2121  tem = min(max(xlamx(i), 7.e-5), 3.e-4)
2122  tem = 0.2 / tem
2123  tem1 = 3.14 * tem * tem
2124  sigmagfm(i) = tem1 / garea(i)
2125  sigmagfm(i) = max(sigmagfm(i), 0.001)
2126  sigmagfm(i) = min(sigmagfm(i), 0.999)
2127  endif
2128  enddo
2129 !
2131  do i = 1, im
2132  if(cnvflg(i)) then
2133  if (gdx(i) < dxcrtuf) then
2134  scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
2135  scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
2136  else
2137  scaldfunc(i) = 1.0
2138  endif
2139  xmb(i) = xmb(i) * scaldfunc(i)
2140  xmb(i) = min(xmb(i),xmbmax(i))
2141  endif
2142  enddo
2143 c
2144 c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
2145 c
2146  do k = 1, km
2147  do i = 1, im
2148  if (cnvflg(i) .and. k <= kmax(i)) then
2149  to(i,k) = t1(i,k)
2150  qo(i,k) = q1(i,k)
2151  uo(i,k) = u1(i,k)
2152  vo(i,k) = v1(i,k)
2153  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
2154  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
2155  val = 1.e-8
2156  qeso(i,k) = max(qeso(i,k), val )
2157  endif
2158  enddo
2159  enddo
2160 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2161 c
2162 c--- feedback: simply the changes from the cloud with unit mass flux
2163 c--- multiplied by the mass flux necessary to keep the
2164 c--- equilibrium with the larger-scale.
2165 c
2170  do i = 1, im
2171  delhbar(i) = 0.
2172  delqbar(i) = 0.
2173  deltbar(i) = 0.
2174  delubar(i) = 0.
2175  delvbar(i) = 0.
2176  qcond(i) = 0.
2177  enddo
2178  do k = 1, km
2179  do i = 1, im
2180  if (cnvflg(i) .and. k <= kmax(i)) then
2181  if(k <= ktcon(i)) then
2182  dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
2183  t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
2184  q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
2185 ! tem = 1./rcs(i)
2186 ! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
2187 ! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
2188  u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
2189  v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
2190  dp = 1000. * del(i,k)
2191  delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
2192  delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
2193  deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
2194  delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
2195  delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
2196  endif
2197  endif
2198  enddo
2199  enddo
2201  do k = 1, km
2202  do i = 1, im
2203  if (cnvflg(i) .and. k <= kmax(i)) then
2204  if(k <= ktcon(i)) then
2205  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
2206  qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
2207  val = 1.e-8
2208  qeso(i,k) = max(qeso(i,k), val )
2209  endif
2210  endif
2211  enddo
2212  enddo
2213 c
2215  do i = 1, im
2216  rntot(i) = 0.
2217  delqev(i) = 0.
2218  delq2(i) = 0.
2219  flg(i) = cnvflg(i)
2220  enddo
2221  do k = km, 1, -1
2222  do i = 1, im
2223  if (cnvflg(i) .and. k <= kmax(i)) then
2224  if(k < ktcon(i)) then
2225  aup = 1.
2226  if(k <= kb(i)) aup = 0.
2227  adw = 1.
2228  if(k >= jmin(i)) adw = 0.
2229  rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
2230  rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
2231  endif
2232  endif
2233  enddo
2234  enddo
2238  do k = km, 1, -1
2239  do i = 1, im
2240  if (k <= kmax(i)) then
2241  deltv(i) = 0.
2242  delq(i) = 0.
2243  qevap(i) = 0.
2244  if(cnvflg(i) .and. k < ktcon(i)) then
2245  aup = 1.
2246  if(k <= kb(i)) aup = 0.
2247  adw = 1.
2248  if(k >= jmin(i)) adw = 0.
2249  rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
2250  rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
2251  endif
2252  if(flg(i) .and. k < ktcon(i)) then
2253  evef = edt(i) * evfact
2254  if(islimsk(i) == 1) evef=edt(i) * evfactl
2255 ! if(islimsk(i) == 1) evef=.07
2256 c if(islimsk(i) == 1) evef = 0.
2257  qcond(i) = evef * (q1(i,k) - qeso(i,k))
2258  & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
2259  dp = 1000. * del(i,k)
2260  if(rn(i) > 0. .and. qcond(i) < 0.) then
2261  qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
2262  qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
2263  delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
2264  endif
2265  if(rn(i) > 0. .and. qcond(i) < 0. .and.
2266  & delq2(i) > rntot(i)) then
2267  qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
2268  flg(i) = .false.
2269  endif
2270  if(rn(i) > 0. .and. qevap(i) > 0.) then
2271  q1(i,k) = q1(i,k) + qevap(i)
2272  t1(i,k) = t1(i,k) - elocp * qevap(i)
2273  rn(i) = rn(i) - .001 * qevap(i) * dp / g
2274  deltv(i) = - elocp*qevap(i)/dt2
2275  delq(i) = + qevap(i)/dt2
2276  delqev(i) = delqev(i) + .001*dp*qevap(i)/g
2277  endif
2278  delqbar(i) = delqbar(i) + delq(i)*dp/g
2279  deltbar(i) = deltbar(i) + deltv(i)*dp/g
2280  endif
2281  endif
2282  enddo
2283  enddo
2284 cj
2285 ! do i = 1, im
2286 ! if(me == 31 .and. cnvflg(i)) then
2287 ! if(cnvflg(i)) then
2288 ! print *, ' deep delhbar, delqbar, deltbar = ',
2289 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
2290 ! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
2291 ! print *, ' precip =', hvap*rn(i)*1000./dt2
2292 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
2293 ! endif
2294 ! enddo
2295 c
2296 c precipitation rate converted to actual precip
2297 c in unit of m instead of kg
2298 c
2299  do i = 1, im
2300  if(cnvflg(i)) then
2301 c
2302 c in the event of upper level rain evaporation and lower level downdraft
2303 c moistening, rn can become negative, in this case, we back out of the
2304 c heating and the moistening
2305 c
2306  if(rn(i) < 0. .and. .not.flg(i)) rn(i) = 0.
2307  if(rn(i) <= 0.) then
2308  rn(i) = 0.
2309  else
2310  ktop(i) = ktcon(i)
2311  kbot(i) = kbcon(i)
2312  kcnv(i) = 1
2313  cldwrk(i) = aa1(i)
2314  endif
2315  endif
2316  enddo
2317 c
2318 c convective cloud water
2319 c
2321  do k = 1, km
2322  do i = 1, im
2323  if (cnvflg(i) .and. rn(i) > 0.) then
2324  if (k >= kbcon(i) .and. k < ktcon(i)) then
2325  cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
2326  endif
2327  endif
2328  enddo
2329  enddo
2330 c
2331 c convective cloud cover
2332 c
2334  do k = 1, km
2335  do i = 1, im
2336  if (cnvflg(i) .and. rn(i) > 0.) then
2337  if (k >= kbcon(i) .and. k < ktcon(i)) then
2338  cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
2339  cnvc(i,k) = min(cnvc(i,k), 0.6)
2340  cnvc(i,k) = max(cnvc(i,k), 0.0)
2341  endif
2342  endif
2343  enddo
2344  enddo
2345 c
2346 c cloud water
2347 c
2349  if (ncloud > 0) then
2350 !
2351  do k = 1, km
2352  do i = 1, im
2353  if (cnvflg(i) .and. rn(i) > 0.) then
2354 ! if (k > kb(i) .and. k <= ktcon(i)) then
2355  if (k >= kbcon(i) .and. k <= ktcon(i)) then
2356  tem = dellal(i,k) * xmb(i) * dt2
2357  tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2358  if (ql(i,k,2) > -999.0) then
2359  ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
2360  ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
2361  else
2362  ql(i,k,1) = ql(i,k,1) + tem
2363  endif
2364  endif
2365  endif
2366  enddo
2367  enddo
2368 !
2369  endif
2370 c
2372  do k = 1, km
2373  do i = 1, im
2374  if(cnvflg(i) .and. rn(i) <= 0.) then
2375  if (k <= kmax(i)) then
2376  t1(i,k) = to(i,k)
2377  q1(i,k) = qo(i,k)
2378  u1(i,k) = uo(i,k)
2379  v1(i,k) = vo(i,k)
2380  endif
2381  endif
2382  enddo
2383  enddo
2384 !
2385 ! hchuang code change
2386 !
2388 !
2390  do k = 1, km
2391  do i = 1, im
2392  if(cnvflg(i) .and. rn(i) > 0.) then
2393  if(k >= kb(i) .and. k < ktop(i)) then
2394  ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
2395  endif
2396  endif
2397  enddo
2398  enddo
2400  do i = 1, im
2401  if(cnvflg(i) .and. rn(i) > 0.) then
2402  k = ktop(i)-1
2403  dt_mf(i,k) = ud_mf(i,k)
2404  endif
2405  enddo
2407  do k = 1, km
2408  do i = 1, im
2409  if(cnvflg(i) .and. rn(i) > 0.) then
2410  if(k >= 1 .and. k <= jmin(i)) then
2411  dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
2412  endif
2413  endif
2414  enddo
2415  enddo
2416 !!
2417  return
2418  end
samfdeepcnv
subroutine samfdeepcnv(im, ix, km, delt, delp, prslp, psp, phil, ql,
This subroutine contains the entirety of the SAMF deep convection scheme.
Definition: samfdeepcnv.f:76