GFS Physics Documentation
samfshalcnv.f
Go to the documentation of this file.
1 
11 
14 
61  subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql,
62  & q1,t1,u1,v1,rn,kbot,ktop,kcnv,islimsk,garea,
63  & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc,
64 ! & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc,me)
65  & clam,c0s,c1,pgcon,asolfac)
66 !
67  use machine , only : kind_phys
68  use funcphys , only : fpvs
69  use physcons, grav => con_g, cp => con_cp, hvap => con_hvap
70  &, rv => con_rv, fv => con_fvirt, t0c => con_t0c
71  &, rd => con_rd, cvap => con_cvap, cliq => con_cliq
72  &, eps => con_eps, epsm1 => con_epsm1
73  implicit none
74 !
75  integer im, ix, km, ncloud,
76  & kbot(im), ktop(im), kcnv(im)
77 ! &, me
78  real(kind=kind_phys) delt
79  real(kind=kind_phys) psp(im), delp(ix,km), prslp(ix,km)
80  real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km),
81  & ql(ix,km,2),q1(ix,km), t1(ix,km),
82  & u1(ix,km), v1(ix,km),
83 ! & u1(ix,km), v1(ix,km), rcs(im),
84  & rn(im), garea(im),
85  & dot(ix,km), phil(ix,km), hpbl(im),
86  & cnvw(ix,km),cnvc(ix,km)
87 ! hchuang code change mass flux output
88  &, ud_mf(im,km),dt_mf(im,km)
89 !
90  integer i,j,indx, k, kk, km1, n
91  integer kpbl(im)
92  integer, dimension(im), intent(in) :: islimsk
93 !
94  real(kind=kind_phys) dellat, delta,
95  & c0l, c0s, d0,
96  & c1, asolfac,
97  & desdt, dp,
98  & dq, dqsdp, dqsdt, dt,
99  & dt2, dtmax, dtmin, dxcrt,
100  & dv1h, dv2h, dv3h,
101  & dv1q, dv2q, dv3q,
102  & dz, dz1, e1, clam,
103  & el2orc, elocp, aafac, cm,
104  & es, etah, h1,
105  & evef, evfact, evfactl, fact1,
106  & fact2, factor, dthk,
107  & g, gamma, pprime, betaw,
108  & qlk, qrch, qs,
109  & rfact, shear, tfac,
110  & val, val1, val2,
111  & w1, w1l, w1s, w2,
112  & w2l, w2s, w3, w3l,
113  & w3s, w4, w4l, w4s,
114  & rho, tem, tem1, tem2,
115  & ptem, ptem1,
116  & pgcon
117 !
118  integer kb(im), kbcon(im), kbcon1(im),
119  & ktcon(im), ktcon1(im), ktconn(im),
120  & kbm(im), kmax(im)
121 !
122  real(kind=kind_phys) aa1(im), cina(im),
123  & umean(im), tauadv(im), gdx(im),
124  & delhbar(im), delq(im), delq2(im),
125  & delqbar(im), delqev(im), deltbar(im),
126  & deltv(im), dtconv(im), edt(im),
127  & pdot(im), po(im,km),
128  & qcond(im), qevap(im), hmax(im),
129  & rntot(im), vshear(im),
130  & xlamud(im), xmb(im), xmbmax(im),
131  & delubar(im), delvbar(im)
132 !
133  real(kind=kind_phys) c0(im)
134 c
135  real(kind=kind_phys) crtlamd
136 !
137  real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn,
138  & cinacr, cinacrmx, cinacrmn
139 !
140 ! parameters for updraft velocity calculation
141  real(kind=kind_phys) bet1, cd1, f1, gam1,
142  & bb1, bb2
143 ! & bb1, bb2, wucb
144 cc
145 c physical parameters
146 ! parameter(g=grav,asolfac=0.89)
147  parameter(g=grav)
148  parameter(elocp=hvap/cp,
149  & el2orc=hvap*hvap/(rv*cp))
150 ! parameter(c0s=0.002,c1=5.e-4,d0=.01)
151  parameter(d0=.01)
152 ! parameter(c0l=c0s*asolfac)
153 !
154 ! asolfac: aerosol-aware parameter based on Lim & Hong (2012)
155 ! asolfac= cx / c0s(=.002)
156 ! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s)
157 ! Nccn: CCN number concentration in cm^(-3)
158 ! Until a realistic Nccn is provided, Nccns are assumed
159 ! as Nccn=100 for sea and Nccn=7000 for land
160 !
161  parameter(cm=1.0,delta=fv)
162  parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
163  parameter(dthk=25.)
164  parameter(cinpcrmx=180.,cinpcrmn=120.)
165 ! parameter(cinacrmx=-120.,cinacrmn=-120.)
166  parameter(cinacrmx=-120.,cinacrmn=-80.)
167  parameter(crtlamd=3.e-4)
168  parameter(dtmax=10800.,dtmin=600.)
169  parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
170  parameter(betaw=.03,dxcrt=15.e3)
171  parameter(h1=0.33333333)
172 c local variables and arrays
173  real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
174  & uo(im,km), vo(im,km), qeso(im,km)
175 ! for updraft velocity calculation
176  real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km)
177  real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im)
178 !
179 c cloud water
180 ! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
181  real(kind=kind_phys) qlko_ktcon(im), dellal(im,km),
182  & dbyo(im,km), zo(im,km), xlamue(im,km),
183  & heo(im,km), heso(im,km),
184  & dellah(im,km), dellaq(im,km),
185  & dellau(im,km), dellav(im,km), hcko(im,km),
186  & ucko(im,km), vcko(im,km), qcko(im,km),
187  & qrcko(im,km), eta(im,km),
188  & zi(im,km), pwo(im,km), c0t(im,km),
189  & sumx(im), tx1(im), cnvwt(im,km)
190 !
191  logical totflg, cnvflg(im), flg(im)
192 !
193  real(kind=kind_phys) tf, tcr, tcrf
194  parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
195 !
196 c-----------------------------------------------------------------------
197 !
198 !************************************************************************
199 ! convert input Pa terms to Cb terms -- Moorthi
202  ps = psp * 0.001
203  prsl = prslp * 0.001
204  del = delp * 0.001
205 !************************************************************************
206 !
207  km1 = km - 1
208 c
209 c initialize arrays
210 c
212  do i=1,im
213  cnvflg(i) = .true.
214  if(kcnv(i) == 1) cnvflg(i) = .false.
215  if(cnvflg(i)) then
216  kbot(i)=km+1
217  ktop(i)=0
218  endif
219  rn(i)=0.
220  kbcon(i)=km
221  ktcon(i)=1
222  ktconn(i)=1
223  kb(i)=km
224  pdot(i) = 0.
225  qlko_ktcon(i) = 0.
226  edt(i) = 0.
227  aa1(i) = 0.
228  cina(i) = 0.
229  vshear(i) = 0.
230  gdx(i) = sqrt(garea(i))
231  enddo
232 !!
234  totflg = .true.
235  do i=1,im
236  totflg = totflg .and. (.not. cnvflg(i))
237  enddo
238  if(totflg) return
239 !!
241  do i=1,im
242  if(islimsk(i) == 1) then
243  c0(i) = c0s*asolfac
244  else
245  c0(i) = c0s
246  endif
247  enddo
248 !
250  do k = 1, km
251  do i = 1, im
252  if(t1(i,k) > 273.16) then
253  c0t(i,k) = c0(i)
254  else
255  tem = d0 * (t1(i,k) - 273.16)
256  tem1 = exp(tem)
257  c0t(i,k) = c0(i) * tem1
258  endif
259  enddo
260  enddo
261 !
263  do k = 1, km
264  do i = 1, im
265  cnvw(i,k) = 0.
266  cnvc(i,k) = 0.
267  enddo
268  enddo
269 ! hchuang code change
271  do k = 1, km
272  do i = 1, im
273  ud_mf(i,k) = 0.
274  dt_mf(i,k) = 0.
275  enddo
276  enddo
277 c
278  dt2 = delt
279 !
280 c model tunable parameters are all here
281 ! clam = .3
282  aafac = .1
283 c evef = 0.07
284  evfact = 0.3
285  evfactl = 0.3
286 !
287 ! pgcon = 0.7 ! Gregory et al. (1997, QJRMS)
288 ! pgcon = 0.55 ! Zhang & Wu (2003,JAS)
289  w1l = -8.e-3
290  w2l = -4.e-2
291  w3l = -5.e-3
292  w4l = -5.e-4
293  w1s = -2.e-4
294  w2s = -2.e-3
295  w3s = -1.e-3
296  w4s = -2.e-5
297 c
298 c define top layer for search of the downdraft originating layer
299 c and the maximum thetae for updraft
300 c
302  do i=1,im
303  kbm(i) = km
304  kmax(i) = km
305  tx1(i) = 1.0 / ps(i)
306  enddo
307 !
308  do k = 1, km
309  do i=1,im
310  if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1
311  if (prsl(i,k)*tx1(i) > 0.60) kmax(i) = k + 1
312  enddo
313  enddo
314  do i=1,im
315  kbm(i) = min(kbm(i),kmax(i))
316  enddo
317 c
318 c hydrostatic height assume zero terr and compute
319 c updraft entrainment rate as an inverse function of height
320 c
322  do k = 1, km
323  do i=1,im
324  zo(i,k) = phil(i,k) / g
325  enddo
326  enddo
328  do k = 1, km1
329  do i=1,im
330  zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
331  xlamue(i,k) = clam / zi(i,k)
332  enddo
333  enddo
334  do i=1,im
335  xlamue(i,km) = xlamue(i,km1)
336  enddo
337 c
338 c pbl height
339 c
341  do i=1,im
342  flg(i) = cnvflg(i)
343  kpbl(i)= 1
344  enddo
345  do k = 2, km1
346  do i=1,im
347  if (flg(i) .and. zo(i,k) <= hpbl(i)) then
348  kpbl(i) = k
349  else
350  flg(i) = .false.
351  endif
352  enddo
353  enddo
354  do i=1,im
355  kpbl(i)= min(kpbl(i),kbm(i))
356  enddo
357 c
358 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
359 c convert surface pressure to mb from cb
360 c
362  do k = 1, km
363  do i = 1, im
364  if (cnvflg(i) .and. k <= kmax(i)) then
365  pfld(i,k) = prsl(i,k) * 10.0
366  eta(i,k) = 1.
367  hcko(i,k) = 0.
368  qcko(i,k) = 0.
369  qrcko(i,k)= 0.
370  ucko(i,k) = 0.
371  vcko(i,k) = 0.
372  dbyo(i,k) = 0.
373  pwo(i,k) = 0.
374  dellal(i,k) = 0.
375  to(i,k) = t1(i,k)
376  qo(i,k) = q1(i,k)
377  uo(i,k) = u1(i,k)
378  vo(i,k) = v1(i,k)
379 ! uo(i,k) = u1(i,k) * rcs(i)
380 ! vo(i,k) = v1(i,k) * rcs(i)
381  wu2(i,k) = 0.
382  buo(i,k) = 0.
383  drag(i,k) = 0.
384  cnvwt(i,k) = 0.
385  endif
386  enddo
387  enddo
389  do k = 1, km
390  do i=1,im
391  if (cnvflg(i) .and. k <= kmax(i)) then
392  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
393  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
394  val1 = 1.e-8
395  qeso(i,k) = max(qeso(i,k), val1)
396  val2 = 1.e-10
397  qo(i,k) = max(qo(i,k), val2 )
398 ! qo(i,k) = min(qo(i,k),qeso(i,k))
399 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
400  endif
401  enddo
402  enddo
403 c
404 c compute moist static energy
405 c
407  do k = 1, km
408  do i=1,im
409  if (cnvflg(i) .and. k <= kmax(i)) then
410 ! tem = g * zo(i,k) + cp * to(i,k)
411  tem = phil(i,k) + cp * to(i,k)
412  heo(i,k) = tem + hvap * qo(i,k)
413  heso(i,k) = tem + hvap * qeso(i,k)
414 c heo(i,k) = min(heo(i,k),heso(i,k))
415  endif
416  enddo
417  enddo
418 c
419 c determine level with largest moist static energy within pbl
420 c this is the level where updraft starts
421 c
424  do i=1,im
425  if (cnvflg(i)) then
426  hmax(i) = heo(i,1)
427  kb(i) = 1
428  endif
429  enddo
430  do k = 2, km
431  do i=1,im
432  if (cnvflg(i) .and. k <= kpbl(i)) then
433  if(heo(i,k) > hmax(i)) then
434  kb(i) = k
435  hmax(i) = heo(i,k)
436  endif
437  endif
438  enddo
439  enddo
440 c
442  do k = 1, km1
443  do i=1,im
444  if (cnvflg(i) .and. k <= kmax(i)-1) then
445  dz = .5 * (zo(i,k+1) - zo(i,k))
446  dp = .5 * (pfld(i,k+1) - pfld(i,k))
447  es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
448  pprime = pfld(i,k+1) + epsm1 * es
449  qs = eps * es / pprime
450  dqsdp = - qs / pprime
451  desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
452  dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
453  gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
454  dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
455  dq = dqsdt * dt + dqsdp * dp
456  to(i,k) = to(i,k+1) + dt
457  qo(i,k) = qo(i,k+1) + dq
458  po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
459  endif
460  enddo
461  enddo
462 !
464  do k = 1, km1
465  do i=1,im
466  if (cnvflg(i) .and. k <= kmax(i)-1) then
467  qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
468  qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
469  val1 = 1.e-8
470  qeso(i,k) = max(qeso(i,k), val1)
471  val2 = 1.e-10
472  qo(i,k) = max(qo(i,k), val2 )
473 ! qo(i,k) = min(qo(i,k),qeso(i,k))
474  heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
475  & cp * to(i,k) + hvap * qo(i,k)
476  heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
477  & cp * to(i,k) + hvap * qeso(i,k)
478  uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
479  vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
480  endif
481  enddo
482  enddo
483 c
484 c look for the level of free convection as cloud base
485 c
487  do i=1,im
488  flg(i) = cnvflg(i)
489  if(flg(i)) kbcon(i) = kmax(i)
490  enddo
491  do k = 2, km1
492  do i=1,im
493  if (flg(i) .and. k < kbm(i)) then
494  if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
495  kbcon(i) = k
496  flg(i) = .false.
497  endif
498  endif
499  enddo
500  enddo
501 c
502  do i=1,im
503  if(cnvflg(i)) then
504  if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
505  endif
506  enddo
507 !!
509  totflg = .true.
510  do i=1,im
511  totflg = totflg .and. (.not. cnvflg(i))
512  enddo
513  if(totflg) return
514 !!
516  do i=1,im
517  if(cnvflg(i)) then
518 ! pdot(i) = 10.* dot(i,kbcon(i))
519  pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
520  endif
521  enddo
522 c
523 c turn off convection if pressure depth between parcel source level
524 c and cloud base is larger than a critical value, cinpcr
525 c
526  do i=1,im
527  if(cnvflg(i)) then
528  if(islimsk(i) == 1) then
529  w1 = w1l
530  w2 = w2l
531  w3 = w3l
532  w4 = w4l
533  else
534  w1 = w1s
535  w2 = w2s
536  w3 = w3s
537  w4 = w4s
538  endif
539  if(pdot(i) <= w4) then
540  tem = (pdot(i) - w4) / (w3 - w4)
541  elseif(pdot(i) >= -w4) then
542  tem = - (pdot(i) + w4) / (w4 - w3)
543  else
544  tem = 0.
545  endif
546  val1 = -1.
547  tem = max(tem,val1)
548  val2 = 1.
549  tem = min(tem,val2)
550  ptem = 1. - tem
551  ptem1= .5*(cinpcrmx-cinpcrmn)
552  cinpcr = cinpcrmx - ptem * ptem1
553  tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
554  if(tem1 > cinpcr) then
555  cnvflg(i) = .false.
556  endif
557  endif
558  enddo
559 !!
560  totflg = .true.
561  do i=1,im
562  totflg = totflg .and. (.not. cnvflg(i))
563  enddo
564  if(totflg) return
565 !!
566 c
567 c specify the detrainment rate for the updrafts
568 c
570  do i = 1, im
571  if(cnvflg(i)) then
572  xlamud(i) = xlamue(i,kbcon(i))
573 ! xlamud(i) = crtlamd
574  endif
575  enddo
576 c
577 c determine updraft mass flux for the subcloud layers
578 c
584  do k = km1, 1, -1
585  do i = 1, im
586  if (cnvflg(i)) then
587  if(k < kbcon(i) .and. k >= kb(i)) then
588  dz = zi(i,k+1) - zi(i,k)
589  ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
590  eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
591  endif
592  endif
593  enddo
594  enddo
595 c
596 c compute mass flux above cloud base
597 c
598  do i = 1, im
599  flg(i) = cnvflg(i)
600  enddo
601  do k = 2, km1
602  do i = 1, im
603  if(flg(i))then
604  if(k > kbcon(i) .and. k < kmax(i)) then
605  dz = zi(i,k) - zi(i,k-1)
606  ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
607  eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
608  if(eta(i,k) <= 0.) then
609  kmax(i) = k
610  ktconn(i) = k
611  kbm(i) = min(kbm(i),kmax(i))
612  flg(i) = .false.
613  endif
614  endif
615  endif
616  enddo
617  enddo
618 c
619 c compute updraft cloud property
620 c
622  do i = 1, im
623  if(cnvflg(i)) then
624  indx = kb(i)
625  hcko(i,indx) = heo(i,indx)
626  ucko(i,indx) = uo(i,indx)
627  vcko(i,indx) = vo(i,indx)
628  endif
629  enddo
630 c
631 ! cm is an enhancement factor in entrainment rates for momentum
632 !
634  do k = 2, km1
635  do i = 1, im
636  if (cnvflg(i)) then
637  if(k > kb(i) .and. k < kmax(i)) then
638  dz = zi(i,k) - zi(i,k-1)
639  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
640  tem1 = 0.5 * xlamud(i) * dz
641  factor = 1. + tem - tem1
642  hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
643  & (heo(i,k)+heo(i,k-1)))/factor
644  dbyo(i,k) = hcko(i,k) - heso(i,k)
645 !
646  tem = 0.5 * cm * tem
647  factor = 1. + tem
648  ptem = tem + pgcon
649  ptem1= tem - pgcon
650  ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*uo(i,k)
651  & +ptem1*uo(i,k-1))/factor
652  vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*vo(i,k)
653  & +ptem1*vo(i,k-1))/factor
654  endif
655  endif
656  enddo
657  enddo
658 c
659 c taking account into convection inhibition due to existence of
660 c dry layers below cloud base
661 c
663  do i=1,im
664  flg(i) = cnvflg(i)
665  kbcon1(i) = kmax(i)
666  enddo
667  do k = 2, km1
668  do i=1,im
669  if (flg(i) .and. k < kbm(i)) then
670  if(k >= kbcon(i) .and. dbyo(i,k) > 0.) then
671  kbcon1(i) = k
672  flg(i) = .false.
673  endif
674  endif
675  enddo
676  enddo
677  do i=1,im
678  if(cnvflg(i)) then
679  if(kbcon1(i) == kmax(i)) cnvflg(i) = .false.
680  endif
681  enddo
682  do i=1,im
683  if(cnvflg(i)) then
684  tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
685  if(tem > dthk) then
686  cnvflg(i) = .false.
687  endif
688  endif
689  enddo
690 !!
691  totflg = .true.
692  do i = 1, im
693  totflg = totflg .and. (.not. cnvflg(i))
694  enddo
695  if(totflg) return
696 !!
697 c
698 c calculate convective inhibition
699 c
701  do k = 2, km1
702  do i = 1, im
703  if (cnvflg(i)) then
704  if(k > kb(i) .and. k < kbcon1(i)) then
705  dz1 = zo(i,k+1) - zo(i,k)
706  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
707  rfact = 1. + delta * cp * gamma
708  & * to(i,k) / hvap
709  cina(i) = cina(i) +
710 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
711  & dz1 * (g / (cp * to(i,k)))
712  & * dbyo(i,k) / (1. + gamma)
713  & * rfact
714  val = 0.
715  cina(i) = cina(i) +
716 ! & dz1 * eta(i,k) * g * delta *
717  & dz1 * g * delta *
718  & max(val,(qeso(i,k) - qo(i,k)))
719  endif
720  endif
721  enddo
722  enddo
724  do i = 1, im
725  if(cnvflg(i)) then
726 !
727  if(islimsk(i) == 1) then
728  w1 = w1l
729  w2 = w2l
730  w3 = w3l
731  w4 = w4l
732  else
733  w1 = w1s
734  w2 = w2s
735  w3 = w3s
736  w4 = w4s
737  endif
738  if(pdot(i) <= w4) then
739  tem = (pdot(i) - w4) / (w3 - w4)
740  elseif(pdot(i) >= -w4) then
741  tem = - (pdot(i) + w4) / (w4 - w3)
742  else
743  tem = 0.
744  endif
745 
746  val1 = -1.
747  tem = max(tem,val1)
748  val2 = 1.
749  tem = min(tem,val2)
750  tem = 1. - tem
751  tem1= .5*(cinacrmx-cinacrmn)
752  cinacr = cinacrmx - tem * tem1
753 !
754 ! cinacr = cinacrmx
755  if(cina(i) < cinacr) cnvflg(i) = .false.
756  endif
757  enddo
758 !!
759  totflg = .true.
760  do i=1,im
761  totflg = totflg .and. (.not. cnvflg(i))
762  enddo
763  if(totflg) return
764 !!
765 c
766 c determine first guess cloud top as the level of zero buoyancy
767 c limited to the level of p/ps=0.7
768 c
770  do i = 1, im
771  flg(i) = cnvflg(i)
772  if(flg(i)) ktcon(i) = kbm(i)
773  enddo
774  do k = 2, km1
775  do i=1,im
776  if (flg(i) .and. k < kbm(i)) then
777  if(k > kbcon1(i) .and. dbyo(i,k) < 0.) then
778  ktcon(i) = k
779  flg(i) = .false.
780  endif
781  endif
782  enddo
783  enddo
784 c
785 c specify upper limit of mass flux at cloud base
786 c
788  do i = 1, im
789  if(cnvflg(i)) then
790 ! xmbmax(i) = .1
791 !
792  k = kbcon(i)
793  dp = 1000. * del(i,k)
794  xmbmax(i) = dp / (g * dt2)
795 !
796 ! tem = dp / (g * dt2)
797 ! xmbmax(i) = min(tem, xmbmax(i))
798  endif
799  enddo
800 c
801 c compute cloud moisture property and precipitation
802 c
804  do i = 1, im
805  if (cnvflg(i)) then
806  aa1(i) = 0.
807  qcko(i,kb(i)) = qo(i,kb(i))
808  qrcko(i,kb(i)) = qo(i,kb(i))
809  endif
810  enddo
812  do k = 2, km1
813  do i = 1, im
814  if (cnvflg(i)) then
815  if(k > kb(i) .and. k < ktcon(i)) then
816  dz = zi(i,k) - zi(i,k-1)
817  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
818  qrch = qeso(i,k)
819  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
820 cj
821  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
822  tem1 = 0.5 * xlamud(i) * dz
823  factor = 1. + tem - tem1
824  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
825  & (qo(i,k)+qo(i,k-1)))/factor
826  qrcko(i,k) = qcko(i,k)
827 cj
828  dq = eta(i,k) * (qcko(i,k) - qrch)
829 c
830 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
831 c
832 c below lfc check if there is excess moisture to release latent heat
833 c
834  if(k >= kbcon(i) .and. dq > 0.) then
835  etah = .5 * (eta(i,k) + eta(i,k-1))
836  dp = 1000. * del(i,k)
837  if(ncloud > 0) then
838  ptem = c0t(i,k) + c1
839  qlk = dq / (eta(i,k) + etah * ptem * dz)
840  dellal(i,k) = etah * c1 * dz * qlk * g / dp
841  else
842  qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
843  endif
844  buo(i,k) = buo(i,k) - g * qlk
845  qcko(i,k)= qlk + qrch
846  pwo(i,k) = etah * c0t(i,k) * dz * qlk
847  cnvwt(i,k) = etah * qlk * g / dp
848  endif
849 !
850 ! compute buoyancy and drag for updraft velocity
851 !
852  if(k >= kbcon(i)) then
853  rfact = 1. + delta * cp * gamma
854  & * to(i,k) / hvap
855  buo(i,k) = buo(i,k) + (g / (cp * to(i,k)))
856  & * dbyo(i,k) / (1. + gamma)
857  & * rfact
858  val = 0.
859  buo(i,k) = buo(i,k) + g * delta *
860  & max(val,(qeso(i,k) - qo(i,k)))
861  drag(i,k) = max(xlamue(i,k),xlamud(i))
862  endif
863 !
864  endif
865  endif
866  enddo
867  enddo
868 c
869 c calculate cloud work function
870 c
871 ! do k = 2, km1
872 ! do i = 1, im
873 ! if (cnvflg(i)) then
874 ! if(k >= kbcon(i) .and. k < ktcon(i)) then
875 ! dz1 = zo(i,k+1) - zo(i,k)
876 ! gamma = el2orc * qeso(i,k) / (to(i,k)**2)
877 ! rfact = 1. + delta * cp * gamma
878 ! & * to(i,k) / hvap
879 ! aa1(i) = aa1(i) +
880 !! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
881 ! & dz1 * (g / (cp * to(i,k)))
882 ! & * dbyo(i,k) / (1. + gamma)
883 ! & * rfact
884 ! val = 0.
885 ! aa1(i) = aa1(i) +
886 !! & dz1 * eta(i,k) * g * delta *
887 ! & dz1 * g * delta *
888 ! & max(val,(qeso(i,k) - qo(i,k)))
889 ! endif
890 ! endif
891 ! enddo
892 ! enddo
893 ! do i = 1, im
894 ! if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
895 ! enddo
896 !
897 ! calculate cloud work function
898 !
904  do i = 1, im
905  if (cnvflg(i)) then
906  aa1(i) = 0.
907  endif
908  enddo
909  do k = 2, km1
910  do i = 1, im
911  if (cnvflg(i)) then
912  if(k >= kbcon(i) .and. k < ktcon(i)) then
913  dz1 = zo(i,k+1) - zo(i,k)
914  aa1(i) = aa1(i) + buo(i,k) * dz1
915  endif
916  endif
917  enddo
918  enddo
919  do i = 1, im
920  if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
921  enddo
922 !!
924  totflg = .true.
925  do i=1,im
926  totflg = totflg .and. (.not. cnvflg(i))
927  enddo
928  if(totflg) return
929 !!
930 c
931 c estimate the onvective overshooting as the level
932 c where the [aafac * cloud work function] becomes zero,
933 c which is the final cloud top
934 c limited to the level of p/ps=0.7
935 c
937  do i = 1, im
938  if (cnvflg(i)) then
939  aa1(i) = aafac * aa1(i)
940  endif
941  enddo
942 c
943  do i = 1, im
944  flg(i) = cnvflg(i)
945  ktcon1(i) = kbm(i)
946  enddo
947  do k = 2, km1
948  do i = 1, im
949  if (flg(i)) then
950  if(k >= ktcon(i) .and. k < kbm(i)) then
951  dz1 = zo(i,k+1) - zo(i,k)
952  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
953  rfact = 1. + delta * cp * gamma
954  & * to(i,k) / hvap
955  aa1(i) = aa1(i) +
956 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
957  & dz1 * (g / (cp * to(i,k)))
958  & * dbyo(i,k) / (1. + gamma)
959  & * rfact
960 ! val = 0.
961 ! aa1(i) = aa1(i) +
962 !! & dz1 * eta(i,k) * g * delta *
963 ! & dz1 * g * delta *
964 ! & max(val,(qeso(i,k) - qo(i,k)))
965  if(aa1(i) < 0.) then
966  ktcon1(i) = k
967  flg(i) = .false.
968  endif
969  endif
970  endif
971  enddo
972  enddo
973 c
974 c compute cloud moisture property, detraining cloud water
975 c and precipitation in overshooting layers
976 c
978  do k = 2, km1
979  do i = 1, im
980  if (cnvflg(i)) then
981  if(k >= ktcon(i) .and. k < ktcon1(i)) then
982  dz = zi(i,k) - zi(i,k-1)
983  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
984  qrch = qeso(i,k)
985  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
986 cj
987  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
988  tem1 = 0.5 * xlamud(i) * dz
989  factor = 1. + tem - tem1
990  qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
991  & (qo(i,k)+qo(i,k-1)))/factor
992  qrcko(i,k) = qcko(i,k)
993 cj
994  dq = eta(i,k) * (qcko(i,k) - qrch)
995 c
996 c check if there is excess moisture to release latent heat
997 c
998  if(dq > 0.) then
999  etah = .5 * (eta(i,k) + eta(i,k-1))
1000  dp = 1000. * del(i,k)
1001  if(ncloud > 0) then
1002  ptem = c0t(i,k) + c1
1003  qlk = dq / (eta(i,k) + etah * ptem * dz)
1004  dellal(i,k) = etah * c1 * dz * qlk * g / dp
1005  else
1006  qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1007  endif
1008  qcko(i,k) = qlk + qrch
1009  pwo(i,k) = etah * c0t(i,k) * dz * qlk
1010  cnvwt(i,k) = etah * qlk * g / dp
1011  endif
1012  endif
1013  endif
1014  enddo
1015  enddo
1016 !
1017 ! compute updraft velocity square(wu2)
1019 !
1020 ! bb1 = 2. * (1.+bet1*cd1)
1021 ! bb2 = 2. / (f1*(1.+gam1))
1022 !
1023 ! bb1 = 3.9
1024 ! bb2 = 0.67
1025 !
1026 ! bb1 = 2.0
1027 ! bb2 = 4.0
1028 !
1029  bb1 = 4.0
1030  bb2 = 0.8
1031 !
1032 ! do i = 1, im
1033 ! if (cnvflg(i)) then
1034 ! k = kbcon1(i)
1035 ! tem = po(i,k) / (rd * to(i,k))
1036 ! wucb = -0.01 * dot(i,k) / (tem * g)
1037 ! if(wucb > 0.) then
1038 ! wu2(i,k) = wucb * wucb
1039 ! else
1040 ! wu2(i,k) = 0.
1041 ! endif
1042 ! endif
1043 ! enddo
1044  do k = 2, km1
1045  do i = 1, im
1046  if (cnvflg(i)) then
1047  if(k > kbcon1(i) .and. k < ktcon(i)) then
1048  dz = zi(i,k) - zi(i,k-1)
1049  tem = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz
1050  tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz
1051  ptem = (1. - tem) * wu2(i,k-1)
1052  ptem1 = 1. + tem
1053  wu2(i,k) = (ptem + tem1) / ptem1
1054  wu2(i,k) = max(wu2(i,k), 0.)
1055  endif
1056  endif
1057  enddo
1058  enddo
1059 !
1060 ! compute updraft velocity averaged over the whole cumulus
1061 !
1063  do i = 1, im
1064  wc(i) = 0.
1065  sumx(i) = 0.
1066  enddo
1067  do k = 2, km1
1068  do i = 1, im
1069  if (cnvflg(i)) then
1070  if(k > kbcon1(i) .and. k < ktcon(i)) then
1071  dz = zi(i,k) - zi(i,k-1)
1072  tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1073  wc(i) = wc(i) + tem * dz
1074  sumx(i) = sumx(i) + dz
1075  endif
1076  endif
1077  enddo
1078  enddo
1079  do i = 1, im
1080  if(cnvflg(i)) then
1081  if(sumx(i) == 0.) then
1082  cnvflg(i)=.false.
1083  else
1084  wc(i) = wc(i) / sumx(i)
1085  endif
1086  val = 1.e-4
1087  if (wc(i) < val) cnvflg(i)=.false.
1088  endif
1089  enddo
1090 c
1091 c exchange ktcon with ktcon1
1092 c
1093  do i = 1, im
1094  if(cnvflg(i)) then
1095  kk = ktcon(i)
1096  ktcon(i) = ktcon1(i)
1097  ktcon1(i) = kk
1098  endif
1099  enddo
1100 c
1101 c this section is ready for cloud water
1102 c
1103  if(ncloud > 0) then
1104 c
1105 c compute liquid and vapor separation at cloud top
1106 c
1108  do i = 1, im
1109  if(cnvflg(i)) then
1110  k = ktcon(i) - 1
1111  gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1112  qrch = qeso(i,k)
1113  & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1114  dq = qcko(i,k) - qrch
1115 c
1116 c check if there is excess moisture to release latent heat
1117 c
1118  if(dq > 0.) then
1119  qlko_ktcon(i) = dq
1120  qcko(i,k) = qrch
1121  endif
1122  endif
1123  enddo
1124  endif
1125 c
1126 c--- compute precipitation efficiency in terms of windshear
1127 c
1128 !! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 :
1129 !! \f[
1130 !! 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
1131 !! \f]
1132 !! 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$.
1133  do i = 1, im
1134  if(cnvflg(i)) then
1135  vshear(i) = 0.
1136  endif
1137  enddo
1138  do k = 2, km
1139  do i = 1, im
1140  if (cnvflg(i)) then
1141  if(k > kb(i) .and. k <= ktcon(i)) then
1142  shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
1143  & + (vo(i,k)-vo(i,k-1)) ** 2)
1144  vshear(i) = vshear(i) + shear
1145  endif
1146  endif
1147  enddo
1148  enddo
1149  do i = 1, im
1150  if(cnvflg(i)) then
1151  vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
1152  e1=1.591-.639*vshear(i)
1153  & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1154  edt(i)=1.-e1
1155  val = .9
1156  edt(i) = min(edt(i),val)
1157  val = .0
1158  edt(i) = max(edt(i),val)
1159  endif
1160  enddo
1161 c
1162 c--- what would the change be, that a cloud with unit mass
1163 c--- will do to the environment?
1164 c
1167  do k = 1, km
1168  do i = 1, im
1169  if(cnvflg(i) .and. k <= kmax(i)) then
1170  dellah(i,k) = 0.
1171  dellaq(i,k) = 0.
1172  dellau(i,k) = 0.
1173  dellav(i,k) = 0.
1174  endif
1175  enddo
1176  enddo
1177 c
1178 c--- changed due to subsidence and entrainment
1179 c
1180  do k = 2, km1
1181  do i = 1, im
1182  if (cnvflg(i)) then
1183  if(k > kb(i) .and. k < ktcon(i)) then
1184  dp = 1000. * del(i,k)
1185  dz = zi(i,k) - zi(i,k-1)
1186 c
1187  dv1h = heo(i,k)
1188  dv2h = .5 * (heo(i,k) + heo(i,k-1))
1189  dv3h = heo(i,k-1)
1190  dv1q = qo(i,k)
1191  dv2q = .5 * (qo(i,k) + qo(i,k-1))
1192  dv3q = qo(i,k-1)
1193 c
1194  tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1195  tem1 = xlamud(i)
1196 cj
1197  dellah(i,k) = dellah(i,k) +
1198  & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h
1199  & - tem*eta(i,k-1)*dv2h*dz
1200  & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1201  & ) *g/dp
1202 cj
1203  dellaq(i,k) = dellaq(i,k) +
1204  & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q
1205  & - tem*eta(i,k-1)*dv2q*dz
1206  & + tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz
1207  & ) *g/dp
1208 cj
1209  tem1=eta(i,k)*(uo(i,k)-ucko(i,k))
1210  tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1))
1211  dellau(i,k) = dellau(i,k) + (tem1-tem2) * g/dp
1212 cj
1213  tem1=eta(i,k)*(vo(i,k)-vcko(i,k))
1214  tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1))
1215  dellav(i,k) = dellav(i,k) + (tem1-tem2) * g/dp
1216 cj
1217  endif
1218  endif
1219  enddo
1220  enddo
1221 c
1222 c------- cloud top
1223 c
1224  do i = 1, im
1225  if(cnvflg(i)) then
1226  indx = ktcon(i)
1227  dp = 1000. * del(i,indx)
1228  dv1h = heo(i,indx-1)
1229  dellah(i,indx) = eta(i,indx-1) *
1230  & (hcko(i,indx-1) - dv1h) * g / dp
1231  dv1q = qo(i,indx-1)
1232  dellaq(i,indx) = eta(i,indx-1) *
1233  & (qcko(i,indx-1) - dv1q) * g / dp
1234  dellau(i,indx) = eta(i,indx-1) *
1235  & (ucko(i,indx-1) - uo(i,indx-1)) * g / dp
1236  dellav(i,indx) = eta(i,indx-1) *
1237  & (vcko(i,indx-1) - vo(i,indx-1)) * g / dp
1238 c
1239 c cloud water
1240 c
1241  dellal(i,indx) = eta(i,indx-1) *
1242  & qlko_ktcon(i) * g / dp
1243  endif
1244  enddo
1245 !
1246 !
1247 ! compute convective turn-over time
1248 !
1250  do i= 1, im
1251  if(cnvflg(i)) then
1252  tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
1253  dtconv(i) = tem / wc(i)
1254  tfac = 1. + gdx(i) / 75000.
1255  dtconv(i) = tfac * dtconv(i)
1256  dtconv(i) = max(dtconv(i),dtmin)
1257  dtconv(i) = max(dtconv(i),dt2)
1258  dtconv(i) = min(dtconv(i),dtmax)
1259  endif
1260  enddo
1261 !
1263  do i= 1, im
1264  if(cnvflg(i)) then
1265  sumx(i) = 0.
1266  umean(i) = 0.
1267  endif
1268  enddo
1269  do k = 2, km1
1270  do i = 1, im
1271  if(cnvflg(i)) then
1272  if(k >= kbcon1(i) .and. k < ktcon1(i)) then
1273  dz = zi(i,k) - zi(i,k-1)
1274  tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
1275  umean(i) = umean(i) + tem * dz
1276  sumx(i) = sumx(i) + dz
1277  endif
1278  endif
1279  enddo
1280  enddo
1281  do i= 1, im
1282  if(cnvflg(i)) then
1283  umean(i) = umean(i) / sumx(i)
1284  umean(i) = max(umean(i), 1.)
1285  tauadv(i) = gdx(i) / umean(i)
1286  endif
1287  enddo
1288 c
1289 c compute cloud base mass flux as a function of the mean
1290 c updraft velcoity
1291 c
1294  do i= 1, im
1295  if(cnvflg(i)) then
1296  k = kbcon(i)
1297  rho = po(i,k)*100. / (rd*to(i,k))
1298  tfac = tauadv(i) / dtconv(i)
1299  tfac = min(tfac, 1.)
1300  xmb(i) = tfac*betaw*rho*wc(i)
1301  endif
1302  enddo
1303 !
1305  do i = 1, im
1306  if(cnvflg(i)) then
1307  tem = min(max(xlamue(i,kbcon(i)), 2.e-4), 6.e-4)
1308  tem = 0.2 / tem
1309  tem1 = 3.14 * tem * tem
1310  sigmagfm(i) = tem1 / garea(i)
1311  sigmagfm(i) = max(sigmagfm(i), 0.001)
1312  sigmagfm(i) = min(sigmagfm(i), 0.999)
1313  endif
1314  enddo
1315 !
1317  do i = 1, im
1318  if(cnvflg(i)) then
1319  if (gdx(i) < dxcrt) then
1320  scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
1321  scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
1322  else
1323  scaldfunc(i) = 1.0
1324  endif
1325  xmb(i) = xmb(i) * scaldfunc(i)
1326  xmb(i) = min(xmb(i),xmbmax(i))
1327  endif
1328  enddo
1331 c
1332 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1333 c
1334  do k = 1, km
1335  do i = 1, im
1336  if (cnvflg(i) .and. k <= kmax(i)) then
1337  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1338  qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
1339  val = 1.e-8
1340  qeso(i,k) = max(qeso(i,k), val )
1341  endif
1342  enddo
1343  enddo
1344 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1345 c
1349  do i = 1, im
1350  delhbar(i) = 0.
1351  delqbar(i) = 0.
1352  deltbar(i) = 0.
1353  delubar(i) = 0.
1354  delvbar(i) = 0.
1355  qcond(i) = 0.
1356  enddo
1357  do k = 1, km
1358  do i = 1, im
1359  if (cnvflg(i)) then
1360  if(k > kb(i) .and. k <= ktcon(i)) then
1361  dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1362  t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1363  q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1364 ! tem = 1./rcs(i)
1365 ! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1366 ! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1367  u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
1368  v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
1369  dp = 1000. * del(i,k)
1370  delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
1371  delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
1372  deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
1373  delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
1374  delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
1375  endif
1376  endif
1377  enddo
1378  enddo
1379 !
1381  do k = 1, km
1382  do i = 1, im
1383  if (cnvflg(i)) then
1384  if(k > kb(i) .and. k <= ktcon(i)) then
1385  qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1386  qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
1387  val = 1.e-8
1388  qeso(i,k) = max(qeso(i,k), val )
1389  endif
1390  endif
1391  enddo
1392  enddo
1393 c
1395  do i = 1, im
1396  rntot(i) = 0.
1397  delqev(i) = 0.
1398  delq2(i) = 0.
1399  flg(i) = cnvflg(i)
1400  enddo
1401  do k = km, 1, -1
1402  do i = 1, im
1403  if (cnvflg(i)) then
1404  if(k < ktcon(i) .and. k > kb(i)) then
1405  rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
1406  endif
1407  endif
1408  enddo
1409  enddo
1410 c
1411 c evaporating rain
1412 c
1416  do k = km, 1, -1
1417  do i = 1, im
1418  if (k <= kmax(i)) then
1419  deltv(i) = 0.
1420  delq(i) = 0.
1421  qevap(i) = 0.
1422  if(cnvflg(i)) then
1423  if(k < ktcon(i) .and. k > kb(i)) then
1424  rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
1425  endif
1426  endif
1427  if(flg(i) .and. k < ktcon(i)) then
1428  evef = edt(i) * evfact
1429  if(islimsk(i) == 1) evef=edt(i) * evfactl
1430 ! if(islimsk(i) == 1) evef=.07
1431 c if(islimsk(i) == 1) evef = 0.
1432  qcond(i) = evef * (q1(i,k) - qeso(i,k))
1433  & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
1434  dp = 1000. * del(i,k)
1435  if(rn(i) > 0. .and. qcond(i) < 0.) then
1436  qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
1437  qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
1438  delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
1439  endif
1440  if(rn(i) > 0. .and. qcond(i) < 0. .and.
1441  & delq2(i) > rntot(i)) then
1442  qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
1443  flg(i) = .false.
1444  endif
1445  if(rn(i) > 0. .and. qevap(i) > 0.) then
1446  tem = .001 * dp / g
1447  tem1 = qevap(i) * tem
1448  if(tem1 > rn(i)) then
1449  qevap(i) = rn(i) / tem
1450  rn(i) = 0.
1451  else
1452  rn(i) = rn(i) - tem1
1453  endif
1454  q1(i,k) = q1(i,k) + qevap(i)
1455  t1(i,k) = t1(i,k) - elocp * qevap(i)
1456  deltv(i) = - elocp*qevap(i)/dt2
1457  delq(i) = + qevap(i)/dt2
1458  delqev(i) = delqev(i) + .001*dp*qevap(i)/g
1459  endif
1460  delqbar(i) = delqbar(i) + delq(i)*dp/g
1461  deltbar(i) = deltbar(i) + deltv(i)*dp/g
1462  endif
1463  endif
1464  enddo
1465  enddo
1466 cj
1467 ! do i = 1, im
1468 ! if(me == 31 .and. cnvflg(i)) then
1469 ! if(cnvflg(i)) then
1470 ! print *, ' shallow delhbar, delqbar, deltbar = ',
1471 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
1472 ! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
1473 ! print *, ' precip =', hvap*rn(i)*1000./dt2
1474 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
1475 ! endif
1476 ! enddo
1477 cj
1478  do i = 1, im
1479  if(cnvflg(i)) then
1480  if(rn(i) < 0. .or. .not.flg(i)) rn(i) = 0.
1481  ktop(i) = ktcon(i)
1482  kbot(i) = kbcon(i)
1483  kcnv(i) = 2
1484  endif
1485  enddo
1486 c
1487 c convective cloud water
1488 c
1490  do k = 1, km
1491  do i = 1, im
1492  if (cnvflg(i)) then
1493  if (k >= kbcon(i) .and. k < ktcon(i)) then
1494  cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
1495  endif
1496  endif
1497  enddo
1498  enddo
1499 
1500 c
1501 c convective cloud cover
1502 c
1504  do k = 1, km
1505  do i = 1, im
1506  if (cnvflg(i)) then
1507  if (k >= kbcon(i) .and. k < ktcon(i)) then
1508  cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
1509  cnvc(i,k) = min(cnvc(i,k), 0.2)
1510  cnvc(i,k) = max(cnvc(i,k), 0.0)
1511  endif
1512  endif
1513  enddo
1514  enddo
1515 c
1516 c cloud water
1517 c
1519  if (ncloud > 0) then
1520 !
1521  do k = 1, km1
1522  do i = 1, im
1523  if (cnvflg(i)) then
1524 ! if (k > kb(i) .and. k <= ktcon(i)) then
1525  if (k >= kbcon(i) .and. k <= ktcon(i)) then
1526  tem = dellal(i,k) * xmb(i) * dt2
1527  tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
1528  if (ql(i,k,2) > -999.0) then
1529  ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
1530  ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
1531  else
1532  ql(i,k,1) = ql(i,k,1) + tem
1533  endif
1534  endif
1535  endif
1536  enddo
1537  enddo
1538 !
1539  endif
1540 !
1541 ! hchuang code change
1542 !
1544 !
1546  do k = 1, km
1547  do i = 1, im
1548  if(cnvflg(i)) then
1549  if(k >= kb(i) .and. k < ktop(i)) then
1550  ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
1551  endif
1552  endif
1553  enddo
1554  enddo
1556  do i = 1, im
1557  if(cnvflg(i)) then
1558  k = ktop(i)-1
1559  dt_mf(i,k) = ud_mf(i,k)
1560  endif
1561  enddo
1562 !!
1563  return
1564  end
samfshalcnv
subroutine samfshalcnv(im, ix, km, delt, delp, prslp, psp, phil, ql,
This subroutine contains the entirety of the SAMF shallow convection scheme.
Definition: samfshalcnv.f:62