GFS Operational Physics Documentation  Revision: 81451
shalcnv.f
Go to the documentation of this file.
1 
11 
14 
55  subroutine shalcnv(im,ix,km,jcap,delt,delp,prslp,psp,phil,ql,
56  & q1,t1,u1,v1,rn,kbot,ktop,kcnv,islimsk,
57  & dot,ncloud,hpbl,heat,evap,ud_mf,dt_mf,cnvw,cnvc)
58 ! & q1,t1,u1,v1,rcs,rn,kbot,ktop,kcnv,islimsk,
59 ! & dot,ncloud,hpbl,heat,evap,ud_mf,dt_mf,me)
60 !
61  use machine , only : kind_phys
62  use funcphys , only : fpvs
63  use physcons, grav => con_g, cp => con_cp, hvap => con_hvap
64  &, rv => con_rv, fv => con_fvirt, t0c => con_t0c
65  &, rd => con_rd, cvap => con_cvap, cliq => con_cliq
66  &, eps => con_eps, epsm1 => con_epsm1
67  implicit none
68 !
69  integer im, ix, km, jcap, ncloud,
70  & kbot(im), ktop(im), kcnv(im)
71 ! &, me
72  real(kind=kind_phys) delt
73  real(kind=kind_phys) psp(im), delp(ix,km), prslp(ix,km)
74  real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km),
75  & ql(ix,km,2),q1(ix,km), t1(ix,km),
76  & u1(ix,km), v1(ix,km), ! & u1(ix,km), v1(ix,km), rcs(im),
77  & rn(im),
78  & dot(ix,km), phil(ix,km), hpbl(im),
79  & heat(im), evap(im), cnvw(ix,km),
80  & cnvc(ix,km), ! hchuang code change mass flux output
81  & ud_mf(im,km),dt_mf(im,km)
82 !
83  integer i,j,indx, k, kk, km1
84  integer kpbl(im)
85  integer, dimension(im), intent(in) :: islimsk
86 !
87  real(kind=kind_phys) c0, dellat, delta,
88  & desdt,
89  & dp,
90  & dq, dqsdp, dqsdt, dt,
91  & dt2, dtmax, dtmin, dv1h,
92  & dv1q, dv2h, dv2q, dv1u,
93  & dv1v, dv2u, dv2v, dv3q,
94  & dv3h, dv3u, dv3v, clam,
95  & dz, dz1, e1,
96  & el2orc, elocp, aafac,
97  & es, etah, h1, dthk,
98  & evef, evfact, evfactl, fact1,
99  & fact2, factor, fjcap,
100  & g, gamma, pprime, betaw,
101  & qlk, qrch, qs, c1,
102  & rfact, shear, tem1,
103  & val, val1,
104  & val2, w1, w1l, w1s,
105  & w2, w2l, w2s, w3,
106  & w3l, w3s, w4, w4l,
107  & w4s, tem, ptem, ptem1,
108  & pgcon
109 !
110  integer kb(im), kbcon(im), kbcon1(im),
111  & ktcon(im), ktcon1(im),
112  & kbm(im), kmax(im)
113 !
114  real(kind=kind_phys) aa1(im),
115  & delhbar(im), delq(im), delq2(im),
116  & delqbar(im), delqev(im), deltbar(im),
117  & deltv(im), edt(im),
118  & wstar(im), sflx(im),
119  & pdot(im), po(im,km),
120  & qcond(im), qevap(im), hmax(im),
121  & rntot(im), vshear(im),
122  & xlamud(im), xmb(im), xmbmax(im),
123  & delubar(im), delvbar(im)
124 c
125  real(kind=kind_phys) cincr, cincrmax, cincrmin
126 cc
127 c physical parameters
128  parameter(g=grav)
129  parameter(elocp=hvap/cp,
130  & el2orc=hvap*hvap/(rv*cp))
131  parameter(c0=.002,c1=5.e-4,delta=fv)
132  parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
133  parameter(cincrmax=180.,cincrmin=120.,dthk=25.)
134  parameter(h1=0.33333333)
135 c local variables and arrays
136  real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
137  & uo(im,km), vo(im,km), qeso(im,km)
138 c cloud water
139 ! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
140  real(kind=kind_phys) qlko_ktcon(im), dellal(im,km),
141  & dbyo(im,km), zo(im,km), xlamue(im,km),
142  & heo(im,km), heso(im,km),
143  & dellah(im,km), dellaq(im,km),
144  & dellau(im,km), dellav(im,km), hcko(im,km),
145  & ucko(im,km), vcko(im,km), qcko(im,km),
146  & qrcko(im,km), eta(im,km),
147  & zi(im,km), pwo(im,km),
148  & tx1(im), cnvwt(im,km)
149 !
150  logical totflg, cnvflg(im), flg(im)
151 !
152  real(kind=kind_phys) tf, tcr, tcrf
153  parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
154 !
155 c-----------------------------------------------------------------------
156 !
157 !************************************************************************
158 ! convert input pa terms to cb terms -- moorthi
161  ps = psp * 0.001
162  prsl = prslp * 0.001
163  del = delp * 0.001
164 !************************************************************************
165 !
166  km1 = km - 1
167 c
168 c compute surface buoyancy flux
169 c
175  do i=1,im
176  sflx(i) = heat(i)+fv*t1(i,1)*evap(i)
177  enddo
178 c
179 c initialize arrays
180 c
182  do i=1,im
183  cnvflg(i) = .true.
184  if(kcnv(i).eq.1) cnvflg(i) = .false.
185  if(sflx(i).le.0.) cnvflg(i) = .false.
186  if(cnvflg(i)) then
187  kbot(i)=km+1
188  ktop(i)=0
189  endif
190  rn(i)=0.
191  kbcon(i)=km
192  ktcon(i)=1
193  kb(i)=km
194  pdot(i) = 0.
195  qlko_ktcon(i) = 0.
196  edt(i) = 0.
197  aa1(i) = 0.
198  vshear(i) = 0.
199  enddo
201 ! hchuang code change
202  do k = 1, km
203  do i = 1, im
204  ud_mf(i,k) = 0.
205  dt_mf(i,k) = 0.
206  enddo
207  enddo
208 !!
210  totflg = .true.
211  do i=1,im
212  totflg = totflg .and. (.not. cnvflg(i))
213  enddo
214  if(totflg) return
215 !!
216 c
218  dt2 = delt
219  val = 1200.
220  dtmin = max(dt2, val )
221  val = 3600.
222  dtmax = max(dt2, val )
223 c model tunable parameters are all here
224  clam = .3
225  aafac = .1
226  betaw = .03
227 c evef = 0.07
228  evfact = 0.3
229  evfactl = 0.3
230 !
231 ! pgcon = 0.7 ! gregory et al. (1997, qjrms)
232  pgcon = 0.55 ! zhang & wu (2003,jas)
233  fjcap = (float(jcap) / 126.) ** 2
234  val = 1.
235  fjcap = max(fjcap,val)
236  w1l = -8.e-3
237  w2l = -4.e-2
238  w3l = -5.e-3
239  w4l = -5.e-4
240  w1s = -2.e-4
241  w2s = -2.e-3
242  w3s = -1.e-3
243  w4s = -2.e-5
244 c
245 c define top layer for search of the downdraft originating layer
246 c and the maximum thetae for updraft
247 c
249  do i=1,im
250  kbm(i) = km
251  kmax(i) = km
252  tx1(i) = 1.0 / ps(i)
253  enddo
254 !
255  do k = 1, km
256  do i=1,im
257  if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
258  if (prsl(i,k)*tx1(i) .gt. 0.60) kmax(i) = k + 1
259  enddo
260  enddo
261  do i=1,im
262  kbm(i) = min(kbm(i),kmax(i))
263  enddo
264 c
265 c hydrostatic height assume zero terr and compute
266 c updraft entrainment rate as an inverse function of height
267 c
269  do k = 1, km
270  do i=1,im
271  zo(i,k) = phil(i,k) / g
272  enddo
273  enddo
275  do k = 1, km1
276  do i=1,im
277  zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
278  xlamue(i,k) = clam / zi(i,k)
279  enddo
280  enddo
281  do i=1,im
282  xlamue(i,km) = xlamue(i,km1)
283  enddo
284 c
285 c pbl height
286 c
288  do i=1,im
289  flg(i) = cnvflg(i)
290  kpbl(i)= 1
291  enddo
292  do k = 2, km1
293  do i=1,im
294  if (flg(i).and.zo(i,k).le.hpbl(i)) then
295  kpbl(i) = k
296  else
297  flg(i) = .false.
298  endif
299  enddo
300  enddo
301  do i=1,im
302  kpbl(i)= min(kpbl(i),kbm(i))
303  enddo
304 c
305 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
306 c convert surface pressure to mb from cb
307 c
309  do k = 1, km
310  do i = 1, im
311  if (cnvflg(i) .and. k .le. kmax(i)) then
312  pfld(i,k) = prsl(i,k) * 10.0
313  eta(i,k) = 1.
314  hcko(i,k) = 0.
315  qcko(i,k) = 0.
316  qrcko(i,k)= 0.
317  ucko(i,k) = 0.
318  vcko(i,k) = 0.
319  dbyo(i,k) = 0.
320  pwo(i,k) = 0.
321  dellal(i,k) = 0.
322  to(i,k) = t1(i,k)
323  qo(i,k) = q1(i,k)
324  uo(i,k) = u1(i,k)
325  vo(i,k) = v1(i,k)
326 ! uo(i,k) = u1(i,k) * rcs(i)
327 ! vo(i,k) = v1(i,k) * rcs(i)
328  cnvwt(i,k) = 0.
329  endif
330  enddo
331  enddo
332 c
333 c column variables
334 c p is pressure of the layer (mb)
335 c t is temperature at t-dt (k)..tn
336 c q is mixing ratio at t-dt (kg/kg)..qn
337 c to is temperature at t+dt (k)... this is after advection and turbulan
338 c qo is mixing ratio at t+dt (kg/kg)..q1
339 c
341  do k = 1, km
342  do i=1,im
343  if (cnvflg(i) .and. k .le. kmax(i)) then
344  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
345  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
346  val1 = 1.e-8
347  qeso(i,k) = max(qeso(i,k), val1)
348  val2 = 1.e-10
349  qo(i,k) = max(qo(i,k), val2 )
350 ! qo(i,k) = min(qo(i,k),qeso(i,k))
351 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
352  endif
353  enddo
354  enddo
355 c
356 c compute moist static energy
357 c
359  do k = 1, km
360  do i=1,im
361  if (cnvflg(i) .and. k .le. kmax(i)) then
362 ! tem = g * zo(i,k) + cp * to(i,k)
363  tem = phil(i,k) + cp * to(i,k)
364  heo(i,k) = tem + hvap * qo(i,k)
365  heso(i,k) = tem + hvap * qeso(i,k)
366 c heo(i,k) = min(heo(i,k),heso(i,k))
367  endif
368  enddo
369  enddo
370 c
371 c determine level with largest moist static energy within pbl
372 c this is the level where updraft starts
373 c
376  do i=1,im
377  if (cnvflg(i)) then
378  hmax(i) = heo(i,1)
379  kb(i) = 1
380  endif
381  enddo
382  do k = 2, km
383  do i=1,im
384  if (cnvflg(i).and.k.le.kpbl(i)) then
385  if(heo(i,k).gt.hmax(i)) then
386  kb(i) = k
387  hmax(i) = heo(i,k)
388  endif
389  endif
390  enddo
391  enddo
392 c
394  do k = 1, km1
395  do i=1,im
396  if (cnvflg(i) .and. k .le. kmax(i)-1) then
397  dz = .5 * (zo(i,k+1) - zo(i,k))
398  dp = .5 * (pfld(i,k+1) - pfld(i,k))
399  es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
400  pprime = pfld(i,k+1) + epsm1 * es
401  qs = eps * es / pprime
402  dqsdp = - qs / pprime
403  desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
404  dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
405  gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
406  dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
407  dq = dqsdt * dt + dqsdp * dp
408  to(i,k) = to(i,k+1) + dt
409  qo(i,k) = qo(i,k+1) + dq
410  po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
411  endif
412  enddo
413  enddo
414 !
416  do k = 1, km1
417  do i=1,im
418  if (cnvflg(i) .and. k .le. kmax(i)-1) then
419  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
420  qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
421  val1 = 1.e-8
422  qeso(i,k) = max(qeso(i,k), val1)
423  val2 = 1.e-10
424  qo(i,k) = max(qo(i,k), val2 )
425 ! qo(i,k) = min(qo(i,k),qeso(i,k))
426  heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
427  & cp * to(i,k) + hvap * qo(i,k)
428  heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
429  & cp * to(i,k) + hvap * qeso(i,k)
430  uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
431  vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
432  endif
433  enddo
434  enddo
435 c
436 c look for the level of free convection as cloud base
437 c!> - Search below the index "kbm" for the level of free convection (LFC) where the condition \f$h_b > h^*\f$ is first met, where \f$h_b, h^*\f$ are the state moist static energy at the parcel's starting level and saturation moist static energy, respectively. Set "kbcon" to the index of the LFC.
438  do i=1,im
439  flg(i) = cnvflg(i)
440  if(flg(i)) kbcon(i) = kmax(i)
441  enddo
442  do k = 2, km1
443  do i=1,im
444  if (flg(i).and.k.lt.kbm(i)) then
445  if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
446  kbcon(i) = k
447  flg(i) = .false.
448  endif
449  endif
450  enddo
451  enddo
452 c
454  do i=1,im
455  if(cnvflg(i)) then
456  if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
457  endif
458  enddo
459 !!
460  totflg = .true.
461  do i=1,im
462  totflg = totflg .and. (.not. cnvflg(i))
463  enddo
464  if(totflg) return
465 !!
466 c
467 c determine critical convective inhibition
468 c as a function of vertical velocity at cloud base.
469 c
471  do i=1,im
472  if(cnvflg(i)) then
473 ! pdot(i) = 10.* dot(i,kbcon(i))
474  pdot(i) = 0.01 * dot(i,kbcon(i)) ! now dot is in pa/s
475  endif
476  enddo
477  do i=1,im
478  if(cnvflg(i)) then
479  if(islimsk(i) == 1) then
480  w1 = w1l
481  w2 = w2l
482  w3 = w3l
483  w4 = w4l
484  else
485  w1 = w1s
486  w2 = w2s
487  w3 = w3s
488  w4 = w4s
489  endif
490  if(pdot(i).le.w4) then
491  ptem = (pdot(i) - w4) / (w3 - w4)
492  elseif(pdot(i).ge.-w4) then
493  ptem = - (pdot(i) + w4) / (w4 - w3)
494  else
495  ptem = 0.
496  endif
497  val1 = -1.
498  ptem = max(ptem,val1)
499  val2 = 1.
500  ptem = min(ptem,val2)
501  ptem = 1. - ptem
502  ptem1= .5*(cincrmax-cincrmin)
503  cincr = cincrmax - ptem * ptem1
504  tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
505  if(tem1.gt.cincr) then
506  cnvflg(i) = .false.
507  endif
508  endif
509  enddo
510 !!
511  totflg = .true.
512  do i=1,im
513  totflg = totflg .and. (.not. cnvflg(i))
514  enddo
515  if(totflg) return
516 !!
517 c
518 c assume the detrainment rate for the updrafts to be same as
519 c the entrainment rate at cloud base
520 c
522  do i = 1, im
523  if(cnvflg(i)) then
524  xlamud(i) = xlamue(i,kbcon(i))
525  endif
526  enddo
527 c
528 c determine updraft mass flux for the subcloud layers
529 c
535  do k = km1, 1, -1
536  do i = 1, im
537  if (cnvflg(i)) then
538  if(k.lt.kbcon(i).and.k.ge.kb(i)) then
539  dz = zi(i,k+1) - zi(i,k)
540  ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
541  eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
542  endif
543  endif
544  enddo
545  enddo
546 c
547 c compute mass flux above cloud base
548 c
549  do k = 2, km1
550  do i = 1, im
551  if(cnvflg(i))then
552  if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
553  dz = zi(i,k) - zi(i,k-1)
554  ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
555  eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
556  endif
557  endif
558  enddo
559  enddo
560 c
561 c compute updraft cloud property
562 c
564  do i = 1, im
565  if(cnvflg(i)) then
566  indx = kb(i)
567  hcko(i,indx) = heo(i,indx)
568  ucko(i,indx) = uo(i,indx)
569  vcko(i,indx) = vo(i,indx)
570  endif
571  enddo
572 c
574  do k = 2, km1
575  do i = 1, im
576  if (cnvflg(i)) then
577  if(k.gt.kb(i).and.k.lt.kmax(i)) then
578  dz = zi(i,k) - zi(i,k-1)
579  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
580  tem1 = 0.5 * xlamud(i) * dz
581  factor = 1. + tem - tem1
582  ptem = 0.5 * tem + pgcon
583  ptem1= 0.5 * tem - pgcon
584  hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
585  & (heo(i,k)+heo(i,k-1)))/factor
586  ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)
587  & +ptem1*uo(i,k-1))/factor
588  vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)
589  & +ptem1*vo(i,k-1))/factor
590  dbyo(i,k) = hcko(i,k) - heso(i,k)
591  endif
592  endif
593  enddo
594  enddo
595 c
596 c taking account into convection inhibition due to existence of
597 c dry layers below cloud base
598 c
600  do i=1,im
601  flg(i) = cnvflg(i)
602  kbcon1(i) = kmax(i)
603  enddo
604  do k = 2, km1
605  do i=1,im
606  if (flg(i).and.k.lt.kbm(i)) then
607  if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
608  kbcon1(i) = k
609  flg(i) = .false.
610  endif
611  endif
612  enddo
613  enddo
614  do i=1,im
615  if(cnvflg(i)) then
616  if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
617  endif
618  enddo
619  do i=1,im
620  if(cnvflg(i)) then
621  tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
622  if(tem.gt.dthk) then
623  cnvflg(i) = .false.
624  endif
625  endif
626  enddo
627 !!
628  totflg = .true.
629  do i = 1, im
630  totflg = totflg .and. (.not. cnvflg(i))
631  enddo
632  if(totflg) return
633 !!
634 c
635 c determine first guess cloud top as the level of zero buoyancy
636 c limited to the level of sigma=0.7
637 c
639  do i = 1, im
640  flg(i) = cnvflg(i)
641  if(flg(i)) ktcon(i) = kbm(i)
642  enddo
643  do k = 2, km1
644  do i=1,im
645  if (flg(i).and.k .lt. kbm(i)) then
646  if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
647  ktcon(i) = k
648  flg(i) = .false.
649  endif
650  endif
651  enddo
652  enddo
653 c
654 c turn off shallow convection if cloud top is less than pbl top
655 c
656 ! do i=1,im
657 ! if(cnvflg(i)) then
658 ! kk = kpbl(i)+1
659 ! if(ktcon(i).le.kk) cnvflg(i) = .false.
660 ! endif
661 ! enddo
662 !!
663 ! totflg = .true.
664 ! do i = 1, im
665 ! totflg = totflg .and. (.not. cnvflg(i))
666 ! enddo
667 ! if(totflg) return
668 !!
669 c
670 c specify upper limit of mass flux at cloud base
671 c
673  do i = 1, im
674  if(cnvflg(i)) then
675 ! xmbmax(i) = .1
676 !
677  k = kbcon(i)
678  dp = 1000. * del(i,k)
679  xmbmax(i) = dp / (g * dt2)
680 !
681 ! tem = dp / (g * dt2)
682 ! xmbmax(i) = min(tem, xmbmax(i))
683  endif
684  enddo
685 c
686 c compute cloud moisture property and precipitation
687 c
689  do i = 1, im
690  if (cnvflg(i)) then
691  aa1(i) = 0.
692  qcko(i,kb(i)) = qo(i,kb(i))
693  qrcko(i,kb(i)) = qo(i,kb(i))
694  endif
695  enddo
697  do k = 2, km1
698  do i = 1, im
699  if (cnvflg(i)) then
700  if(k.gt.kb(i).and.k.lt.ktcon(i)) then
701  dz = zi(i,k) - zi(i,k-1)
702  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
703  qrch = qeso(i,k)
704  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
705 cj
706  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
707  tem1 = 0.5 * xlamud(i) * dz
708  factor = 1. + tem - tem1
709  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
710  & (qo(i,k)+qo(i,k-1)))/factor
711  qrcko(i,k) = qcko(i,k)
712 cj
713  dq = eta(i,k) * (qcko(i,k) - qrch)
714 c
715 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
716 c
717 c below lfc check if there is excess moisture to release latent heat
718 c
719  if(k.ge.kbcon(i).and.dq.gt.0.) then
720  etah = .5 * (eta(i,k) + eta(i,k-1))
721  if(ncloud.gt.0.) then
722  dp = 1000. * del(i,k)
723  qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
724  dellal(i,k) = etah * c1 * dz * qlk * g / dp
725  else
726  qlk = dq / (eta(i,k) + etah * c0 * dz)
727  endif
728  aa1(i) = aa1(i) - dz * g * qlk
729  qcko(i,k)= qlk + qrch
730  pwo(i,k) = etah * c0 * dz * qlk
731  cnvwt(i,k) = etah * qlk * g / dp
732  endif
733  endif
734  endif
735  enddo
736  enddo
737 c
738 c calculate cloud work function
739 c
745  do k = 2, km1
746  do i = 1, im
747  if (cnvflg(i)) then
748  if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
749  dz1 = zo(i,k+1) - zo(i,k)
750  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
751  rfact = 1. + delta * cp * gamma
752  & * to(i,k) / hvap
753  aa1(i) = aa1(i) +
754  & dz1 * (g / (cp * to(i,k)))
755  & * dbyo(i,k) / (1. + gamma)
756  & * rfact
757  val = 0.
758  aa1(i)=aa1(i)+
759  & dz1 * g * delta *
760  & max(val,(qeso(i,k) - qo(i,k)))
761  endif
762  endif
763  enddo
764  enddo
766  do i = 1, im
767  if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
768  enddo
769 !!
770  totflg = .true.
771  do i=1,im
772  totflg = totflg .and. (.not. cnvflg(i))
773  enddo
774  if(totflg) return
775 !!
776 c
777 c estimate the onvective overshooting as the level
778 c where the [aafac * cloud work function] becomes zero,
779 c which is the final cloud top
780 c limited to the level of sigma=0.7
781 c
783  do i = 1, im
784  if (cnvflg(i)) then
785  aa1(i) = aafac * aa1(i)
786  endif
787  enddo
788 c
789  do i = 1, im
790  flg(i) = cnvflg(i)
791  ktcon1(i) = kbm(i)
792  enddo
793  do k = 2, km1
794  do i = 1, im
795  if (flg(i)) then
796  if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
797  dz1 = zo(i,k+1) - zo(i,k)
798  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
799  rfact = 1. + delta * cp * gamma
800  & * to(i,k) / hvap
801  aa1(i) = aa1(i) +
802  & dz1 * (g / (cp * to(i,k)))
803  & * dbyo(i,k) / (1. + gamma)
804  & * rfact
805  if(aa1(i).lt.0.) then
806  ktcon1(i) = k
807  flg(i) = .false.
808  endif
809  endif
810  endif
811  enddo
812  enddo
813 c
814 c compute cloud moisture property, detraining cloud water
815 c and precipitation in overshooting layers
816 c
818  do k = 2, km1
819  do i = 1, im
820  if (cnvflg(i)) then
821  if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
822  dz = zi(i,k) - zi(i,k-1)
823  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
824  qrch = qeso(i,k)
825  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
826 cj
827  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
828  tem1 = 0.5 * xlamud(i) * dz
829  factor = 1. + tem - tem1
830  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
831  & (qo(i,k)+qo(i,k-1)))/factor
832  qrcko(i,k) = qcko(i,k)
833 cj
834  dq = eta(i,k) * (qcko(i,k) - qrch)
835 c
836 c check if there is excess moisture to release latent heat
837 c
838  if(dq.gt.0.) then
839  etah = .5 * (eta(i,k) + eta(i,k-1))
840  if(ncloud.gt.0.) then
841  dp = 1000. * del(i,k)
842  qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
843  dellal(i,k) = etah * c1 * dz * qlk * g / dp
844  else
845  qlk = dq / (eta(i,k) + etah * c0 * dz)
846  endif
847  qcko(i,k) = qlk + qrch
848  pwo(i,k) = etah * c0 * dz * qlk
849  cnvwt(i,k) = etah * qlk * g / dp
850  endif
851  endif
852  endif
853  enddo
854  enddo
855 c
856 c exchange ktcon with ktcon1
857 c
858  do i = 1, im
859  if(cnvflg(i)) then
860  kk = ktcon(i)
861  ktcon(i) = ktcon1(i)
862  ktcon1(i) = kk
863  endif
864  enddo
865 c
866 c this section is ready for cloud water
867 c
868  if(ncloud.gt.0) then
869 c
870 c compute liquid and vapor separation at cloud top
871 c
873  do i = 1, im
874  if(cnvflg(i)) then
875  k = ktcon(i) - 1
876  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
877  qrch = qeso(i,k)
878  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
879  dq = qcko(i,k) - qrch
880 c
881 c check if there is excess moisture to release latent heat
882 c
883  if(dq.gt.0.) then
884  qlko_ktcon(i) = dq
885  qcko(i,k) = qrch
886  endif
887  endif
888  enddo
889  endif
890 c
891 c--- compute precipitation efficiency in terms of windshear
892 c
893 !! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 :
894 !! \f[
895 !! E = 1.591 - 0.639\frac{\Delta V}{\Delta z} + 0.0953\left(\frac{\Delta V}{\Delta z}\right)^2 - 0.00496\left(\frac{\Delta V}{\Delta z}\right)^3
896 !! \f]
897 !! where \f$\Delta V\f$ is the integrated horizontal shear over the cloud depth, \f$\Delta z\f$, (the ratio is converted to units of \f$10^{-3} s^{-1}\f$). The variable "edt" is \f$1-E\f$ and is constrained to the range \f$[0,0.9]\f$.
898  do i = 1, im
899  if(cnvflg(i)) then
900  vshear(i) = 0.
901  endif
902  enddo
903  do k = 2, km
904  do i = 1, im
905  if (cnvflg(i)) then
906  if(k.gt.kb(i).and.k.le.ktcon(i)) then
907  shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
908  & + (vo(i,k)-vo(i,k-1)) ** 2)
909  vshear(i) = vshear(i) + shear
910  endif
911  endif
912  enddo
913  enddo
914  do i = 1, im
915  if(cnvflg(i)) then
916  vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
917  e1=1.591-.639*vshear(i)
918  & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
919  edt(i)=1.-e1
920  val = .9
921  edt(i) = min(edt(i),val)
922  val = .0
923  edt(i) = max(edt(i),val)
924  endif
925  enddo
926 c
927 c--- what would the change be, that a cloud with unit mass
928 c--- will do to the environment?
929 c
932  do k = 1, km
933  do i = 1, im
934  if(cnvflg(i) .and. k .le. kmax(i)) then
935  dellah(i,k) = 0.
936  dellaq(i,k) = 0.
937  dellau(i,k) = 0.
938  dellav(i,k) = 0.
939  endif
940  enddo
941  enddo
942 c
943 c--- changed due to subsidence and entrainment
944 c
945  do k = 2, km1
946  do i = 1, im
947  if (cnvflg(i)) then
948  if(k.gt.kb(i).and.k.lt.ktcon(i)) then
949  dp = 1000. * del(i,k)
950  dz = zi(i,k) - zi(i,k-1)
951 c
952  dv1h = heo(i,k)
953  dv2h = .5 * (heo(i,k) + heo(i,k-1))
954  dv3h = heo(i,k-1)
955  dv1q = qo(i,k)
956  dv2q = .5 * (qo(i,k) + qo(i,k-1))
957  dv3q = qo(i,k-1)
958  dv1u = uo(i,k)
959  dv2u = .5 * (uo(i,k) + uo(i,k-1))
960  dv3u = uo(i,k-1)
961  dv1v = vo(i,k)
962  dv2v = .5 * (vo(i,k) + vo(i,k-1))
963  dv3v = vo(i,k-1)
964 c
965  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
966  tem1 = xlamud(i)
967 cj
968  dellah(i,k) = dellah(i,k) +
969  & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h
970  & - tem*eta(i,k-1)*dv2h*dz
971  & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
972  & ) *g/dp
973 cj
974  dellaq(i,k) = dellaq(i,k) +
975  & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q
976  & - tem*eta(i,k-1)*dv2q*dz
977  & + tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz
978  & ) *g/dp
979 cj
980  dellau(i,k) = dellau(i,k) +
981  & ( eta(i,k)*dv1u - eta(i,k-1)*dv3u
982  & - tem*eta(i,k-1)*dv2u*dz
983  & + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz
984  & - pgcon*eta(i,k-1)*(dv1u-dv3u)
985  & ) *g/dp
986 cj
987  dellav(i,k) = dellav(i,k) +
988  & ( eta(i,k)*dv1v - eta(i,k-1)*dv3v
989  & - tem*eta(i,k-1)*dv2v*dz
990  & + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz
991  & - pgcon*eta(i,k-1)*(dv1v-dv3v)
992  & ) *g/dp
993 cj
994  endif
995  endif
996  enddo
997  enddo
998 c
999 c------- cloud top
1000 c
1001  do i = 1, im
1002  if(cnvflg(i)) then
1003  indx = ktcon(i)
1004  dp = 1000. * del(i,indx)
1005  dv1h = heo(i,indx-1)
1006  dellah(i,indx) = eta(i,indx-1) *
1007  & (hcko(i,indx-1) - dv1h) * g / dp
1008  dv1q = qo(i,indx-1)
1009  dellaq(i,indx) = eta(i,indx-1) *
1010  & (qcko(i,indx-1) - dv1q) * g / dp
1011  dv1u = uo(i,indx-1)
1012  dellau(i,indx) = eta(i,indx-1) *
1013  & (ucko(i,indx-1) - dv1u) * g / dp
1014  dv1v = vo(i,indx-1)
1015  dellav(i,indx) = eta(i,indx-1) *
1016  & (vcko(i,indx-1) - dv1v) * g / dp
1017 c
1018 c cloud water
1019 c
1020  dellal(i,indx) = eta(i,indx-1) *
1021  & qlko_ktcon(i) * g / dp
1022  endif
1023  enddo
1024 c
1025 c mass flux at cloud base for shallow convection
1026 c (grant, 2001)
1027 c
1033  do i= 1, im
1034  if(cnvflg(i)) then
1035  k = kbcon(i)
1036 ! ptem = g*sflx(i)*zi(i,k)/t1(i,1)
1037  ptem = g*sflx(i)*hpbl(i)/t1(i,1)
1038  wstar(i) = ptem**h1
1039  tem = po(i,k)*100. / (rd*t1(i,k))
1040  xmb(i) = betaw*tem*wstar(i)
1041  xmb(i) = min(xmb(i),xmbmax(i))
1042  endif
1043  enddo
1046 c
1047 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1048 c
1049  do k = 1, km
1050  do i = 1, im
1051  if (cnvflg(i) .and. k .le. kmax(i)) then
1052  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1053  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
1054  val = 1.e-8
1055  qeso(i,k) = max(qeso(i,k), val )
1056  endif
1057  enddo
1058  enddo
1059 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1060 c
1064  do i = 1, im
1065  delhbar(i) = 0.
1066  delqbar(i) = 0.
1067  deltbar(i) = 0.
1068  delubar(i) = 0.
1069  delvbar(i) = 0.
1070  qcond(i) = 0.
1071  enddo
1072  do k = 1, km
1073  do i = 1, im
1074  if (cnvflg(i)) then
1075  if(k.gt.kb(i).and.k.le.ktcon(i)) then
1076  dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1077  t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1078  q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1079 ! tem = 1./rcs(i)
1080 ! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1081 ! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1082  u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
1083  v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
1084  dp = 1000. * del(i,k)
1085  delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
1086  delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
1087  deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
1088  delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
1089  delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
1090  endif
1091  endif
1092  enddo
1093  enddo
1095  do k = 1, km
1096  do i = 1, im
1097  if (cnvflg(i)) then
1098  if(k.gt.kb(i).and.k.le.ktcon(i)) then
1099  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1100  qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
1101  val = 1.e-8
1102  qeso(i,k) = max(qeso(i,k), val )
1103  endif
1104  endif
1105  enddo
1106  enddo
1107 c
1109  do i = 1, im
1110  rntot(i) = 0.
1111  delqev(i) = 0.
1112  delq2(i) = 0.
1113  flg(i) = cnvflg(i)
1114  enddo
1115  do k = km, 1, -1
1116  do i = 1, im
1117  if (cnvflg(i)) then
1118  if(k.lt.ktcon(i).and.k.gt.kb(i)) then
1119  rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
1120  endif
1121  endif
1122  enddo
1123  enddo
1124 c
1125 c evaporating rain
1126 c
1130  do k = km, 1, -1
1131  do i = 1, im
1132  if (k .le. kmax(i)) then
1133  deltv(i) = 0.
1134  delq(i) = 0.
1135  qevap(i) = 0.
1136  if(cnvflg(i)) then
1137  if(k.lt.ktcon(i).and.k.gt.kb(i)) then
1138  rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
1139  endif
1140  endif
1141  if(flg(i).and.k.lt.ktcon(i)) then
1142  evef = edt(i) * evfact
1143  if(islimsk(i) == 1) evef=edt(i) * evfactl
1144 ! if(islimsk(i) == 1) evef=.07
1145 c if(islimsk(i) == 1) evef = 0.
1146  qcond(i) = evef * (q1(i,k) - qeso(i,k))
1147  & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
1148  dp = 1000. * del(i,k)
1149  if(rn(i).gt.0..and.qcond(i).lt.0.) then
1150  qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
1151  qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
1152  delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
1153  endif
1154  if(rn(i).gt.0..and.qcond(i).lt.0..and.
1155  & delq2(i).gt.rntot(i)) then
1156  qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
1157  flg(i) = .false.
1158  endif
1159  if(rn(i).gt.0..and.qevap(i).gt.0.) then
1160  tem = .001 * dp / g
1161  tem1 = qevap(i) * tem
1162  if(tem1.gt.rn(i)) then
1163  qevap(i) = rn(i) / tem
1164  rn(i) = 0.
1165  else
1166  rn(i) = rn(i) - tem1
1167  endif
1168  q1(i,k) = q1(i,k) + qevap(i)
1169  t1(i,k) = t1(i,k) - elocp * qevap(i)
1170  deltv(i) = - elocp*qevap(i)/dt2
1171  delq(i) = + qevap(i)/dt2
1172  delqev(i) = delqev(i) + .001*dp*qevap(i)/g
1173  endif
1174  dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
1175  delqbar(i) = delqbar(i) + delq(i)*dp/g
1176  deltbar(i) = deltbar(i) + deltv(i)*dp/g
1177  endif
1178  endif
1179  enddo
1180  enddo
1181 cj
1182 ! do i = 1, im
1183 ! if(me.eq.31.and.cnvflg(i)) then
1184 ! if(cnvflg(i)) then
1185 ! print *, ' shallow delhbar, delqbar, deltbar = ',
1186 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
1187 ! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
1188 ! print *, ' precip =', hvap*rn(i)*1000./dt2
1189 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
1190 ! endif
1191 ! enddo
1192 cj
1193  do i = 1, im
1194  if(cnvflg(i)) then
1195  if(rn(i).lt.0..or..not.flg(i)) rn(i) = 0.
1196  ktop(i) = ktcon(i)
1197  kbot(i) = kbcon(i)
1198  kcnv(i) = 0
1199  endif
1200  enddo
1201 c
1202 c convective cloud water
1203 c
1205  do k = 1, km
1206  do i = 1, im
1207  if (cnvflg(i) .and. rn(i).gt.0.) then
1208  if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1209  cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
1210  endif
1211  endif
1212  enddo
1213  enddo
1214 
1215 c
1216 c convective cloud cover
1217 c
1219  do k = 1, km
1220  do i = 1, im
1221  if (cnvflg(i) .and. rn(i).gt.0.) then
1222  if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1223  cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
1224  cnvc(i,k) = min(cnvc(i,k), 0.2)
1225  cnvc(i,k) = max(cnvc(i,k), 0.0)
1226  endif
1227  endif
1228  enddo
1229  enddo
1230 
1231 c
1232 c cloud water
1233 c
1235  if (ncloud.gt.0) then
1236 !
1237  do k = 1, km1
1238  do i = 1, im
1239  if (cnvflg(i)) then
1240  if (k.gt.kb(i).and.k.le.ktcon(i)) then
1241  tem = dellal(i,k) * xmb(i) * dt2
1242  tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
1243  if (ql(i,k,2) .gt. -999.0) then
1244  ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
1245  ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
1246  else
1247  ql(i,k,1) = ql(i,k,1) + tem
1248  endif
1249  endif
1250  endif
1251  enddo
1252  enddo
1253 !
1254  endif
1255 !
1256 ! hchuang code change
1257 !
1259  do k = 1, km
1260  do i = 1, im
1261  if(cnvflg(i)) then
1262  if(k.ge.kb(i) .and. k.lt.ktop(i)) then
1263  ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
1264  endif
1265  endif
1266  enddo
1267  enddo
1269  do i = 1, im
1270  if(cnvflg(i)) then
1271  k = ktop(i)-1
1272  dt_mf(i,k) = ud_mf(i,k)
1273  endif
1274  enddo
1275 !!
1276  return
1277  end
1278 
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
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
shalcnv
subroutine shalcnv(im, ix, km, jcap, delt, delp, prslp, psp, phil, ql, q1, t1, u1, v1, rn, kbot, ktop, kcnv, islimsk, dot, ncloud, hpbl, heat, evap, ud_mf, dt_mf, cnvw, cnvc)
This subroutine contains the entirety of the SAS shallow convection scheme.
Definition: shalcnv.f:58
physcons::con_rd
real(kind=kind_phys), parameter con_rd
gas constant air ( )
Definition: physcons.f:76