CCPP SciDoc v7.0.x  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
samfshalcnv.f
1
3
5 module samfshalcnv
6
7 use samfcnv_aerosols, only : samfshalcnv_aerosols
8 use progsigma, only : progsigma_calc
9 use progomega, only : progomega_calc
10 contains
11
12 subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, &
13 & errmsg, errflg)
14
15 integer, intent(in) :: imfshalcnv
16 integer, intent(in) :: imfshalcnv_samf
17
18 ! CCPP error handling
19 character(len=*), intent(out) :: errmsg
20 integer, intent(out) :: errflg
21
22 ! Consistency checks
23 if (imfshalcnv/=imfshalcnv_samf) then
24 write(errmsg,'(*(a))') 'Logic error: namelist choice of', &
25 & ' shallow convection is different from SAMF'
26 errflg = 1
27 return
28 end if
29 end subroutine samfshalcnv_init
30
53 subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
54 & eps,epsm1,fv,grav,hvap,rd,rv, &
55 & t0c,delt,ntk,ntr,delp,first_time_step,restart, &
56 & tmf,qmicro,progsigma,progomega, &
57 & prslp,psp,phil,tkeh,qtr,prevsq,q,q1,t1,u1,v1,fscav, &
58 & rn,kbot,ktop,kcnv,islimsk,garea,cscale, &
59 & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, &
60 & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, &
61 & sigmain,sigmaout,omegain,omegaout,betadcu,betamcu,betascu, &
62 & errmsg,errflg)
63!
64 use machine , only : kind_phys
65 use funcphys , only : fpvs
66
67 implicit none
68!
69 integer, intent(in) :: im, km, itc, ntc, ntk, ntr, ncloud
70 integer, intent(in) :: islimsk(:)
71 real(kind=kind_phys), intent(in) :: cliq, cp, cvap, &
72 & eps, epsm1, fv, grav, hvap, rd, rv, t0c, betascu, betadcu, &
73 & betamcu
74 real(kind=kind_phys), intent(in) :: delt, cscale
75 real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), &
76 & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), &
77 & tmf(:,:,:), q(:,:)
78 real(kind=kind_phys), intent(in), optional :: qmicro(:,:), &
79 & prevsq(:,:)
80 real(kind=kind_phys), intent(in), optional :: sigmain(:,:), &
81 & omegain(:,:)
82!
83 real(kind=kind_phys), dimension(:), intent(in) :: fscav
84 integer, intent(inout) :: kcnv(:)
85 ! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH
86 real(kind=kind_phys), intent(inout) :: qtr(:,:,:), &
87 & q1(:,:), t1(:,:), u1(:,:), v1(:,:), tkeh(:,:)
88!
89 integer, intent(out) :: kbot(:), ktop(:)
90 real(kind=kind_phys), intent(out) :: rn(:), &
91 & cnvw(:,:), cnvc(:,:), dt_mf(:,:)
92!
93 real(kind=kind_phys), intent(out) :: ud_mf(:,:)
94 real(kind=kind_phys), intent(inout), optional :: sigmaout(:,:), &
95 & omegaout(:,:)
96
97 real(kind=kind_phys), intent(in) :: clam, c0s, c1, &
98 & asolfac, evef, pgcon
99 logical, intent(in) :: hwrf_samfshal,first_time_step, &
100 & restart,progsigma,progomega
101 character(len=*), intent(out) :: errmsg
102 integer, intent(out) :: errflg
103!
104! local variables
105 integer i,j,indx, k, kk, km1, n
106 integer kpbl(im)
107!
108 real(kind=kind_phys) clamd, tkemx, tkemn, dtke
109!
110 real(kind=kind_phys) dellat,
111 & c0l, d0,
112 & desdt, dp,
113 & dq, dqsdp, dqsdt, dt,
114 & dt2, dtmax, dtmin,
115 & dxcrt, dxcrtc0,
116 & dv1h, dv2h, dv3h,
117 & dz, dz1, e1,
118 & el2orc, elocp, aafac,
119 & cm, cq,
120 & es, etah, h1, shevf,
121! & evfact, evfactl,
122 & fact1, fact2, factor,
123 & cthk, cthkmn, dthk,
124 & gamma, pprime, betaw, tauadv,
125 & qlk, qrch, qs,
126 & rfact, shear, tfac,
127 & val, val1, val2,
128 & w1, w1l, w1s, w2,
129 & w2l, w2s, w3, w3l,
130 & w3s, w4, w4l, w4s,
131 & rho, tem, tem1, tem2,
132 & ptem, ptem1
133!
134 integer kb(im), kb1(im), kbcon(im), kbcon1(im),
135 & ktcon(im), ktcon1(im),
136 & kbm(im), kmax(im)
137!
138 real(kind=kind_phys) aa1(im), cina(im),
139 & tkemean(im), clamt(im),
140 & ps(im), del(im,km), prsl(im,km),
141 & umean(im), advfac(im), gdx(im),
142 & delhbar(im), delq(im), delq2(im),
143 & delqbar(im), delqev(im), deltbar(im),
144! & deltv(im), dtconv(im), edt(im),
145 & deltv(im), dtconv(im),
146 & pdot(im), po(im,km),
147 & qcond(im), qevap(im), hmax(im),
148! & rntot(im), vshear(im),
149 & rntot(im),
150 & xlamud(im), xmb(im), xmbmax(im),
151 & delebar(im,ntr),
152 & delubar(im), delvbar(im)
153!
154 real(kind=kind_phys) c0(im), sfcpbl(im)
155c
156 real(kind=kind_phys) crtlame, crtlamd
157!
158 real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn,
159 & cinacr, cinacrmx, cinacrmn,
160 & sfclfac, rhcrt,
161 & tkcrt, cmxfac
162!
163! parameters for updraft velocity calculation
164 real(kind=kind_phys) bb1, bb2, csmf, wucb
165cc
166
167! parameters for prognostic sigma closure
168 real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
169 & omegac(im),zeta(im,km),dbyo1(im,km),
170 & sigmab(im),qadv(im,km)
171 real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins,
172 & sigminm
173 logical flag_shallow,flag_mid
174c physical parameters
175! parameter(g=grav,asolfac=0.89)
176! parameter(g=grav)
177! parameter(elocp=hvap/cp,
178! & el2orc=hvap*hvap/(rv*cp))
179! parameter(c0s=0.002,c1=5.e-4,d0=.01)
180! parameter(d0=.01)
181 parameter(d0=.001)
182! parameter(c0l=c0s*asolfac)
183!
184! asolfac: aerosol-aware parameter based on Lim & Hong (2012)
185! asolfac= cx / c0s(=.002)
186! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s)
187! Nccn: CCN number concentration in cm^(-3)
188! Until a realistic Nccn is provided, Nccns are assumed
189! as Nccn=100 for sea and Nccn=1000 for land
190!
191 parameter(cm=1.0,cq=1.0)
192! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
193 parameter(clamd=0.1,tkemx=0.65,tkemn=0.05)
194 parameter(dtke=tkemx-tkemn)
195 parameter(cthk=200.,cthkmn=0.,dthk=25.)
196 parameter(sfclfac=0.2,rhcrt=0.75)
197 parameter(cinpcrmx=180.,cinpcrmn=120.)
198! shevf is an enhancing evaporation factor for shallow convection
199 parameter(cinacrmx=-120.,shevf=2.0)
200 parameter(dtmax=10800.,dtmin=600.)
201 parameter(bb1=4.0,bb2=0.8,csmf=0.2)
202 parameter(tkcrt=2.,cmxfac=10.)
203! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
204 parameter(betaw=.03,dxcrtc0=9.e3)
205 parameter(h1=0.33333333)
206! progsigma
207 parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01)
208c local variables and arrays
209 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
210 & uo(im,km), vo(im,km), qeso(im,km),
211 & ctr(im,km,ntr), ctro(im,km,ntr)
212! for aerosol transport
213! real(kind=kind_phys) qaero(im,km,ntc)
214c variables for tracer wet deposition,
215 real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw,
216 & wet_dep
217 real(kind=kind_phys), parameter :: escav = 0.8 ! wet scavenging efficiency
218!
219! for updraft velocity calculation
220 real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km),
221 & wush(im,km), wc(im)
222!
223! for updraft fraction & scale-aware function
224 real(kind=kind_phys) scaldfunc(im), sigmagfm(im)
225!
226c cloud water
227! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
228 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km),
229 & dbyo(im,km), zo(im,km), xlamue(im,km),
230 & rh(im,km),
231 & heo(im,km), heso(im,km),
232 & dellah(im,km), dellaq(im,km),
233 & dellae(im,km,ntr),
234 & dellau(im,km), dellav(im,km), hcko(im,km),
235 & ucko(im,km), vcko(im,km), qcko(im,km),
236 & qrcko(im,km), ecko(im,km,ntr),
237 & ercko(im,km,ntr), eta(im,km),
238 & zi(im,km), pwo(im,km), c0t(im,km),
239 & sumx(im), tx1(im), cnvwt(im,km)
240 &, rhbar(im)
241!
242! variables for Total Variation Diminishing (TVD) flux-limiter scheme
243! on environmental subsidence and uplifting
244!
245 real(kind=kind_phys) q_diff(im,0:km-1), e_diff(im,0:km-1,ntr),
246 & flxtvd(im,km-1)
247 real(kind=kind_phys) rrkp, phkp
248 real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im)
249!
250 logical do_aerosols, totflg, cnvflg(im), flg(im)
251!
252 real(kind=kind_phys) tf, tcr, tcrf
253 parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
254
255
256c-----------------------------------------------------------------------
257!
258! Initialize CCPP error handling variables
259 errmsg = ''
260 errflg = 0
261
262 gravinv = 1./grav
263 invdelt = 1./delt
264
265 elocp = hvap/cp
266 el2orc = hvap*hvap/(rv*cp)
267
268 fact1 = (cvap-cliq)/rv
269 fact2 = hvap/rv-fact1*t0c
270
271 if (.not.hwrf_samfshal) then
272 cinacrmn=-80.
273 endif
274
275 if (progsigma) then
276 dxcrt=10.e3
277 else
278 dxcrt=15.e3
279 endif
280
281c-----------------------------------------------------------------------
282 if (.not.hwrf_samfshal) then
284 do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0)
285 if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3)
286 endif
287!
288!************************************************************************
289! convert input Pa terms to Cb terms -- Moorthi
292 ps = psp * 0.001
293 prsl = prslp * 0.001
294 del = delp * 0.001
295!************************************************************************
296!
297 km1 = km - 1
298c
299c initialize arrays
300c
302!
303 chem_c = 0.
304 chem_pw = 0.
305 wet_dep = 0.
306!
307 if(hwrf_samfshal) then
308 do i=1,im
309 cnvflg(i) = .true.
310 if(kcnv(i) == 1) cnvflg(i) = .false.
311 if(cnvflg(i)) then
312 kbot(i)=km+1
313 ktop(i)=0
314 endif
315 sfcpbl(i) = sfclfac * hpbl(i)
316 rn(i)=0.
317 kbcon(i)=km
318 ktcon(i)=1
319 kb(i)=km
320 pdot(i) = 0.
321 qlko_ktcon(i) = 0.
322! edt(i) = 0.
323 aa1(i) = 0.
324 cina(i) = 0.
325! vshear(i) = 0.
326 advfac(i) = 0.
327 gdx(i) = sqrt(garea(i))
328 xmb(i) = 0.
329 scaldfunc(i)=-1.0 ! wang initialized
330 sigmagfm(i)=-1.0
331 enddo
332
333 else !gfs_samfshal
334 do i=1,im
335 cnvflg(i) = .true.
336 if(kcnv(i) == 1) cnvflg(i) = .false.
337 if(cnvflg(i)) then
338 kbot(i)=km+1
339 ktop(i)=0
340 endif
341 sfcpbl(i) = sfclfac * hpbl(i)
342 rn(i)=0.
343 kbcon(i)=km
344 ktcon(i)=1
345 kb(i)=km
346 pdot(i) = 0.
347 qlko_ktcon(i) = 0.
348! edt(i) = 0.0
349 aa1(i) = 0.
350 cina(i) = 0.
351! vshear(i) = 0.
352 gdx(i) = sqrt(garea(i))
353 xmb(i) = 0.
354 enddo
355 endif
356!!
358 totflg = .true.
359 do i=1,im
360 totflg = totflg .and. (.not. cnvflg(i))
361 enddo
362 if(totflg) return
363!!
365 do i=1,im
366 if(islimsk(i) == 1) then
367 c0(i) = c0s*asolfac
368 else
369 c0(i) = c0s
370 endif
371 enddo
372!
374 do i=1,im
375 if(gdx(i) < dxcrtc0) then
376 tem = gdx(i) / dxcrtc0
377 tem1 = tem**3
378 c0(i) = c0(i) * tem1
379 endif
380 enddo
381!
383 do k = 1, km
384 do i = 1, im
385 if(t1(i,k) > 273.16) then
386 c0t(i,k) = c0(i)
387 else
388 tem = d0 * (t1(i,k) - 273.16)
389 tem1 = exp(tem)
390 c0t(i,k) = c0(i) * tem1
391 endif
392 enddo
393 enddo
394!
396 do k = 1, km
397 do i = 1, im
398 cnvw(i,k) = 0.
399 cnvc(i,k) = 0.
400 enddo
401 enddo
402! hchuang code change
404 do k = 1, km
405 do i = 1, im
406 ud_mf(i,k) = 0.
407 dt_mf(i,k) = 0.
408 enddo
409 enddo
410c
411 dt2 = delt
412!
413c model tunable parameters are all here
414 aafac = .1
415! evfact = 0.3
416! evfactl = 0.3
417!
418 crtlame = 1.0e-4
419 crtlamd = 3.0e-4
420!
421 w1l = -8.e-3
422 w2l = -4.e-2
423 w3l = -5.e-3
424 w4l = -5.e-4
425 w1s = -2.e-4
426 w2s = -2.e-3
427 w3s = -1.e-3
428 w4s = -2.e-5
429c
430c define top layer for search of the downdraft originating layer
431c and the maximum thetae for updraft
432c
434 do i=1,im
435 kbm(i) = km
436 kmax(i) = km
437 tx1(i) = 1.0 / ps(i)
438 enddo
439!
440 do k = 1, km
441 do i=1,im
442 if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1
443 if (prsl(i,k)*tx1(i) > 0.60) kmax(i) = k + 1
444 enddo
445 enddo
446 do i=1,im
447 kbm(i) = min(kbm(i),kmax(i))
448 enddo
449c
450c hydrostatic height assume zero terr and compute
451c updraft entrainment rate as an inverse function of height
452c
454 do k = 1, km
455 do i=1,im
456 zo(i,k) = phil(i,k) / grav
457 enddo
458 enddo
460 if(hwrf_samfshal) then
461 do k = 1, km1
462 do i=1,im
463 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
464 xlamue(i,k) = clam / zi(i,k)
465 enddo
466 enddo
467 do i=1,im
468 xlamue(i,km) = xlamue(i,km1)
469 enddo
470 else
471 do k = 1, km1
472 do i=1,im
473 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
474 enddo
475 enddo
476 endif
477c
478c pbl height
479c
481 do i=1,im
482 flg(i) = cnvflg(i)
483 kpbl(i)= 1
484 enddo
485 do k = 2, km1
486 do i=1,im
487 if (flg(i) .and. zo(i,k) <= hpbl(i)) then
488 kpbl(i) = k
489 else
490 flg(i) = .false.
491 endif
492 enddo
493 enddo
494 do i=1,im
495 kpbl(i)= min(kpbl(i),kbm(i))
496 enddo
497c
498c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
499c convert surface pressure to mb from cb
500c
502 do k = 1, km
503 do i = 1, im
504 if (cnvflg(i) .and. k <= kmax(i)) then
505 pfld(i,k) = prsl(i,k) * 10.0
506 eta(i,k) = 1.
507 rh(i,k) = 0.
508 hcko(i,k) = 0.
509 qcko(i,k) = 0.
510 qrcko(i,k)= 0.
511 ucko(i,k) = 0.
512 vcko(i,k) = 0.
513 dbyo(i,k) = 0.
514 pwo(i,k) = 0.
515 dellal(i,k) = 0.
516 to(i,k) = t1(i,k)
517 qo(i,k) = q1(i,k)
518 uo(i,k) = u1(i,k)
519 vo(i,k) = v1(i,k)
520! uo(i,k) = u1(i,k) * rcs(i)
521! vo(i,k) = v1(i,k) * rcs(i)
522 wu2(i,k) = 0.
523 buo(i,k) = 0.
524 wush(i,k) = 0.
525 drag(i,k) = 0.
526 cnvwt(i,k) = 0.
527 endif
528 enddo
529 enddo
530
531 do i = 1,im
532 omegac(i)=0.
533 enddo
534
535 do k = 1, km
536 do i = 1, im
537 dbyo1(i,k)=0.
538 zdqca(i,k)=0.
539 omega_u(i,k)=0.
540 zeta(i,k)=1.0
541 enddo
542 enddo
543!
544! initialize tracer variables
545!
546 if (.not.hwrf_samfshal) then
547 do n = 3, ntr+2
548 kk = n-2
549 do k = 1, km
550 do i = 1, im
551 if (cnvflg(i) .and. k <= kmax(i)) then
552 ctr(i,k,kk) = qtr(i,k,n)
553 ctro(i,k,kk) = qtr(i,k,n)
554 ecko(i,k,kk) = 0.
555 ercko(i,k,kk) = 0.
556 endif
557 enddo
558 enddo
559 enddo
560 endif
562 do k = 1, km
563 do i=1,im
564 if (cnvflg(i) .and. k <= kmax(i)) then
565 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
566 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
567 val1 = 1.e-8
568 qeso(i,k) = max(qeso(i,k), val1)
569 val2 = 1.e-10
570 qo(i,k) = max(qo(i,k), val2 )
571! qo(i,k) = min(qo(i,k),qeso(i,k))
572! tvo(i,k) = to(i,k) + fv * to(i,k) * qo(i,k)
573 endif
574 enddo
575 enddo
576c
577c compute moist static energy
578c
580 do k = 1, km
581 do i=1,im
582 if (cnvflg(i) .and. k <= kmax(i)) then
583! tem = grav * zo(i,k) + cp * to(i,k)
584 tem = phil(i,k) + cp * to(i,k)
585 heo(i,k) = tem + hvap * qo(i,k)
586 heso(i,k) = tem + hvap * qeso(i,k)
587c heo(i,k) = min(heo(i,k),heso(i,k))
588 endif
589 enddo
590 enddo
591c
592c determine level with largest moist static energy within pbl
593c this is the level where updraft starts
594c
597 do i=1,im
598 flg(i) = cnvflg(i)
599 kb1(i) = 1
600 enddo
601 do k = 1, km1
602 do i=1,im
603 if (flg(i) .and. zo(i,k) <= sfcpbl(i)) then
604 kb1(i) = k
605 else
606 flg(i) = .false.
607 endif
608 enddo
609 enddo
610 do i=1,im
611 kb1(i) = min(kb1(i),kpbl(i))
612 enddo
613c
615 do i=1,im
616 if (cnvflg(i)) then
617 hmax(i) = heo(i,kb1(i))
618 kb(i) = kb1(i)
619 endif
620 enddo
621 do k = 2, km
622 do i=1,im
623 if(cnvflg(i) .and. (k > kb1(i) .and. k <= kpbl(i))) then
624 if(heo(i,k) > hmax(i)) then
625 kb(i) = k
626 hmax(i) = heo(i,k)
627 endif
628 endif
629 enddo
630 enddo
631c
633 do k = 1, km1
634 do i=1,im
635 if (cnvflg(i) .and. k <= kmax(i)-1) then
636 dz = .5 * (zo(i,k+1) - zo(i,k))
637 dp = .5 * (pfld(i,k+1) - pfld(i,k))
638 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
639 pprime = pfld(i,k+1) + epsm1 * es
640 qs = eps * es / pprime
641 dqsdp = - qs / pprime
642 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
643 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
644 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
645 dt = (grav*dz + hvap*dqsdp*dp) / (cp*(1. + gamma))
646 dq = dqsdt * dt + dqsdp * dp
647 to(i,k) = to(i,k+1) + dt
648 qo(i,k) = qo(i,k+1) + dq
649 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
650 endif
651 enddo
652 enddo
653!
655 do k = 1, km1
656 do i=1,im
657 if (cnvflg(i) .and. k <= kmax(i)-1) then
658 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
659 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
660 val1 = 1.e-8
661 qeso(i,k) = max(qeso(i,k), val1)
662 val2 = 1.e-10
663 qo(i,k) = max(qo(i,k), val2 )
664! qo(i,k) = min(qo(i,k),qeso(i,k))
665 rh(i,k) = min(qo(i,k)/qeso(i,k), 1.)
666 heo(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
667 & cp * to(i,k) + hvap * qo(i,k)
668 heso(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
669 & cp * to(i,k) + hvap * qeso(i,k)
670 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
671 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
672 endif
673 enddo
674 enddo
675
676 if (.not.hwrf_samfshal) then
677 do n = 1, ntr
678 do k = 1, km1
679 do i=1,im
680 if (cnvflg(i) .and. k <= kmax(i)-1) then
681 ctro(i,k,n) = .5 * (ctro(i,k,n) + ctro(i,k+1,n))
682 endif
683 enddo
684 enddo
685 enddo
686 endif
687c
688c look for the level of free convection as cloud base
689c
691 do i=1,im
692 flg(i) = cnvflg(i)
693 if(flg(i)) kbcon(i) = kmax(i)
694 enddo
695 do k = 2, km1
696 do i=1,im
697 if (flg(i) .and. k < kbm(i)) then
698 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
699 kbcon(i) = k
700 flg(i) = .false.
701 endif
702 endif
703 enddo
704 enddo
705c
706 do i=1,im
707 if(cnvflg(i)) then
708 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
709 endif
710 enddo
711!!
713 totflg = .true.
714 do i=1,im
715 totflg = totflg .and. (.not. cnvflg(i))
716 enddo
717 if(totflg) return
718!!
720 do i=1,im
721 if(cnvflg(i)) then
722! pdot(i) = 10.* dot(i,kbcon(i))
723 pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
724 endif
725 enddo
726c
727c turn off convection if pressure depth between parcel source level
728c and cloud base is larger than a critical value, cinpcr
729c
730 do i=1,im
731 if(cnvflg(i)) then
732 if(islimsk(i) == 1) then
733 w1 = w1l
734 w2 = w2l
735 w3 = w3l
736 w4 = w4l
737 else
738 w1 = w1s
739 w2 = w2s
740 w3 = w3s
741 w4 = w4s
742 endif
743 if(pdot(i) <= w4) then
744 tem = (pdot(i) - w4) / (w3 - w4)
745 elseif(pdot(i) >= -w4) then
746 tem = - (pdot(i) + w4) / (w4 - w3)
747 else
748 tem = 0.
749 endif
750 val1 = -1.
751 tem = max(tem,val1)
752 val2 = 1.
753 tem = min(tem,val2)
754 ptem = 1. - tem
755 ptem1= .5*(cinpcrmx-cinpcrmn)
756 cinpcr = cinpcrmx - ptem * ptem1
757 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
758
759 if(tem1 > cinpcr .and.
760 & zi(i,kbcon(i)) > hpbl(i)) then
761 cnvflg(i) = .false.
762 endif
763 endif
764 enddo
765!!
766 totflg = .true.
767 do i=1,im
768 totflg = totflg .and. (.not. cnvflg(i))
769 enddo
770 if(totflg) return
771!
772! re-define kb & kbcon
773!
774 do i=1,im
775 if (cnvflg(i)) then
776 hmax(i) = heo(i,1)
777 kb(i) = 1
778 endif
779 enddo
780 do k = 2, km
781 do i=1,im
782 if (cnvflg(i) .and. k <= kpbl(i)) then
783 if(heo(i,k) > hmax(i)) then
784 kb(i) = k
785 hmax(i) = heo(i,k)
786 endif
787 endif
788 enddo
789 enddo
790!
791 do i=1,im
792 flg(i) = cnvflg(i)
793 if(flg(i)) kbcon(i) = kmax(i)
794 enddo
795 do k = 2, km1
796 do i=1,im
797 if (flg(i) .and. k < kbm(i)) then
798 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
799 kbcon(i) = k
800 flg(i) = .false.
801 endif
802 endif
803 enddo
804 enddo
805!
806 do i=1,im
807 if(cnvflg(i)) then
808 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
809 endif
810 enddo
811!!
812 totflg = .true.
813 do i=1,im
814 totflg = totflg .and. (.not. cnvflg(i))
815 enddo
816 if(totflg) return
817!!
818 do i=1,im
819 if(cnvflg(i)) then
820! pdot(i) = 10.* dot(i,kbcon(i))
821 pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
822 endif
823 enddo
824!
826!
827 do i = 1, im
828 rhbar(i) = 0.
829 sumx(i) = 0.
830 enddo
831 do k = 1, km1
832 do i = 1, im
833 if (cnvflg(i)) then
834 if(k >= kb(i) .and. k < kbcon(i)) then
835 dz = zo(i,k+1) - zo(i,k)
836 rhbar(i) = rhbar(i) + rh(i,k) * dz
837 sumx(i) = sumx(i) + dz
838 endif
839 endif
840 enddo
841 enddo
842 do i= 1, im
843 if(cnvflg(i)) then
844 rhbar(i) = rhbar(i) / sumx(i)
845 if(rhbar(i) < rhcrt) then
846 cnvflg(i) = .false.
847 endif
848 endif
849 enddo
850!!
851 totflg = .true.
852 do i=1,im
853 totflg = totflg .and. (.not. cnvflg(i))
854 enddo
855 if(totflg) return
856!!
857!
858! turbulent entrainment rate assumed to be proportional
859! to subcloud mean TKE
860!
861!c
862!c specify the detrainment rate for the updrafts
863!c
864 if (hwrf_samfshal) then
865 do i = 1, im
866 if(cnvflg(i)) then
867 xlamud(i) = xlamue(i,kbcon(i))
868! xlamud(i) = crtlamd
869 endif
870 enddo
871 else
872 if(ntk > 0) then
873 do i= 1, im
874 if(cnvflg(i)) then
875 sumx(i) = 0.
876 tkemean(i) = 0.
877 endif
878 enddo
879!
880 do k = 1, km1
881 do i = 1, im
882 if(cnvflg(i)) then
883 if(k >= kb(i) .and. k < kbcon(i)) then
884 dz = zo(i,k+1) - zo(i,k)
885 tkemean(i) = tkemean(i) + tkeh(i,k) * dz
886 sumx(i) = sumx(i) + dz
887 endif
888 endif
889 enddo
890 enddo
891!
892 do i= 1, im
893 if(cnvflg(i)) then
894 tkemean(i) = tkemean(i) / sumx(i)
895 if(tkemean(i) > tkemx) then
896 clamt(i) = clam + clamd
897 else if(tkemean(i) < tkemn) then
898 clamt(i) = clam - clamd
899 else
900 tem = tkemx - tkemean(i)
901 tem1 = 1. - 2. * tem / dtke
902 clamt(i) = clam + clamd * tem1
903 endif
904 endif
905 enddo
906!
907 do i=1,im
908 if(cnvflg(i)) then
909 if(tkemean(i) > tkcrt) then
910 tem = 1. + tkemean(i)/tkcrt
911 tem1 = min(tem, cmxfac)
912 clamt(i) = tem1 * clam
913 endif
914 endif
915 enddo
916!
917 else
918!
919 do i= 1, im
920 if(cnvflg(i)) then
921 clamt(i) = clam
922 endif
923 enddo
924!
925 endif
926!!
927!
928! assume updraft entrainment rate
929! is an inverse function of height
930!
931 do k = 1, km1
932 do i=1,im
933 if(cnvflg(i)) then
934 dz = zo(i,k+1) - zo(i,k)
935 xlamue(i,k) = clamt(i) / (zi(i,k) + dz)
936 xlamue(i,k) = max(xlamue(i,k), crtlame)
937 endif
938 enddo
939 enddo
940 do i=1,im
941 if(cnvflg(i)) then
942 xlamue(i,km) = xlamue(i,km1)
943 endif
944 enddo
945c
946c specify the detrainment rate for the updrafts
947c
948!! (The updraft detrainment rate is set constant and equal to the entrainment rate at cloud base.)
949!!
951 do i = 1, im
952 if(cnvflg(i)) then
953! xlamud(i) = xlamue(i,kbcon(i))
954! xlamud(i) = crtlamd
955 xlamud(i) = 0.001 * clamt(i)
956 endif
957 enddo
958 endif ! hwrf_samfshal
959c
960c determine updraft mass flux for the subcloud layers
961c
967 do k = km1, 1, -1
968 do i = 1, im
969 if (cnvflg(i)) then
970 if(k < kbcon(i) .and. k >= kb(i)) then
971 dz = zi(i,k+1) - zi(i,k)
972 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
973 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
974 endif
975 endif
976 enddo
977 enddo
978c
979c compute mass flux above cloud base
980c
981 do i = 1, im
982 flg(i) = cnvflg(i)
983 enddo
984 do k = 2, km1
985 do i = 1, im
986 if(flg(i))then
987 if(k > kbcon(i) .and. k < kmax(i)) then
988 dz = zi(i,k) - zi(i,k-1)
989 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
990 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
991 if(eta(i,k) <= 0.) then
992 kmax(i) = k
993 kbm(i) = min(kbm(i),kmax(i))
994 flg(i) = .false.
995 endif
996 endif
997 endif
998 enddo
999 enddo
1000c
1001c compute updraft cloud property
1002c
1004 do i = 1, im
1005 if(cnvflg(i)) then
1006 indx = kb(i)
1007 hcko(i,indx) = heo(i,indx)
1008 ucko(i,indx) = uo(i,indx)
1009 vcko(i,indx) = vo(i,indx)
1010 endif
1011 enddo
1012! for tracers
1013 if (.not. hwrf_samfshal) then
1014 do n = 1, ntr
1015 do i = 1, im
1016 if(cnvflg(i)) then
1017 indx = kb(i)
1018 ecko(i,indx,n) = ctro(i,indx,n)
1019 ercko(i,indx,n) = ctro(i,indx,n)
1020 endif
1021 enddo
1022 enddo
1023 endif
1024c
1025! cm is an enhancement factor in entrainment rates for momentum
1026!
1028 do k = 2, km1
1029 do i = 1, im
1030 if (cnvflg(i)) then
1031 if(k > kb(i) .and. k < kmax(i)) then
1032 dz = zi(i,k) - zi(i,k-1)
1033 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1034 tem1 = 0.5 * xlamud(i) * dz
1035 factor = 1. + tem - tem1
1036 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
1037 & (heo(i,k)+heo(i,k-1)))/factor
1038 dbyo(i,k) = hcko(i,k) - heso(i,k)
1039!
1040 tem = 0.5 * cm * tem
1041 factor = 1. + tem
1042 ptem = tem + pgcon
1043 ptem1= tem - pgcon
1044 ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*uo(i,k)
1045 & +ptem1*uo(i,k-1))/factor
1046 vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*vo(i,k)
1047 & +ptem1*vo(i,k-1))/factor
1048 endif
1049 endif
1050 enddo
1051 enddo
1052
1053 if (.not.hwrf_samfshal) then
1054 if (do_aerosols) then
1055 kk = itc -3
1056 else
1057 kk = ntr
1058 endif
1059 do n = 1, kk
1060 do k = 2, km1
1061 do i = 1, im
1062 if (cnvflg(i)) then
1063 if(k > kb(i) .and. k < kmax(i)) then
1064 dz = zi(i,k) - zi(i,k-1)
1065 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1066 tem = cq * tem
1067 factor = 1. + tem
1068 ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem*
1069 & (ctro(i,k,n)+ctro(i,k-1,n)))/factor
1070 ercko(i,k,n) = ecko(i,k,n)
1071 endif
1072 endif
1073 enddo
1074 enddo
1075 enddo
1076 if (do_aerosols) then
1077 do n = 1, ntc
1078 kk = n + itc -3
1079 do k = 2, km1
1080 do i = 1, im
1081 if (cnvflg(i)) then
1082 if(k > kb(i) .and. k < kmax(i)) then
1083 dz = zi(i,k) - zi(i,k-1)
1084 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1085 tem = cq * tem
1086 factor = 1. + tem
1087 ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *
1088 & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
1089 ercko(i,k,kk) = ecko(i,k,kk)
1090 chem_c(i,k,n) = escav * fscav(n) * ecko(i,k,kk)
1091 tem = chem_c(i,k,n) / (1. + c0t(i,k) * dz)
1092 chem_pw(i,k,n) = c0t(i,k) * dz * tem * eta(i,k-1)
1093 ecko(i,k,kk) = tem + ecko(i,k,kk) - chem_c(i,k,n)
1094 endif
1095 endif
1096 enddo
1097 enddo
1098 enddo
1099 if(ntk > 2) then
1100 kk = ntk -2
1101 do k = 2, km1
1102 do i = 1, im
1103 if (cnvflg(i)) then
1104 if(k > kb(i) .and. k < kmax(i)) then
1105 dz = zi(i,k) - zi(i,k-1)
1106 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1107 tem = cq * tem
1108 factor = 1. + tem
1109 ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *
1110 & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
1111 ercko(i,k,kk) = ecko(i,k,kk)
1112 endif
1113 endif
1114 enddo
1115 enddo
1116 endif
1117 endif
1118 endif
1119c
1120c taking account into convection inhibition due to existence of
1121c dry layers below cloud base
1122c
1124 do i=1,im
1125 flg(i) = cnvflg(i)
1126 kbcon1(i) = kmax(i)
1127 enddo
1128 do k = 2, km1
1129 do i=1,im
1130 if (flg(i) .and. k < kbm(i)) then
1131 if(k >= kbcon(i) .and. dbyo(i,k) > 0.) then
1132 kbcon1(i) = k
1133 flg(i) = .false.
1134 endif
1135 endif
1136 enddo
1137 enddo
1138 do i=1,im
1139 if(cnvflg(i)) then
1140 if(kbcon1(i) == kmax(i)) cnvflg(i) = .false.
1141 endif
1142 enddo
1143 do i=1,im
1144 if(cnvflg(i)) then
1145 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
1146 if(tem > dthk) then
1147 cnvflg(i) = .false.
1148 endif
1149 endif
1150 enddo
1151!!
1152 totflg = .true.
1153 do i = 1, im
1154 totflg = totflg .and. (.not. cnvflg(i))
1155 enddo
1156 if(totflg) return
1157!!
1158c
1159c calculate convective inhibition
1160c
1162 do k = 2, km1
1163 do i = 1, im
1164 if (cnvflg(i)) then
1165 if(k > kb(i) .and. k < kbcon1(i)) then
1166 dz1 = zo(i,k+1) - zo(i,k)
1167 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1168 rfact = 1. + fv * cp * gamma
1169 & * to(i,k) / hvap
1170 cina(i) = cina(i) +
1171! & dz1 * eta(i,k) * (grav / (cp * to(i,k)))
1172 & dz1 * (grav / (cp * to(i,k)))
1173 & * dbyo(i,k) / (1. + gamma)
1174 & * rfact
1175 val = 0.
1176 cina(i) = cina(i) +
1177! & dz1 * eta(i,k) * grav * fv *
1178 & dz1 * grav * fv *
1179 & max(val,(qeso(i,k) - qo(i,k)))
1180 endif
1181 endif
1182 enddo
1183 enddo
1185
1186 if (hwrf_samfshal) then
1187 do i = 1, im
1188 if(cnvflg(i)) then
1189 cinacr = cinacrmx
1190 if(cina(i) < cinacr) cnvflg(i) = .false.
1191 endif
1192 enddo
1193 else
1194 do i = 1, im
1195 if(cnvflg(i)) then
1196 if(islimsk(i) == 1) then
1197 w1 = w1l
1198 w2 = w2l
1199 w3 = w3l
1200 w4 = w4l
1201 else
1202 w1 = w1s
1203 w2 = w2s
1204 w3 = w3s
1205 w4 = w4s
1206 endif
1207 if(pdot(i) <= w4) then
1208 tem = (pdot(i) - w4) / (w3 - w4)
1209 elseif(pdot(i) >= -w4) then
1210 tem = - (pdot(i) + w4) / (w4 - w3)
1211 else
1212 tem = 0.
1213 endif
1214
1215 val1 = -1.
1216 tem = max(tem,val1)
1217 val2 = 1.
1218 tem = min(tem,val2)
1219 tem = 1. - tem
1220 tem1= .5*(cinacrmx-cinacrmn)
1221 cinacr = cinacrmx - tem * tem1
1222 if(cina(i) < cinacr) cnvflg(i) = .false.
1223 endif
1224 enddo
1225 endif
1226!!
1227 totflg = .true.
1228 do i=1,im
1229 totflg = totflg .and. (.not. cnvflg(i))
1230 enddo
1231 if(totflg) return
1232!!
1233c
1234c determine first guess cloud top as the level of zero buoyancy
1235c limited to the level of p/ps=0.7
1236c
1238 do i = 1, im
1239 flg(i) = cnvflg(i)
1240 if(flg(i)) ktcon(i) = 1
1241 enddo
1242 do k = 2, km1
1243 do i=1,im
1244 if (flg(i) .and. k < kbm(i)) then
1245 if(k > kbcon1(i) .and. dbyo(i,k) < 0.) then
1246 ktcon(i) = k
1247 flg(i) = .false.
1248 endif
1249 endif
1250 enddo
1251 enddo
1252c
1253c turn off shallow convection if cloud depth is larger than cthk or less than cthkmn
1254c
1255 do i = 1, im
1256 if(cnvflg(i)) then
1257 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
1258 if(tem > cthk .or. tem < cthkmn) then
1259 cnvflg(i) = .false.
1260 endif
1261 endif
1262 enddo
1263
1264c
1265c specify upper limit of mass flux at cloud base
1266c
1268 do i = 1, im
1269 if(cnvflg(i)) then
1270 k = kbcon(i)
1271 dp = 1000. * del(i,k)
1272 xmbmax(i) = dp / (grav * dt2)
1273 endif
1274 enddo
1275c
1276c compute cloud moisture property and precipitation
1277c
1279 do i = 1, im
1280 if (cnvflg(i)) then
1281 aa1(i) = 0.
1282 qcko(i,kb(i)) = qo(i,kb(i))
1283 qrcko(i,kb(i)) = qo(i,kb(i))
1284 endif
1285 enddo
1287 do k = 2, km1
1288 do i = 1, im
1289 if (cnvflg(i)) then
1290 if(k > kb(i) .and. k < ktcon(i)) then
1291 dz = zi(i,k) - zi(i,k-1)
1292 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1293 qrch = qeso(i,k)
1294 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1295cj
1296 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1297 tem = cq * tem
1298 factor = 1. + tem
1299 qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem*
1300 & (qo(i,k)+qo(i,k-1)))/factor
1301 qrcko(i,k) = qcko(i,k)
1302cj
1303 dq = eta(i,k) * (qcko(i,k) - qrch)
1304c
1305! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
1306c
1307c below lfc check if there is excess moisture to release latent heat
1308c
1309 if(k >= kbcon(i) .and. dq > 0.) then
1310 etah = .5 * (eta(i,k) + eta(i,k-1))
1311 dp = 1000. * del(i,k)
1312 if(ncloud > 0) then
1313 ptem = c0t(i,k) + c1
1314 qlk = dq / (eta(i,k) + etah * ptem * dz)
1315 dellal(i,k) = etah * c1 * dz * qlk * grav / dp
1316 else
1317 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1318 endif
1319 buo(i,k) = buo(i,k) - grav * qlk
1320 qcko(i,k)= qlk + qrch
1321 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1322 cnvwt(i,k) = etah * qlk * grav / dp
1323 zdqca(i,k)=dq/eta(i,k)
1324 endif
1325!
1326! compute buoyancy and drag for updraft velocity
1327!
1328 if(k >= kbcon(i)) then
1329 rfact = 1. + fv * cp * gamma
1330 & * to(i,k) / hvap
1331 buo(i,k) = buo(i,k) + (grav / (cp * to(i,k)))
1332 & * dbyo(i,k) / (1. + gamma)
1333 & * rfact
1334 val = 0.
1335 buo(i,k) = buo(i,k) + grav * fv *
1336 & max(val,(qeso(i,k) - qo(i,k)))
1337 drag(i,k) = max(xlamue(i,k),xlamud(i))
1338!
1339 tem = ((uo(i,k)-uo(i,k-1))/dz)**2
1340 tem = tem+((vo(i,k)-vo(i,k-1))/dz)**2
1341 wush(i,k) = csmf * sqrt(tem)
1342!
1343 endif
1344!
1345 endif
1346 endif
1347 enddo
1348 enddo
1349c
1350c calculate cloud work function
1351c
1352! do k = 2, km1
1353! do i = 1, im
1354! if (cnvflg(i)) then
1355! if(k >= kbcon(i) .and. k < ktcon(i)) then
1356! dz1 = zo(i,k+1) - zo(i,k)
1357! gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1358! rfact = 1. + fv * cp * gamma
1359! & * to(i,k) / hvap
1360! aa1(i) = aa1(i) +
1361!! & dz1 * eta(i,k) * (grav / (cp * to(i,k)))
1362! & dz1 * (grav / (cp * to(i,k)))
1363! & * dbyo(i,k) / (1. + gamma)
1364! & * rfact
1365! val = 0.
1366! aa1(i) = aa1(i) +
1367!! & dz1 * eta(i,k) * grav * fv *
1368! & dz1 * grav * fv *
1369! & max(val,(qeso(i,k) - qo(i,k)))
1370! endif
1371! endif
1372! enddo
1373! enddo
1374! do i = 1, im
1375! if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
1376! enddo
1377!
1378! calculate cloud work function
1379!
1385 do i = 1, im
1386 if (cnvflg(i)) then
1387 aa1(i) = 0.
1388 endif
1389 enddo
1390 do k = 2, km1
1391 do i = 1, im
1392 if (cnvflg(i)) then
1393 if(k >= kbcon(i) .and. k < ktcon(i)) then
1394 dz1 = zo(i,k+1) - zo(i,k)
1395 aa1(i) = aa1(i) + buo(i,k) * dz1
1396 dbyo1(i,k) = hcko(i,k) - heso(i,k)
1397 endif
1398 endif
1399 enddo
1400 enddo
1401 do i = 1, im
1402 if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
1403 enddo
1404!!
1406 totflg = .true.
1407 do i=1,im
1408 totflg = totflg .and. (.not. cnvflg(i))
1409 enddo
1410 if(totflg) return
1411!!
1412c
1413c estimate the onvective overshooting as the level
1414c where the [aafac * cloud work function] becomes zero,
1415c which is the final cloud top
1416c limited to the level of p/ps=0.7
1417c
1419 do i = 1, im
1420 if (cnvflg(i)) then
1421 aa1(i) = aafac * aa1(i)
1422 endif
1423 enddo
1424c
1425 do i = 1, im
1426 flg(i) = cnvflg(i)
1427 ktcon1(i) = kbm(i)
1428 enddo
1429 do k = 2, km1
1430 do i = 1, im
1431 if (flg(i)) then
1432 if(k >= ktcon(i) .and. k < kbm(i)) then
1433 dz1 = zo(i,k+1) - zo(i,k)
1434 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1435 rfact = 1. + fv * cp * gamma
1436 & * to(i,k) / hvap
1437 aa1(i) = aa1(i) +
1438! & dz1 * eta(i,k) * (grav / (cp * to(i,k)))
1439 & dz1 * (grav / (cp * to(i,k)))
1440 & * dbyo(i,k) / (1. + gamma)
1441 & * rfact
1442! val = 0.
1443! aa1(i) = aa1(i) +
1444!! & dz1 * eta(i,k) * grav * fv *
1445! & dz1 * grav * fv *
1446! & max(val,(qeso(i,k) - qo(i,k)))
1447 if(aa1(i) < 0.) then
1448 ktcon1(i) = k
1449 flg(i) = .false.
1450 endif
1451 endif
1452 endif
1453 enddo
1454 enddo
1455c
1456c compute cloud moisture property, detraining cloud water
1457c and precipitation in overshooting layers
1458c
1460 do k = 2, km1
1461 do i = 1, im
1462 if (cnvflg(i)) then
1463 if(k >= ktcon(i) .and. k < ktcon1(i)) then
1464 dz = zi(i,k) - zi(i,k-1)
1465 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1466 qrch = qeso(i,k)
1467 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1468cj
1469 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1470 tem = cq * tem
1471 factor = 1. + tem
1472 qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem*
1473 & (qo(i,k)+qo(i,k-1)))/factor
1474 qrcko(i,k) = qcko(i,k)
1475cj
1476 dq = eta(i,k) * (qcko(i,k) - qrch)
1477c
1478c check if there is excess moisture to release latent heat
1479c
1480 if(dq > 0.) then
1481 etah = .5 * (eta(i,k) + eta(i,k-1))
1482 dp = 1000. * del(i,k)
1483 if(ncloud > 0) then
1484 ptem = c0t(i,k) + c1
1485 qlk = dq / (eta(i,k) + etah * ptem * dz)
1486 dellal(i,k) = etah * c1 * dz * qlk * grav / dp
1487 else
1488 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1489 endif
1490 qcko(i,k) = qlk + qrch
1491 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1492 cnvwt(i,k) = etah * qlk * grav / dp
1493 zdqca(i,k)=dq/eta(i,k)
1494 endif
1495 endif
1496 endif
1497 enddo
1498 enddo
1499!
1500! compute updraft velocity square(wu2)
1503 if (hwrf_samfshal) then
1504 do i = 1, im
1505 if (cnvflg(i)) then
1506 k = kbcon1(i)
1507 tem = po(i,k) / (rd * to(i,k))
1508 wucb = -0.01 * dot(i,k) / (tem * grav)
1509 if(wucb > 0.) then
1510 wu2(i,k) = wucb * wucb
1511 else
1512 wu2(i,k) = 0.
1513 endif
1514 endif
1515 enddo
1516 endif
1517!
1518 if (progomega) then
1519 call progomega_calc(first_time_step,restart,im,km,
1520 & kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
1521 & grav,buo,drag,wush,xlamue,bb1,bb2)
1522 do k = 1, km
1523 do i = 1, im
1524 if (cnvflg(i)) then
1525 if(k > kbcon1(i) .and. k < ktcon(i)) then
1526 omega_u(i,k)=omegaout(i,k)
1527 omega_u(i,k)=max(omega_u(i,k),-80.)
1528! Convert to m/s for use in convective time-scale:
1529 rho = po(i,k)*100. / (rd * to(i,k))
1530 tem = (-omega_u(i,k)) / ((rho * grav))
1531 wu2(i,k) = tem**2
1532 wu2(i,k) = max(wu2(i,k), 0.)
1533 endif
1534 endif
1535 enddo
1536 enddo
1537
1538 else
1539! diagnostic updraft velocity
1540 do k = 2, km1
1541 do i = 1, im
1542 if (cnvflg(i)) then
1543 if(k > kbcon1(i) .and. k < ktcon(i)) then
1544 dz = zi(i,k) - zi(i,k-1)
1545 tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
1546 tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
1547 tem2 = wush(i,k) * sqrt(wu2(i,k-1))
1548 tem2 = (tem1 - tem2) * dz
1549 ptem = (1. - tem) * wu2(i,k-1)
1550 ptem1 = 1. + tem
1551 wu2(i,k) = (ptem + tem2) / ptem1
1552 wu2(i,k) = max(wu2(i,k), 0.)
1553 endif
1554 endif
1555 enddo
1556 enddo
1557!convert to Pa/s for use in closure
1558 do k = 2, km1
1559 do i = 1, im
1560 if (cnvflg(i)) then
1561 if(k > kbcon1(i) .and. k < ktcon(i)) then
1562 rho = po(i,k)*100. / (rd * to(i,k))
1563 omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav
1564 omega_u(i,k)=max(omega_u(i,k),-80.)
1565 endif
1566 endif
1567 enddo
1568 enddo
1569
1570 endif !progomega
1571
1572! compute updraft velocity averaged over the whole cumulus
1573!
1575 do i = 1, im
1576 wc(i) = 0.
1577 sumx(i) = 0.
1578 enddo
1579 do k = 2, km1
1580 do i = 1, im
1581 if (cnvflg(i)) then
1582 if(k > kbcon1(i) .and. k < ktcon(i)) then
1583 dz = zi(i,k) - zi(i,k-1)
1584 tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1585 wc(i) = wc(i) + tem * dz
1586 sumx(i) = sumx(i) + dz
1587 endif
1588 endif
1589 enddo
1590 enddo
1591 do i = 1, im
1592 if(cnvflg(i)) then
1593 if(sumx(i) == 0.) then
1594 cnvflg(i)=.false.
1595 else
1596 wc(i) = wc(i) / sumx(i)
1597 endif
1598 val = 1.e-4
1599 if (wc(i) < val) cnvflg(i)=.false.
1600 endif
1601 enddo
1602c
1604 if(progsigma)then
1605 do i = 1, im
1606 omegac(i) = 0.
1607 sumx(i) = 0.
1608 enddo
1609 do k = 2, km1
1610 do i = 1, im
1611 if (cnvflg(i)) then
1612 if(k > kbcon1(i) .and. k < ktcon(i)) then
1613 dp = 1000. * del(i,k)
1614 tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
1615 omegac(i) = omegac(i) + tem * dp
1616 sumx(i) = sumx(i) + dp
1617 endif
1618 endif
1619 enddo
1620 enddo
1621 do i = 1, im
1622 if(cnvflg(i)) then
1623 if(sumx(i) == 0.) then
1624 cnvflg(i)=.false.
1625 else
1626 omegac(i) = omegac(i) / sumx(i)
1627 endif
1628 val = -1.2
1629 if (omegac(i) > val) cnvflg(i)=.false.
1630 endif
1631 enddo
1632
1634 do k = 2, km1
1635 do i = 1, im
1636 if (cnvflg(i)) then
1637 if(k > kbcon1(i) .and. k < ktcon(i)) then
1638 if(omega_u(i,k) .ne. 0.)then
1639 zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k))
1640 else
1641 zeta(i,k)=0.
1642 endif
1643 zeta(i,k)=max(0.,zeta(i,k))
1644 zeta(i,k)=min(1.,zeta(i,k))
1645 endif
1646 endif
1647 enddo
1648 enddo
1649 endif !if progsigma
1650
1651c exchange ktcon with ktcon1
1652c
1653 do i = 1, im
1654 if(cnvflg(i)) then
1655 kk = ktcon(i)
1656 ktcon(i) = ktcon1(i)
1657 ktcon1(i) = kk
1658 endif
1659 enddo
1660c
1661c this section is ready for cloud water
1662c
1663 if(ncloud > 0) then
1664c
1665c compute liquid and vapor separation at cloud top
1666c
1668 do i = 1, im
1669 if(cnvflg(i)) then
1670 k = ktcon(i) - 1
1671 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1672 qrch = qeso(i,k)
1673 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1674 dq = qcko(i,k) - qrch
1675c
1676c check if there is excess moisture to release latent heat
1677c
1678 if(dq > 0.) then
1679 qlko_ktcon(i) = dq
1680 qcko(i,k) = qrch
1681 zdqca(i,k) = dq
1682 endif
1683 endif
1684 enddo
1685 endif
1686c
1687
1688c--- compute precipitation efficiency in terms of windshear
1689c
1690!! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 :
1691!! \f[
1692!! 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
1693!! \f]
1694!! 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$.
1695! do i = 1, im
1696! if(cnvflg(i)) then
1697! vshear(i) = 0.
1698! endif
1699! enddo
1700! do k = 2, km
1701! do i = 1, im
1702! if (cnvflg(i)) then
1703! if(k > kb(i) .and. k <= ktcon(i)) then
1704! shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
1705! & + (vo(i,k)-vo(i,k-1)) ** 2)
1706! vshear(i) = vshear(i) + shear
1707! endif
1708! endif
1709! enddo
1710! enddo
1711! do i = 1, im
1712! if(cnvflg(i)) then
1713! vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
1714! e1=1.591-.639*vshear(i)
1715! & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1716! edt(i)=1.-e1
1717! val = .9
1718! edt(i) = min(edt(i),val)
1719! val = .0
1720! edt(i) = max(edt(i),val)
1721! endif
1722! enddo
1723c
1724c--- what would the change be, that a cloud with unit mass
1725c--- will do to the environment?
1726c
1729 do k = 1, km
1730 do i = 1, im
1731 if(cnvflg(i) .and. k <= kmax(i)) then
1732 dellah(i,k) = 0.
1733 dellaq(i,k) = 0.
1734 dellau(i,k) = 0.
1735 dellav(i,k) = 0.
1736 endif
1737 enddo
1738 enddo
1739 if (.not.hwrf_samfshal) then
1740 do n = 1, ntr
1741 do k = 1, km
1742 do i = 1, im
1743 if(cnvflg(i) .and. k <= kmax(i)) then
1744 dellae(i,k,n) = 0.
1745 endif
1746 enddo
1747 enddo
1748 enddo
1749 endif
1750c
1751c--- changed due to subsidence and entrainment
1752c
1753 do k = 2, km1
1754 do i = 1, im
1755 if (cnvflg(i)) then
1756 if(k > kb(i) .and. k < ktcon(i)) then
1757 dp = 1000. * del(i,k)
1758 dz = zi(i,k) - zi(i,k-1)
1759c
1760 dv1h = heo(i,k)
1761 dv2h = .5 * (heo(i,k) + heo(i,k-1))
1762 dv3h = heo(i,k-1)
1763c
1764 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1765 tem1 = xlamud(i)
1766
1767 factor = grav / dp
1768cj
1769 dellah(i,k) = dellah(i,k) +
1770 & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h
1771 & - tem*eta(i,k-1)*dv2h*dz
1772 & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1773 & ) * factor
1774cj
1775 tem1 = -eta(i,k) * qrcko(i,k)
1776 tem2 = -eta(i,k-1) * qcko(i,k-1)
1777 dellaq(i,k) = dellaq(i,k) + (tem1-tem2) * factor
1778cj
1779 tem1=eta(i,k)*(uo(i,k)-ucko(i,k))
1780 tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1))
1781 dellau(i,k) = dellau(i,k) + (tem1-tem2) * factor
1782cj
1783 tem1=eta(i,k)*(vo(i,k)-vcko(i,k))
1784 tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1))
1785 dellav(i,k) = dellav(i,k) + (tem1-tem2) * factor
1786cj
1787 endif
1788 endif
1789 enddo
1790 enddo
1791 if(.not.hwrf_samfshal) then
1792 do n = 1, ntr
1793 do k = 2, km1
1794 do i = 1, im
1795 if (cnvflg(i)) then
1796 if(k > kb(i) .and. k < ktcon(i)) then
1797 dp = 1000. * del(i,k)
1798cj
1799 tem1 = -eta(i,k) * ercko(i,k,n)
1800 tem2 = -eta(i,k-1) * ecko(i,k-1,n)
1801 dellae(i,k,n) = dellae(i,k,n) + (tem1-tem2) * grav/dp
1802cj
1803 endif
1804 endif
1805 enddo
1806 enddo
1807 enddo
1808 endif
1809c
1810c------- cloud top
1811c
1812 do i = 1, im
1813 if(cnvflg(i)) then
1814 indx = ktcon(i)
1815 dp = 1000. * del(i,indx)
1816 tem = eta(i,indx-1) * grav / dp
1817 dellah(i,indx) = tem * (hcko(i,indx-1) - heo(i,indx-1))
1818 dellaq(i,indx) = tem * qcko(i,indx-1)
1819 dellau(i,indx) = tem * (ucko(i,indx-1) - uo(i,indx-1))
1820 dellav(i,indx) = tem * (vcko(i,indx-1) - vo(i,indx-1))
1821c
1822c cloud water
1823c
1824 dellal(i,indx) = tem * qlko_ktcon(i)
1825 endif
1826 enddo
1827 if (.not.hwrf_samfshal) then
1828 do n = 1, ntr
1829 do i = 1, im
1830 if(cnvflg(i)) then
1831 indx = ktcon(i)
1832 dp = 1000. * del(i,indx)
1833 dellae(i,indx,n) = eta(i,indx-1) *
1834 & ecko(i,indx-1,n) * grav / dp
1835 endif
1836 enddo
1837 enddo
1838 endif
1839!
1840! compute change rates due to environmental subsidence & uplift
1841! using a positive definite TVD flux-limiter scheme
1842!
1843! for moisture
1844!
1845 do k=1,km1
1846 do i=1,im
1847 if(cnvflg(i) .and. k <= ktcon(i)) then
1848 q_diff(i,k) = q1(i,k) - q1(i,k+1)
1849 endif
1850 enddo
1851 enddo
1852 do i=1,im
1853 if(cnvflg(i)) then
1854 if(q1(i,1) >= 0.) then
1855 q_diff(i,0) = max(0.,2.*q1(i,1)-q1(i,2))-
1856 & q1(i,1)
1857 else
1858 q_diff(i,0) = min(0.,2.*q1(i,1)-q1(i,2))-
1859 & q1(i,1)
1860 endif
1861 endif
1862 enddo
1863!
1864 flxtvd = 0.
1865 do k = 1, km1
1866 do i = 1, im
1867 if(cnvflg(i) .and.
1868 & (k >= kb(i) .and. k < ktcon(i))) then
1869 if(eta(i,k) > 0.) then
1870 rrkp = 0.
1871 if(abs(q_diff(i,k)) > 1.e-22)
1872 & rrkp = q_diff(i,k+1) / q_diff(i,k)
1873 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1874 tem1 = q1(i,k+1) +
1875 & phkp*(qo(i,k)-q1(i,k+1))
1876 flxtvd(i,k) = eta(i,k) * tem1
1877 endif
1878 endif
1879 enddo
1880 enddo
1881!
1882 do k = 2, km1
1883 do i = 1, im
1884 if(cnvflg(i) .and.
1885 & (k > kb(i) .and. k <= ktcon(i))) then
1886 dp = 1000. * del(i,k)
1887 dellaq(i,k) = dellaq(i,k) +
1888 & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp
1889 endif
1890 enddo
1891 enddo
1892!
1893! for tracers including TKE & ozone
1894!
1895 if (.not.hwrf_samfshal) then
1896!
1897 do n=1,ntr
1898 do k=1,km1
1899 do i=1,im
1900 if(cnvflg(i) .and. k <= ktcon(i)) then
1901 e_diff(i,k,n) = ctr(i,k,n) - ctr(i,k+1,n)
1902 endif
1903 enddo
1904 enddo
1905 do i=1,im
1906 if(cnvflg(i)) then
1907 if(ctr(i,1,n) >= 0.) then
1908 e_diff(i,0,n) = max(0.,2.*ctr(i,1,n)-ctr(i,2,n))-
1909 & ctr(i,1,n)
1910 else
1911 e_diff(i,0,n) = min(0.,2.*ctr(i,1,n)-ctr(i,2,n))-
1912 & ctr(i,1,n)
1913 endif
1914 endif
1915 enddo
1916 enddo
1917!
1918 do n=1,ntr
1919!
1920 flxtvd = 0.
1921 do k= 1, km1
1922 do i = 1, im
1923 if(cnvflg(i) .and.
1924 & (k >= kb(i) .and. k < ktcon(i))) then
1925 if(eta(i,k) > 0.) then
1926 rrkp = 0.
1927 if(abs(e_diff(i,k,n)) > 1.e-22)
1928 & rrkp = e_diff(i,k+1,n) / e_diff(i,k,n)
1929 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1930 tem1 = ctr(i,k+1,n) +
1931 & phkp*(ctro(i,k,n)-ctr(i,k+1,n))
1932 flxtvd(i,k) = eta(i,k) * tem1
1933 endif
1934 endif
1935 enddo
1936 enddo
1937!
1938 do k = 2, km1
1939 do i = 1, im
1940 if(cnvflg(i) .and.
1941 & (k > kb(i) .and. k <= ktcon(i))) then
1942 dp = 1000. * del(i,k)
1943 dellae(i,k,n) = dellae(i,k,n) +
1944 & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp
1945 endif
1946 enddo
1947 enddo
1948!
1949 enddo
1950!
1951 endif
1952!
1953! compute convective turn-over time
1954!
1956 do i= 1, im
1957 if(cnvflg(i)) then
1958 tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
1959 dtconv(i) = tem / wc(i)
1960 if (.not.hwrf_samfshal) then
1961 tfac = 1. + gdx(i) / 75000.
1962 dtconv(i) = tfac * dtconv(i)
1963 endif
1964 dtconv(i) = max(dtconv(i),dtmin)
1965 dtconv(i) = max(dtconv(i),dt2)
1966 dtconv(i) = min(dtconv(i),dtmax)
1967 endif
1968 enddo
1969!
1971 do i= 1, im
1972 if(cnvflg(i)) then
1973 sumx(i) = 0.
1974 umean(i) = 0.
1975 endif
1976 enddo
1977 do k = 2, km1
1978 do i = 1, im
1979 if(cnvflg(i)) then
1980 if(k >= kbcon1(i) .and. k < ktcon1(i)) then
1981 dz = zi(i,k) - zi(i,k-1)
1982 tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
1983 umean(i) = umean(i) + tem * dz
1984 sumx(i) = sumx(i) + dz
1985 endif
1986 endif
1987 enddo
1988 enddo
1989 do i= 1, im
1990 if(cnvflg(i)) then
1991 umean(i) = umean(i) / sumx(i)
1992 umean(i) = max(umean(i), 1.)
1993 tauadv = gdx(i) / umean(i)
1994 advfac(i) = tauadv / dtconv(i)
1995 advfac(i) = min(advfac(i), 1.)
1996 endif
1997 enddo
1998c
1999c compute cloud base mass flux as a function of the mean
2000c updraft velcoity
2001c
2003 if(progsigma)then
2004! Initial computations, dynamic q-tendency
2005 if(first_time_step .and. .not.restart)then
2006 do k = 1,km
2007 do i = 1,im
2008 qadv(i,k)=0.
2009 enddo
2010 enddo
2011 else
2012 do k = 1,km
2013 do i = 1,im
2014 qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
2015 enddo
2016 enddo
2017 endif
2018
2019 do k = 1,km
2020 do i = 1,im
2021 tmfq(i,k)=tmf(i,k,1)
2022 enddo
2023 enddo
2024
2025 flag_shallow = .true.
2026 flag_mid = .false.
2027 call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
2028 & flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,
2029 & delt,qadv,kb,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
2030 & sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
2031 endif
2032
2035
2036 do i= 1, im
2037 if(cnvflg(i)) then
2038 k = kbcon(i)
2039 rho = po(i,k)*100. / (rd*to(i,k))
2040 if(progsigma .and. gdx(i) < dxcrtas)then
2041 xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv)
2042 else
2043 xmb(i) = advfac(i)*betaw*rho*wc(i)
2044 endif
2045 endif
2046 enddo
2047!
2049 do i = 1, im
2050 if(cnvflg(i)) then
2051 tem = min(max(xlamue(i,kbcon(i)), 2.e-4), 6.e-4)
2052 tem = 0.2 / tem
2053 tem1 = 3.14 * tem * tem
2054 sigmagfm(i) = tem1 / garea(i)
2055 sigmagfm(i) = max(sigmagfm(i), 0.001)
2056 sigmagfm(i) = min(sigmagfm(i), 0.999)
2057 endif
2058 enddo
2059!
2061 do i = 1, im
2062 if(cnvflg(i)) then
2063 if (gdx(i) < dxcrt) then
2064 if(progsigma)then
2065 scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i))
2066 else
2067 scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
2068 endif
2069 scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
2070 else
2071 scaldfunc(i) = 1.0
2072 endif
2073 xmb(i) = xmb(i) * scaldfunc(i)
2074 xmb(i) = min(xmb(i),xmbmax(i))
2075 endif
2076 enddo
2077!
2079!
2080! if (.not.hwrf_samfshal) then
2081! if (do_aerosols)
2082! & call samfshalcnv_aerosols(im, im, km, itc, ntc, ntr, delt,
2083!! & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kbcon, ktcon, fscav,
2084! & cnvflg, kb, kmax, ktcon, fscav,
2085!! & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp,
2086! & xmb, c0t, eta, zi, xlamue, xlamud, delp,
2087! & qtr, qaero)
2088! endif
2089!
2092c
2093c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2094c
2095 do k = 1, km
2096 do i = 1, im
2097 if (cnvflg(i) .and. k <= kmax(i)) then
2098 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
2099 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
2100 val = 1.e-8
2101 qeso(i,k) = max(qeso(i,k), val )
2102 endif
2103 enddo
2104 enddo
2105c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2106c
2110 do i = 1, im
2111 delhbar(i) = 0.
2112 delqbar(i) = 0.
2113 deltbar(i) = 0.
2114 delubar(i) = 0.
2115 delvbar(i) = 0.
2116 qcond(i) = 0.
2117 enddo
2118 if (.not. hwrf_samfshal) then
2119 do n = 1, ntr
2120 do i = 1, im
2121 delebar(i,n) = 0.
2122 enddo
2123 enddo
2124 endif
2125 do k = 1, km
2126 do i = 1, im
2127 if (cnvflg(i)) then
2128 if(k > kb(i) .and. k <= ktcon(i)) then
2129 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
2130 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
2131 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
2132! tem = 1./rcs(i)
2133! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
2134! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
2135 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
2136 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
2137 dp = 1000. * del(i,k)
2138 tem = xmb(i) * dp / grav
2139 delhbar(i) = delhbar(i) + tem * dellah(i,k)
2140 delqbar(i) = delqbar(i) + tem * dellaq(i,k)
2141 deltbar(i) = deltbar(i) + tem * dellat
2142 delubar(i) = delubar(i) + tem * dellau(i,k)
2143 delvbar(i) = delvbar(i) + tem * dellav(i,k)
2144 endif
2145 endif
2146 enddo
2147 enddo
2148!
2149! Negative moisture is set to zero after borrowing it from
2150! positive values within the mass-flux transport layers
2151!
2152 do i = 1,im
2153 tsumn(i) = 0.
2154 tsump(i) = 0.
2155 rtnp(i) = 1.
2156 enddo
2157 do k = 1,km1
2158 do i = 1,im
2159 if (cnvflg(i)) then
2160 if(k > kb(i) .and. k <= ktcon(i)) then
2161 tem = q1(i,k) * delp(i,k) / grav
2162 if(q1(i,k) < 0.) tsumn(i) = tsumn(i) + tem
2163 if(q1(i,k) > 0.) tsump(i) = tsump(i) + tem
2164 endif
2165 endif
2166 enddo
2167 enddo
2168 do i = 1,im
2169 if(cnvflg(i)) then
2170 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2171 if(tsump(i) > abs(tsumn(i))) then
2172 rtnp(i) = tsumn(i) / tsump(i)
2173 else
2174 rtnp(i) = tsump(i) / tsumn(i)
2175 endif
2176 endif
2177 endif
2178 enddo
2179 do k = 1,km1
2180 do i = 1,im
2181 if (cnvflg(i)) then
2182 if(k > kb(i) .and. k <= ktcon(i)) then
2183 if(rtnp(i) < 0.) then
2184 if(tsump(i) > abs(tsumn(i))) then
2185 if(q1(i,k) < 0.) q1(i,k)= 0.
2186 if(q1(i,k) > 0.) q1(i,k)=(1.+rtnp(i))*q1(i,k)
2187 else
2188 if(q1(i,k) < 0.) q1(i,k)=(1.+rtnp(i))*q1(i,k)
2189 if(q1(i,k) > 0.) q1(i,k)=0.
2190 endif
2191 endif
2192 endif
2193 endif
2194 enddo
2195 enddo
2196!
2197 if (.not.hwrf_samfshal) then
2198!
2199 indx = ntk - 2
2200 do n = 1, ntr
2201!
2202 do k = 1, km
2203 do i = 1, im
2204 if (cnvflg(i)) then
2205 if(k > kb(i) .and. k <= ktcon(i)) then
2206 ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2
2207 dp = 1000. * del(i,k)
2208 delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/grav
2209 endif
2210 endif
2211 enddo
2212 enddo
2213!
2214! Negative TKE, ozone, and aerosols are set to zero after borrowing them
2215! from positive values within the mass-flux transport layers
2216!
2217 do i = 1,im
2218 tsumn(i) = 0.
2219 tsump(i) = 0.
2220 rtnp(i) = 1.
2221 enddo
2222 do k = 1,km1
2223 do i = 1,im
2224 if (cnvflg(i)) then
2225 if(k > kb(i) .and. k <= ktcon(i)) then
2226 if(n == indx) then
2227 if(k > 1) then
2228 dz = zi(i,k) - zi(i,k-1)
2229 else
2230 dz = zi(i,k)
2231 endif
2232 tem = ctr(i,k,n) * dz
2233 else
2234 tem = ctr(i,k,n) * delp(i,k) / grav
2235 endif
2236 if(ctr(i,k,n) < 0.) tsumn(i) = tsumn(i) + tem
2237 if(ctr(i,k,n) > 0.) tsump(i) = tsump(i) + tem
2238 endif
2239 endif
2240 enddo
2241 enddo
2242 do i = 1,im
2243 if(cnvflg(i)) then
2244 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2245 if(tsump(i) > abs(tsumn(i))) then
2246 rtnp(i) = tsumn(i) / tsump(i)
2247 else
2248 rtnp(i) = tsump(i) / tsumn(i)
2249 endif
2250 endif
2251 endif
2252 enddo
2253 do k = 1,km1
2254 do i = 1,im
2255 if (cnvflg(i)) then
2256 if(k > kb(i) .and. k <= ktcon(i)) then
2257 if(rtnp(i) < 0.) then
2258 if(tsump(i) > abs(tsumn(i))) then
2259 if(ctr(i,k,n)<0.) ctr(i,k,n)=0.
2260 if(ctr(i,k,n)>0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n)
2261 else
2262 if(ctr(i,k,n)<0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n)
2263 if(ctr(i,k,n)>0.) ctr(i,k,n)=0.
2264 endif
2265 endif
2266 endif
2267 endif
2268 enddo
2269 enddo
2270!
2271 kk = n+2
2272 do k = 1, km
2273 do i = 1, im
2274 if (cnvflg(i)) then
2275 if(k > kb(i) .and. k <= ktcon(i)) then
2276 qtr(i,k,kk) = ctr(i,k,n)
2277 endif
2278 endif
2279 enddo
2280 enddo
2281!
2282 enddo
2283!
2284 if (do_aerosols) then
2285!
2286 do n = 1, ntc
2287!
2288! convert wet deposition to total mass deposited over dt2 and dp
2289 do k = 2, km1
2290 do i = 1, im
2291 if (cnvflg(i)) then
2292 if(k > kb(i) .and. k < ktcon(i)) then
2293 dp = 1000. * del(i,k)
2294 wet_dep(i,k,n) = chem_pw(i,k,n)*grav/dp
2295 wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*dt2*dp
2296 endif
2297 endif
2298 enddo
2299 enddo
2300!
2301 kk = n + itc - 1
2302 do k = 2, km1
2303 do i = 1, im
2304 if (cnvflg(i)) then
2305 if(k > kb(i) .and. k < ktcon(i)) then
2306 dp = 1000. * del(i,k)
2307 if (qtr(i,k,kk) < 0.) then
2308! borrow negative mass from wet deposition
2309 tem = -qtr(i,k,kk)*dp
2310 if(wet_dep(i,k,n) >= tem) then
2311 wet_dep(i,k,n) = wet_dep(i,k,n) - tem
2312 qtr(i,k,kk) = 0.
2313 else
2314 wet_dep(i,k,n) = 0.
2315 qtr(i,k,kk) = qtr(i,k,kk)+wet_dep(i,k,n)/dp
2316 endif
2317 endif
2318 endif
2319 endif
2320 enddo
2321 enddo
2322!
2323 enddo
2324!
2325 endif
2326!
2327 endif
2328!
2330 do k = 1, km
2331 do i = 1, im
2332 if (cnvflg(i)) then
2333 if(k > kb(i) .and. k <= ktcon(i)) then
2334 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
2335 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
2336 val = 1.e-8
2337 qeso(i,k) = max(qeso(i,k), val )
2338 endif
2339 endif
2340 enddo
2341 enddo
2342c
2344 do i = 1, im
2345 rntot(i) = 0.
2346 delqev(i) = 0.
2347 delq2(i) = 0.
2348 flg(i) = cnvflg(i)
2349 enddo
2350 do k = km, 1, -1
2351 do i = 1, im
2352 if (cnvflg(i)) then
2353 if(k < ktcon(i) .and. k > kb(i)) then
2354 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
2355 endif
2356 endif
2357 enddo
2358 enddo
2359c
2360c evaporating rain
2361c
2365 do k = km, 1, -1
2366 do i = 1, im
2367 if (k <= kmax(i)) then
2368 deltv(i) = 0.
2369 delq(i) = 0.
2370 qevap(i) = 0.
2371 if(cnvflg(i)) then
2372 if(k < ktcon(i) .and. k > kb(i)) then
2373 rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
2374 endif
2375 endif
2376 if(flg(i) .and. k < ktcon(i)) then
2377! evef = edt(i) * evfact
2378! if(islimsk(i) == 1) evef=edt(i) * evfactl
2379! if(islimsk(i) == 1) evef=.07
2380 qcond(i) = shevf * evef * (q1(i,k) - qeso(i,k))
2381 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
2382 dp = 1000. * del(i,k)
2383 factor = dp / grav
2384 if(rn(i) > 0. .and. qcond(i) < 0.) then
2385 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
2386 qevap(i) = min(qevap(i), rn(i)*1000.*grav/dp)
2387 delq2(i) = delqev(i) + .001 * qevap(i) * factor
2388 endif
2389 if(rn(i) > 0. .and. qcond(i) < 0. .and.
2390 & delq2(i) > rntot(i)) then
2391 qevap(i) = 1000.* grav * (rntot(i) - delqev(i)) / dp
2392 flg(i) = .false.
2393 endif
2394 if(rn(i) > 0. .and. qevap(i) > 0.) then
2395 tem = .001 * factor
2396 tem1 = qevap(i) * tem
2397 if(tem1 > rn(i)) then
2398 qevap(i) = rn(i) / tem
2399 rn(i) = 0.
2400 else
2401 rn(i) = rn(i) - tem1
2402 endif
2403 q1(i,k) = q1(i,k) + qevap(i)
2404 t1(i,k) = t1(i,k) - elocp * qevap(i)
2405 deltv(i) = - elocp*qevap(i)/dt2
2406 delq(i) = + qevap(i)/dt2
2407 delqev(i) = delqev(i) + tem * qevap(i)
2408 endif
2409 delqbar(i) = delqbar(i) + delq(i) * factor
2410 deltbar(i) = deltbar(i) + deltv(i) * factor
2411 endif
2412 endif
2413 enddo
2414 enddo
2415cj
2416! do i = 1, im
2417! if(me == 31 .and. cnvflg(i)) then
2418! if(cnvflg(i)) then
2419! print *, ' shallow delhbar, delqbar, deltbar = ',
2420! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
2421! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
2422! print *, ' precip =', hvap*rn(i)*1000./dt2
2423! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
2424! endif
2425! enddo
2426! do n = 1, ntr
2427! do i = 1, im
2428! if(me == 31 .and. cnvflg(i)) then
2429! if(cnvflg(i)) then
2430! print *, ' tracer delebar = ',delebar(i,n)
2431! endif
2432! enddo
2433! enddo
2434cj
2435 do i = 1, im
2436 if(cnvflg(i)) then
2437 if(rn(i) < 0. .or. .not.flg(i)) rn(i) = 0.
2438 ktop(i) = ktcon(i)
2439 kbot(i) = kbcon(i)
2440 kcnv(i) = 2
2441 endif
2442 enddo
2443c
2444c convective cloud water
2445 do k = 1, km
2446 do i = 1, im
2447 if (cnvflg(i)) then
2448 if (k >= kbcon(i) .and. k < ktcon(i)) then
2449 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
2450 if (progsigma) then
2451 cnvw(i,k) = cnvw(i,k) * cscale
2452 else
2453 cnvw(i,k) = cnvw(i,k) * cscale
2454 endif
2455 endif
2456 endif
2457 enddo
2458 enddo
2459c
2460c convective cloud cover
2461c
2463 do k = 1, km
2464 do i = 1, im
2465 if (cnvflg(i)) then
2466 if (k >= kbcon(i) .and. k < ktcon(i)) then
2467 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
2468 cnvc(i,k) = min(cnvc(i,k), 0.2)
2469 cnvc(i,k) = max(cnvc(i,k), 0.0)
2470 endif
2471 endif
2472 enddo
2473 enddo
2474c
2475c cloud water
2476c
2478 if (ncloud > 0) then
2479!
2480 do k = 1, km1
2481 do i = 1, im
2482 if (cnvflg(i)) then
2483! if (k > kb(i) .and. k <= ktcon(i)) then
2484 if (k >= kbcon(i) .and. k <= ktcon(i)) then
2485 tem = dellal(i,k) * xmb(i) * dt2
2486 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2487 if (qtr(i,k,2) > -999.0) then
2488 qtr(i,k,1) = qtr(i,k,1) + tem * tem1 ! ice
2489 qtr(i,k,2) = qtr(i,k,2) + tem *(1.0-tem1) ! water
2490 else
2491 qtr(i,k,1) = qtr(i,k,1) + tem
2492 endif
2493 endif
2494 endif
2495 enddo
2496 enddo
2497!
2498 endif
2500! if (.not. hwrf_samfshal) then
2501! if (do_aerosols) then
2502! do n = 1, ntc
2503! kk = n + itc - 1
2504! do k = 1, km
2505! do i = 1, im
2506! if(cnvflg(i) .and. rn(i) > 0.) then
2507! if (k <= kmax(i)) qtr(i,k,kk) = qaero(i,k,n)
2508! endif
2509! enddo
2510! enddo
2511! enddo
2512! endif
2513! endif
2514!
2515! hchuang code change
2516!
2518!
2520 do k = 1, km
2521 do i = 1, im
2522 if(cnvflg(i)) then
2523 if(k >= kb(i) .and. k < ktop(i)) then
2524 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
2525 endif
2526 endif
2527 enddo
2528 enddo
2530 do i = 1, im
2531 if(cnvflg(i)) then
2532 k = ktop(i)-1
2533 dt_mf(i,k) = ud_mf(i,k)
2534 endif
2535 enddo
2536!
2537! include TKE contribution from shallow convection
2538!
2539 if (.not.hwrf_samfshal) then
2540 if (ntk > 0) then
2541!
2542 do k = 2, km1
2543 do i = 1, im
2544 if(cnvflg(i)) then
2545 if(k > kb(i) .and. k < ktop(i)) then
2546 tem = 0.5 * (eta(i,k-1) + eta(i,k)) * xmb(i)
2547 tem1 = pfld(i,k) * 100. / (rd * t1(i,k))
2548 if(progsigma)then
2549 tem2 = sigmab(i)
2550 else
2551 tem2 = max(sigmagfm(i), betaw)
2552 endif
2553 ptem = tem / (tem2 * tem1)
2554 qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem
2555 endif
2556 endif
2557 enddo
2558 enddo
2559!
2560 endif
2561 endif
2562!!
2563 return
2564 end subroutine samfshalcnv_run
2566 end module samfshalcnv
2567
energy
subroutine energy(parameters, ice,vegtyp,ist,nsnow,nsoil, isnow,dt,rhoair,sfcprs,qair, sfctmp,thair,lwdn,uu,vv,zref, co2air,o2air,solad,solai,cosz,igs, eair,tbot,zsnso,zsoil, elai,esai,fwet,foln, fveg,shdfac, pahv,pahg,pahb, qsnow,dzsnso,lat,canliq,canice,iloc, jloc, thsfc_loc, prslkix, prsik1x, prslk1x, garea1, pblhx, iz0tlnd, itime, psi_opt, ep_1, ep_2, epsm1, cp, z0wrf,z0hwrf, imelt,snicev,snliqv,epore,t2m,fsno, sav,sag,qmelt,fsa,fsr,taux, tauy,fira,fsh,fcev,fgev,fctr, trad,psn,apar,ssoil,btrani,btran, ponding, ts,latheav, latheag, frozen_canopy, frozen_ground, tv,tg,stc,snowh,eah,tah, sneqvo,sneqv,sh2o,smc,snice,snliq, albold,cm,ch,dx,dz8w,q2, ustarx, ifdef ccpp
We use different approaches to deal with subgrid features of radiation transfer and turbulent transfe...
Definition module_sf_noahmplsm.F90:1693
water
subroutine water(parameters, vegtyp, nsnow, nsoil, imelt, dt, uu, vv, fcev, fctr, qprecc, qprecl, elai, esai, sfctmp, qvap, qdew, zsoil, btrani, ficeold, ponding, tg, ist, fveg, iloc, jloc, smceq, bdfall, fp, rain, snow, qsnow, qrain, snowhin, latheav, latheag, frozen_canopy, frozen_ground, isnow, canliq, canice, tv, snowh, sneqv, snice, snliq, stc, zsnso, sh2o, smc, sice, zwt, wa, wt, dzsnso, wslake, smcwtd, deeprech, rech, cmc, ecan, etran, fwet, runsrf, runsub, qin, qdis, ponding1, ponding2, qsnbot, esnow)
compute water budgets (water storages, et components, and runoff)
Definition module_sf_noahmplsm.F90:7047
progomega::progomega_calc
subroutine, public progomega_calc(first_time_step, flag_restart, im, km, kbcon1, ktcon, omegain, delt, del, zi, cnvflg, omegaout, grav, buo, drag, wush, tentr, bb1, bb2)
This subroutine computes a prognostic updraft vertical velocity used in the closure computations in t...
Definition progomega_calc.f90:24
samfshalcnv::samfshalcnv_run
subroutine samfshalcnv_run(im, km, itc, ntc, cliq, cp, cvap, eps, epsm1, fv, grav, hvap, rd, rv, t0c, delt, ntk, ntr, delp, first_time_step, restart, tmf, qmicro, progsigma, progomega, prslp, psp, phil, tkeh, qtr, prevsq, q, q1, t1, u1, v1, fscav, rn, kbot, ktop, kcnv, islimsk, garea, cscale, dot, ncloud, hpbl, ud_mf, dt_mf, cnvw, cnvc, clam, c0s, c1, evef, pgcon, asolfac, hwrf_samfshal, sigmain, sigmaout, omegain, omegaout, betadcu, betamcu, betascu, errmsg, errflg)
Definition samfshalcnv.f:63
progsigma
This module contains the subroutine that calculates the prognostic updraft area fraction that is used...
Definition progsigma_calc.f90:7
samfshalcnv
This module contains the Scale-Aware mass flux Shallow Convection scheme.
Definition samfshalcnv.f:5