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