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