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