CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
samfdeepcnv.f
1
4
7 module samfdeepcnv
8
9 use samfcnv_aerosols, only : samfdeepcnv_aerosols
10 use progsigma, only : progsigma_calc
11
12 contains
13
14 subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, &
15 & errmsg, errflg)
16
17 integer, intent(in) :: imfdeepcnv
18 integer, intent(in) :: imfdeepcnv_samf
19 character(len=*), intent(out) :: errmsg
20 integer, intent(out) :: errflg
21
22
23 ! Consistency checks
24 if (imfdeepcnv/=imfdeepcnv_samf) then
25 write(errmsg,'(*(a))') 'Logic error: namelist choice of', &
26 & ' deep convection is different from SAMF scheme'
27 errflg = 1
28 return
29 end if
30
31 end subroutine samfdeepcnv_init
32
75 subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
76 & tmf,qmicro,itc,ntc,cliq,cp,cvap, &
77 & eps,epsm1,fv,grav,hvap,rd,rv, &
78 & t0c,delt,ntk,ntr,delp, &
79 & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, &
80 & hwrf_samfdeep,progsigma,cldwrk,rn,kbot,ktop,kcnv, &
81 & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, &
82 & QLCN, QICN, w_upi, cf_upi, CNV_MFD, &
83 & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,&
84 & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, &
85 & do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep, &
86 & rainevap,sigmain,sigmaout,betadcu,betamcu,betascu, &
87 & maxMF, do_mynnedmf,errmsg,errflg)
88!
89 use machine , only : kind_phys
90 use funcphys , only : fpvs
91
92 implicit none
93!
94 integer, intent(in) :: im, km, itc, ntc, ntk, ntr, ncloud
95 integer, intent(in) :: islimsk(:)
96 real(kind=kind_phys), intent(in) :: cliq, cp, cvap, eps, epsm1, &
97 & fv, grav, hvap, rd, rv, t0c
98 real(kind=kind_phys), intent(in) :: delt
99 real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), &
100 & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:)
101 real(kind=kind_phys), dimension(:), intent(in) :: fscav
102 logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, &
103 & progsigma,do_mynnedmf
104 real(kind=kind_phys), intent(in) :: nthresh,betadcu,betamcu, &
105 & betascu
106 real(kind=kind_phys), intent(in), optional :: ca_deep(:)
107 real(kind=kind_phys), intent(in), optional :: sigmain(:,:), &
108 & qmicro(:,:), prevsq(:,:)
109 real(kind=kind_phys), intent(in) :: tmf(:,:,:),q(:,:)
110 real(kind=kind_phys), dimension (:), intent(in), optional :: maxmf
111 real(kind=kind_phys), intent(out) :: rainevap(:)
112 real(kind=kind_phys), intent(out), optional :: sigmaout(:,:)
113 logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger
114 integer, intent(inout) :: kcnv(:)
115 ! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH
116 real(kind=kind_phys), intent(inout) :: qtr(:,:,:), &
117 & q1(:,:), t1(:,:), u1(:,:), v1(:,:), &
118 & cnvw(:,:), cnvc(:,:)
119
120 integer, intent(out) :: kbot(:), ktop(:)
121 real(kind=kind_phys), intent(out) :: cldwrk(:), &
122 & rn(:), &
123 & dd_mf(:,:), dt_mf(:,:)
124 real(kind=kind_phys), intent(out), optional :: ud_mf(:,:)
125 ! GJF* These variables are conditionally allocated depending on whether the
126 ! Morrison-Gettelman microphysics is used, so they must be declared
127 ! using assumed shape.
128 real(kind=kind_phys), dimension(:,:), intent(inout), optional :: &
129 & qlcn, qicn, w_upi, cnv_mfd, cnv_dqldt, clcn &
130 &, cnv_fice, cnv_ndrop, cnv_nice, cf_upi
131 ! *GJF
132 integer, intent(in) :: mp_phys, mp_phys_mg
133
134 real(kind=kind_phys), intent(in) :: clam, c0s, c1, &
135 & betal, betas, asolfac, &
136 & evef, pgcon
137 character(len=*), intent(out) :: errmsg
138 integer, intent(out) :: errflg
139!
140!------local variables
141 integer i, indx, jmn, k, kk, km1, n
142! integer latd,lond
143!
144 real(kind=kind_phys) clamd, tkemx, tkemn, dtke,
145 & beta, clamca,
146 & cxlame, cxlamd,
147 & xlamde, xlamdd,
148 & crtlame, crtlamd
149!
150! real(kind=kind_phys) detad
151 real(kind=kind_phys) adw, aup, aafac, d0,
152 & dellat, desdt, dg,
153 & dh, dhh, dp,
154 & dq, dqsdp, dqsdt, dt,
155 & dt2, dtmax, dtmin,
156 & dxcrtas, dxcrtuf,
157 & dv1h, dv2h, dv3h,
158 & dz, dz1, e1, edtmax,
159 & edtmaxl, edtmaxs, el2orc, elocp,
160 & es, etah,
161 & cthk, dthk,
162! & evfact, evfactl,
163 & fact1, fact2, factor,
164 & gamma, pprime, cm, cq,
165 & qlk, qrch, qs,
166 & rain, rfact, shear, tfac,
167 & val, val1, val2,
168 & w1, w1l, w1s, w2,
169 & w2l, w2s, w3, w3l,
170 & w3s, w4, w4l, w4s,
171 & rho, betaw, tauadv,
172 & xdby, xpw, xpwd,
173! & xqrch, mbdt, tem,
174 & xqrch, tem, tem1, tem2,
175 & ptem, ptem1, ptem2
176!
177 integer kb(im), kb1(im), kbcon(im), kbcon1(im),
178 & ktcon(im), ktcon1(im), ktconn(im),
179 & jmin(im), lmin(im), kbmax(im),
180 & kbm(im), kmax(im), kd94(im)
181!
182! real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im),
183 real(kind=kind_phys) aa1(im), tkemean(im),clamt(im),
184 & ps(im), del(im,km), prsl(im,km),
185 & umean(im), advfac(im), gdx(im),
186 & delhbar(im), delq(im), delq2(im),
187 & delqbar(im), delqev(im), deltbar(im),
188 & deltv(im), dtconv(im), edt(im),
189 & edto(im), edtx(im), fld(im),
190 & hcdo(im,km), hmax(im), hmin(im),
191 & ucdo(im,km), vcdo(im,km),aa2(im),
192 & ecdo(im,km,ntr),
193 & pdot(im), po(im,km),
194 & pwavo(im), pwevo(im), mbdt(im),
195 & qcdo(im,km), qcond(im), qevap(im),
196 & rntot(im), vshear(im), xaa0(im),
197 & xlamd(im), xlamdet(im),xlamddt(im),
198 & cxlamet(im), cxlamdt(im),
199 & xk(im), cina(im),
200 & xmb(im), xmbmax(im), xpwav(im),
201 & xpwev(im), xlamx(im), delebar(im,ntr),
202! & xpwev(im), delebar(im,ntr),
203 & delubar(im), delvbar(im)
204!
205 real(kind=kind_phys) c0(im), sfcpbl(im)
206cj
207 real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn,
208 & cinacr, cinacrmx, cinacrmn,
209 & sfclfac, rhcrt,
210 & tkcrt, cmxfac
211cj
212!
213! parameters for updraft velocity calculation
214 real(kind=kind_phys) bb1, bb2, csmf, wucb
215!
216! parameters for prognostic sigma closure
217 real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
218 & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km)
219 real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
220 parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
221 logical flag_shallow, flag_mid
222c physical parameters
223! parameter(grav=grav,asolfac=0.958)
224! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
225! parameter(c0s=.002,c1=.002,d0=.01)
226! parameter(d0=.01)
227 parameter(d0=.001)
228! parameter(c0l=c0s*asolfac)
229!
230! asolfac: aerosol-aware parameter based on Lim (2011)
231! asolfac= cx / c0s(=.002)
232! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s)
233! Nccn: CCN number concentration in cm^(-3)
234! Until a realistic Nccn is provided, Nccns are assumed
235! as Nccn=100 for sea and Nccn=1000 for land
236!
237 parameter(cm=1.0,cq=1.0)
238! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
239 parameter(clamd=0.03,tkemx=0.65,tkemn=0.05)
240 parameter(clamca=0.03)
241 parameter(dtke=tkemx-tkemn)
242 parameter(cthk=200.,dthk=25.,sfclfac=0.2,rhcrt=0.75)
243 parameter(cinpcrmx=180.,cinpcrmn=120.)
244! parameter(cinacrmx=-120.,cinacrmn=-120.)
245 parameter(cinacrmx=-120.,cinacrmn=-80.)
246 parameter(bb1=4.0,bb2=0.8,csmf=0.2)
247 parameter(tkcrt=2.,cmxfac=15.)
248 parameter(betaw=.03)
249
250!
251! local variables and arrays
252 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
253 & uo(im,km), vo(im,km), qeso(im,km),
254 & ctr(im,km,ntr), ctro(im,km,ntr)
255! for aerosol transport
256! real(kind=kind_phys) qaero(im,km,ntc)
257c variables for tracer wet deposition,
258 real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw,
259 & wet_dep
260!
261! for updraft velocity calculation
262 real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km),
263 & wush(im,km), wc(im)
264!
265! for updraft fraction & scale-aware function
266 real(kind=kind_phys) scaldfunc(im), sigmagfm(im)
267!
268c cloud water
269! real(kind=kind_phys) tvo(im,km)
270 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
271 & dbyo(im,km), zo(im,km),
272 & xlamue(im,km), xlamud(im,km),
273 & fent1(im,km), fent2(im,km),
274 & rh(im,km), frh(im,km),
275 & heo(im,km), heso(im,km),
276 & qrcd(im,km), dellah(im,km), dellaq(im,km),
277 & dellae(im,km,ntr),
278 & dellau(im,km), dellav(im,km), hcko(im,km),
279 & ucko(im,km), vcko(im,km), qcko(im,km),
280 & ecko(im,km,ntr),ercko(im,km,ntr),
281 & eta(im,km), etad(im,km), zi(im,km),
282 & qrcko(im,km), qrcdo(im,km),
283 & pwo(im,km), pwdo(im,km), c0t(im,km),
284 & tx1(im), sumx(im), cnvwt(im,km)
285 &, rhbar(im)
286!
287! variables for Total Variation Diminishing (TVD) flux-limiter scheme
288! on environmental subsidence and uplifting
289!
290 real(kind=kind_phys) q_diff(im,0:km-1), e_diff(im,0:km-1,ntr),
291 & flxtvd(im,0:km-1)
292 real(kind=kind_phys) rrkp, phkp
293 real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im)
294!
295 logical do_aerosols, totflg, cnvflg(im), asqecflg(im), flg(im)
296!
297! asqecflg: flag for the quasi-equilibrium assumption of Arakawa-Schubert
298!
299! real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
300!! save pcrit, acritt
301! data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,
302! & 350.,300.,250.,200.,150./
303! data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,
304! & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
305c gdas derived acrit
306c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
307c & .743,.813,.886,.947,1.138,1.377,1.896/
308 real(kind=kind_phys) tf, tcr, tcrf
309 parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
310
311 ! Initialize CCPP error handling variables
312 errmsg = ''
313 errflg = 0
314
315 gravinv = 1./grav
316 invdelt = 1./delt
317
318 elocp = hvap/cp
319 el2orc = hvap*hvap/(rv*cp)
320
321 fact1 = (cvap-cliq)/rv
322 fact2 = hvap/rv-fact1*t0c
323c-----------------------------------------------------------------------
325 if(hwrf_samfdeep) then
326 do_aerosols = .false.
327 else
328 do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0)
329 if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3)
330 endif
331!
332c-----------------------------------------------------------------------
335!************************************************************************
336! convert input Pa terms to Cb terms -- Moorthi
337 ps = psp * 0.001
338 prsl = prslp * 0.001
339 del = delp * 0.001
340!************************************************************************
341!
342!
343 km1 = km - 1
345c
346c initialize arrays
347c
348 chem_c = 0.
349 chem_pw = 0.
350 wet_dep = 0.
351!
352 do i=1,im
353 cnvflg(i) = .true.
354 if(do_mynnedmf) then
355 if(maxmf(i).gt.0.)cnvflg(i)=.false.
356 endif
357 sfcpbl(i) = sfclfac * hpbl(i)
358 rn(i)=0.
359 mbdt(i)=10.
360 kbot(i)=km+1
361 ktop(i)=0
362 kbcon(i)=km
363 ktcon(i)=1
364 ktconn(i)=1
365 dtconv(i) = 3600.
366 cldwrk(i) = 0.
367 pdot(i) = 0.
368 lmin(i) = 1
369 jmin(i) = 1
370 qlko_ktcon(i) = 0.
371 edt(i) = 0.
372 edto(i) = 0.
373 edtx(i) = 0.
374! acrt(i) = 0.
375! acrtfct(i) = 1.
376 aa1(i) = 0.
377 aa2(i) = 0.
378 xaa0(i) = 0.
379 cina(i) = 0.
380 pwavo(i)= 0.
381 pwevo(i)= 0.
382 xmb(i) = 0.
383 xpwav(i)= 0.
384 xpwev(i)= 0.
385 vshear(i) = 0.
386 advfac(i) = 0.
387 rainevap(i) = 0.
388 omegac(i)=0.
389 gdx(i) = sqrt(garea(i))
390 enddo
391
392 do k=1,km
393 do i=1,im
394 xlamud(i,k) = 0.
395 xlamue(i,k) = 0.
396 enddo
397 enddo
398!
399 if (hwrf_samfdeep) then
400 do i=1,im
401 scaldfunc(i)=-1.0
402 sigmagfm(i)=-1.0
403! sigmuout(i)=-1.0
404 enddo
405 endif
406!
408 do i=1,im
409 if(islimsk(i) == 1) then
410 c0(i) = c0s*asolfac
411 else
412 c0(i) = c0s
413 endif
414 enddo
415!
417 do k = 1, km
418 do i = 1, im
419 if(t1(i,k) > 273.16) then
420 c0t(i,k) = c0(i)
421 else
422 tem = d0 * (t1(i,k) - 273.16)
423 tem1 = exp(tem)
424 c0t(i,k) = c0(i) * tem1
425 endif
426 enddo
427 enddo
429 do k = 1, km
430 do i = 1, im
431 cnvw(i,k) = 0.
432 cnvc(i,k) = 0.
433 enddo
434 enddo
435! hchuang code change
437 do k = 1, km
438 do i = 1, im
439 ud_mf(i,k) = 0.
440 dd_mf(i,k) = 0.
441 dt_mf(i,k) = 0.
442 enddo
443 enddo
444 if(mp_phys == mp_phys_mg) then
445 do k = 1, km
446 do i = 1, im
447 qlcn(i,k) = qtr(i,k,2)
448 qicn(i,k) = qtr(i,k,1)
449 w_upi(i,k) = 0.0
450 cf_upi(i,k) = 0.0
451 cnv_mfd(i,k) = 0.0
452
453 cnv_dqldt(i,k) = 0.0
454 clcn(i,k) = 0.0
455 cnv_fice(i,k) = 0.0
456 cnv_ndrop(i,k) = 0.0
457 cnv_nice(i,k) = 0.0
458 enddo
459 enddo
460 endif
461c
462! do k = 1, 15
463! acrit(k) = acritt(k) * (975. - pcrit(k))
464! enddo
465!
466 dt2 = delt
467! val = 1200.
468 val = 600.
469 dtmin = max(dt2, val )
470! val = 5400.
471 val = 10800.
472 dtmax = max(dt2, val )
473! model tunable parameters are all here
474 edtmaxl = .3
475 edtmaxs = .3
476! evfact = 0.3
477! evfactl = 0.3
478 aafac = .1
479 if (hwrf_samfdeep) then
480 cxlame = 1.0e-3
481 else
482 cxlame = 1.0e-4
483 endif
484 cxlamd = 0.75e-4
485 crtlame = 1.0e-4
486 crtlamd = 1.0e-4
487 xlamde = 1.0e-4
488 xlamdd = 1.0e-4
489!
490! pgcon = 0.7 ! Gregory et al. (1997, QJRMS)
491! pgcon = 0.55 ! Zhang & Wu (2003,JAS)
492!
493 w1l = -8.e-3
494 w2l = -4.e-2
495 w3l = -5.e-3
496 w4l = -5.e-4
497 w1s = -2.e-4
498 w2s = -2.e-3
499 w3s = -1.e-3
500 w4s = -2.e-5
501c
502c define top layer for search of the downdraft originating layer
503c and the maximum thetae for updraft
504c
506 do i=1,im
507 kbmax(i) = km
508 kbm(i) = km
509 kmax(i) = km
510 kd94(i) = km
511 tx1(i) = 1.0 / ps(i)
512 enddo
513!
514 do k = 1, km
515 do i=1,im
516 if (prsl(i,k)*tx1(i) > 0.04) kmax(i) = k + 1
517 if (prsl(i,k)*tx1(i) > 0.45) kbmax(i) = k + 1
518 if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1
519 if (prsl(i,k)*tx1(i) > 0.94) kd94(i) = k + 1
520 enddo
521 enddo
522 do i=1,im
523 kmax(i) = min(km,kmax(i))
524 kbmax(i) = min(kbmax(i),kmax(i))
525 kbm(i) = min(kbm(i),kmax(i))
526 kd94(i) = min(kd94(i),kmax(i))
527 enddo
528c
529c hydrostatic height assume zero terr and initially assume
530c updraft entrainment rate as an inverse function of height
531c
533 do k = 1, km
534 do i=1,im
535 zo(i,k) = phil(i,k) / grav
536 enddo
537 enddo
539 do k = 1, km1
540 do i=1,im
541 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
542 enddo
543 enddo
544 if (hwrf_samfdeep) then
545 do k = 1, km1
546 do i=1,im
547 xlamue(i,k) = clam / zi(i,k)
548 enddo
549 enddo
550 endif
551c
552c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
553c convert surface pressure to mb from cb
554c
556 do k = 1, km
557 do i = 1, im
558 if (k <= kmax(i)) then
559 pfld(i,k) = prsl(i,k) * 10.0
560 eta(i,k) = 1.
561 fent1(i,k)= 1.
562 fent2(i,k)= 1.
563 rh(i,k) = 0.
564 frh(i,k) = 0.
565 hcko(i,k) = 0.
566 qcko(i,k) = 0.
567 qrcko(i,k)= 0.
568 ucko(i,k) = 0.
569 vcko(i,k) = 0.
570 etad(i,k) = 1.
571 hcdo(i,k) = 0.
572 qcdo(i,k) = 0.
573 ucdo(i,k) = 0.
574 vcdo(i,k) = 0.
575 qrcd(i,k) = 0.
576 qrcdo(i,k)= 0.
577 dbyo(i,k) = 0.
578 pwo(i,k) = 0.
579 pwdo(i,k) = 0.
580 dellal(i,k) = 0.
581 to(i,k) = t1(i,k)
582 qo(i,k) = q1(i,k)
583 uo(i,k) = u1(i,k)
584 vo(i,k) = v1(i,k)
585! uo(i,k) = u1(i,k) * rcs(i)
586! vo(i,k) = v1(i,k) * rcs(i)
587 wu2(i,k) = 0.
588 buo(i,k) = 0.
589 wush(i,k) = 0.
590 drag(i,k) = 0.
591 cnvwt(i,k)= 0.
592 endif
593 enddo
594 enddo
595
596 do k = 1, km
597 do i = 1, im
598 dbyo1(i,k)=0.
599 zdqca(i,k)=0.
600 omega_u(i,k)=0.
601 zeta(i,k)=1.0
602 enddo
603 enddo
604
605!
606! initialize tracer variables
607!
608 if(.not.hwrf_samfdeep) then
609 do n = 3, ntr+2
610 kk = n-2
611 do k = 1, km
612 do i = 1, im
613 if (k <= kmax(i)) then
614 ctr(i,k,kk) = qtr(i,k,n)
615 ctro(i,k,kk) = qtr(i,k,n)
616 ecko(i,k,kk) = 0.
617 ercko(i,k,kk) = 0.
618 ecdo(i,k,kk) = 0.
619 endif
620 enddo
621 enddo
622 enddo
623 endif
624!
626 do k = 1, km
627 do i=1,im
628 if (k <= kmax(i)) then
629 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
630 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
631 val1 = 1.e-8
632 qeso(i,k) = max(qeso(i,k), val1)
633 val2 = 1.e-10
634 qo(i,k) = max(qo(i,k), val2 )
635! qo(i,k) = min(qo(i,k),qeso(i,k))
636! tvo(i,k) = to(i,k) + fv * to(i,k) * qo(i,k)
637 endif
638 enddo
639 enddo
640c
641c compute moist static energy
642c
644 do k = 1, km
645 do i=1,im
646 if (k <= kmax(i)) then
647! tem = grav * zo(i,k) + cp * to(i,k)
648 tem = phil(i,k) + cp * to(i,k)
649 heo(i,k) = tem + hvap * qo(i,k)
650 heso(i,k) = tem + hvap * qeso(i,k)
651c heo(i,k) = min(heo(i,k),heso(i,k))
652 endif
653 enddo
654 enddo
655c
656c determine level with largest moist static energy
657c this is the level where updraft starts
658c
661 do i=1,im
662 flg(i) = .true.
663 kb1(i) = 1
664 enddo
665 do k = 2, km1
666 do i=1,im
667 if (flg(i) .and. zo(i,k) <= sfcpbl(i)) then
668 kb1(i) = k
669 else
670 flg(i) = .false.
671 endif
672 enddo
673 enddo
674 do i=1,im
675 kb1(i) = min(kb1(i),kbm(i))
676 enddo
677c
679 do i=1,im
680 hmax(i) = heo(i,kb1(i))
681 kb(i) = kb1(i)
682 enddo
683 do k = 2, km
684 do i=1,im
685 if (k > kb1(i) .and. k <= kbm(i)) then
686 if(heo(i,k) > hmax(i)) then
687 kb(i) = k
688 hmax(i) = heo(i,k)
689 endif
690 endif
691 enddo
692 enddo
693c
695 do k = 1, km1
696 do i=1,im
697 if (k <= kmax(i)-1) then
698 dz = .5 * (zo(i,k+1) - zo(i,k))
699 dp = .5 * (pfld(i,k+1) - pfld(i,k))
700 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
701 pprime = pfld(i,k+1) + epsm1 * es
702 qs = eps * es / pprime
703 dqsdp = - qs / pprime
704 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
705 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
706 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
707 dt = (grav*dz + hvap*dqsdp*dp) / (cp * (1. + gamma))
708 dq = dqsdt * dt + dqsdp * dp
709 to(i,k) = to(i,k+1) + dt
710 qo(i,k) = qo(i,k+1) + dq
711 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
712 endif
713 enddo
714 enddo
715!
717 do k = 1, km1
718 do i=1,im
719 if (k <= kmax(i)-1) then
720 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
721 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
722 val1 = 1.e-8
723 qeso(i,k) = max(qeso(i,k), val1)
724 val2 = 1.e-10
725 qo(i,k) = max(qo(i,k), val2 )
726! qo(i,k) = min(qo(i,k),qeso(i,k))
727 rh(i,k) = min(qo(i,k)/qeso(i,k), 1.)
728 frh(i,k) = 1. - rh(i,k)
729 heo(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
730 & cp * to(i,k) + hvap * qo(i,k)
731 heso(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
732 & cp * to(i,k) + hvap * qeso(i,k)
733 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
734 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
735 endif
736 enddo
737 enddo
738 if (.not.hwrf_samfdeep) then
739 do n = 1, ntr
740 do k = 1, km1
741 do i=1,im
742 if (k <= kmax(i)-1) then
743 ctro(i,k,n) = .5 * (ctro(i,k,n) + ctro(i,k+1,n))
744 endif
745 enddo
746 enddo
747 enddo
748 endif
749c
750c look for the level of free convection as cloud base
751c
753 do i=1,im
754 flg(i) = .true.
755 kbcon(i) = kmax(i)
756 enddo
757 do k = 1, km1
758 do i=1,im
759 if (flg(i) .and. k <= kbmax(i)) then
760 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
761 kbcon(i) = k
762 flg(i) = .false.
763 endif
764 endif
765 enddo
766 enddo
767c
769 do i=1,im
770 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
771 enddo
772!!
773 totflg = .true.
774 do i=1,im
775 totflg = totflg .and. (.not. cnvflg(i))
776 enddo
777 if(totflg) return
778!!
780 do i=1,im
781 if(cnvflg(i)) then
782! pdot(i) = 10.* dot(i,kbcon(i))
783 pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
784 endif
785 enddo
786c
787c turn off convection if pressure depth between parcel source level
788c and cloud base is larger than a critical value, cinpcr
789c
790 do i=1,im
791 if(cnvflg(i)) then
792 if(islimsk(i) == 1) then
793 w1 = w1l
794 w2 = w2l
795 w3 = w3l
796 w4 = w4l
797 else
798 w1 = w1s
799 w2 = w2s
800 w3 = w3s
801 w4 = w4s
802 endif
803 if(pdot(i) <= w4) then
804 tem = (pdot(i) - w4) / (w3 - w4)
805 elseif(pdot(i) >= -w4) then
806 tem = - (pdot(i) + w4) / (w4 - w3)
807 else
808 tem = 0.
809 endif
810 val1 = -1.
811 tem = max(tem,val1)
812 val2 = 1.
813 tem = min(tem,val2)
814 ptem = 1. - tem
815 ptem1= .5*(cinpcrmx-cinpcrmn)
816 cinpcr = cinpcrmx - ptem * ptem1
817 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
818 if(tem1 > cinpcr .and.
819 & zi(i,kbcon(i)) > hpbl(i)) then
820 cnvflg(i) = .false.
821 endif
822 endif
823 enddo
824!!
825 if(do_ca .and. ca_trigger)then
826 do i=1,im
827 if(ca_deep(i) > nthresh) cnvflg(i) = .true.
828 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
829 enddo
830 endif
831
832 totflg = .true.
833 do i=1,im
834 totflg = totflg .and. (.not. cnvflg(i))
835 enddo
836 if(totflg) return
837!!
838!
839! re-define kb & kbcon
840!
841 do i=1,im
842 if (cnvflg(i)) then
843 hmax(i) = heo(i,1)
844 kb(i) = 1
845 endif
846 enddo
847 do k = 2, km
848 do i=1,im
849 if (cnvflg(i) .and. k <= kbm(i)) then
850 if(heo(i,k) > hmax(i)) then
851 kb(i) = k
852 hmax(i) = heo(i,k)
853 endif
854 endif
855 enddo
856 enddo
857!
858 do i=1,im
859 flg(i) = cnvflg(i)
860 if(flg(i)) kbcon(i) = kmax(i)
861 enddo
862 do k = 1, km1
863 do i=1,im
864 if (flg(i) .and. k <= kbmax(i)) then
865 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
866 kbcon(i) = k
867 flg(i) = .false.
868 endif
869 endif
870 enddo
871 enddo
872!
873 do i=1,im
874 if(cnvflg(i) .and. kbcon(i) == kmax(i)) then
875 cnvflg(i) = .false.
876 endif
877 enddo
878!!
879 if(do_ca .and. ca_trigger)then
880 do i=1,im
881 if(ca_deep(i) > nthresh) cnvflg(i) = .true.
882 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
883 enddo
884 endif
885
886 totflg = .true.
887 do i=1,im
888 totflg = totflg .and. (.not. cnvflg(i))
889 enddo
890 if(totflg) return
891!!
892 do i=1,im
893 if(cnvflg(i)) then
894! pdot(i) = 10.* dot(i,kbcon(i))
895 pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
896 endif
897 enddo
898!
900!
901 do i = 1, im
902 rhbar(i) = 0.
903 sumx(i) = 0.
904 enddo
905 do k = 1, km1
906 do i = 1, im
907 if (cnvflg(i)) then
908 if(k >= kb(i) .and. k < kbcon(i)) then
909 dz = zo(i,k+1) - zo(i,k)
910 rhbar(i) = rhbar(i) + rh(i,k) * dz
911 sumx(i) = sumx(i) + dz
912 endif
913 endif
914 enddo
915 enddo
916 do i= 1, im
917 if(cnvflg(i)) then
918 rhbar(i) = rhbar(i) / sumx(i)
919 if(rhbar(i) < rhcrt) then
920 cnvflg(i) = .false.
921 endif
922 endif
923 enddo
924!!
925 if(do_ca .and. ca_trigger)then
926 do i=1,im
927 if(ca_deep(i) > nthresh) cnvflg(i) = .true.
928 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
929 enddo
930 endif
931
932 totflg = .true.
933 do i=1,im
934 totflg = totflg .and. (.not. cnvflg(i))
935 enddo
936 if(totflg) return
937!!
938
939! turbulent entrainment rate assumed to be proportional
940! to subcloud mean TKE
941!
942 if(.not. hwrf_samfdeep .and. ntk > 0) then
943!
944 do i= 1, im
945 if(cnvflg(i)) then
946 sumx(i) = 0.
947 tkemean(i) = 0.
948 endif
949 enddo
950 do k = 1, km1
951 do i = 1, im
952 if(cnvflg(i)) then
953 if(k >= kb(i) .and. k < kbcon(i)) then
954 dz = zo(i,k+1) - zo(i,k)
955 tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk))
956 tkemean(i) = tkemean(i) + tem * dz
957 sumx(i) = sumx(i) + dz
958 endif
959 endif
960 enddo
961 enddo
962!
963 do i= 1, im
964 if(cnvflg(i)) then
965 tkemean(i) = tkemean(i) / sumx(i)
966 if(tkemean(i) > tkemx) then
967 clamt(i) = clam + clamd
968 else if(tkemean(i) < tkemn) then
969 clamt(i) = clam - clamd
970 else
971 tem = tkemx - tkemean(i)
972 tem1 = 1. - 2. * tem / dtke
973 clamt(i) = clam + clamd * tem1
974 endif
975 endif
976 enddo
977!
978 if(do_ca .and. ca_entr)then
979 do i=1,im
980 if(cnvflg(i)) then
981 if(ca_deep(i) > nthresh)then
982 clamt(i) = clam - clamca
983 else
984 clamt(i) = clam
985 endif
986 endif
987 enddo
988 endif
989!
990 do i=1,im
991 if(cnvflg(i)) then
992 xlamdet(i) = xlamde
993 xlamddt(i) = xlamdd
994 cxlamet(i) = cxlame
995 cxlamdt(i) = cxlamd
996 if(tkemean(i) > tkcrt) then
997 tem = 1. + tkemean(i)/tkcrt
998 tem1 = min(tem, cmxfac)
999 clamt(i) = tem1 * clam
1000 xlamdet(i) = tem1 * xlamdet(i)
1001 xlamddt(i) = tem1 * xlamddt(i)
1002 cxlamet(i) = tem1 * cxlamet(i)
1003 cxlamdt(i) = tem1 * cxlamdt(i)
1004 endif
1005 endif
1006 enddo
1007!
1008 else
1009!
1010 if(do_ca .and. ca_entr)then
1011 do i=1,im
1012 if(cnvflg(i)) then
1013 if(ca_deep(i) > nthresh)then
1014 clamt(i) = clam - clamca
1015 else
1016 clamt(i) = clam
1017 endif
1018 endif
1019 enddo
1020 else
1021 do i=1,im
1022 if(cnvflg(i))then
1023 clamt(i) = clam
1024 endif
1025 enddo
1026 endif
1027!
1028 do i=1,im
1029 if(cnvflg(i)) then
1030 xlamdet(i) = xlamde
1031 xlamddt(i) = xlamdd
1032 cxlamet(i) = cxlame
1033 cxlamdt(i) = cxlamd
1034 endif
1035 enddo
1036!
1037 endif !(.not. hwrf_samfdeep .and. ntk > 0)
1038!
1039! also initially assume updraft entrainment rate
1040! is an inverse function of height
1041!
1042 do k = 1, km1
1043 do i=1,im
1044 if(cnvflg(i)) then
1045 dz =zo(i,k+1) - zo(i,k)
1046 xlamue(i,k) = clamt(i) / (zi(i,k) + dz)
1047 xlamue(i,k) = max(xlamue(i,k), crtlame)
1048 endif
1049 enddo
1050 enddo
1051c
1052c assume that updraft entrainment rate above cloud base is
1053c same as that at cloud base
1054c
1060 if (hwrf_samfdeep) then
1061 do i=1,im
1062 if(cnvflg(i)) then
1063 xlamx(i) = xlamue(i,kbcon(i))
1064 endif
1065 enddo
1066 do k = 2, km1
1067 do i=1,im
1068 if(cnvflg(i).and. &
1069 & (k > kbcon(i) .and. k < kmax(i))) then
1070 xlamue(i,k) = xlamx(i)
1071 endif
1072 enddo
1073 enddo
1074 endif
1075c
1076c specify detrainment rate for the updrafts
1077c
1078!! (The updraft detrainment rate is set constant and equal to the entrainment rate at cloud base.)
1079!!
1081 if (hwrf_samfdeep) then
1082 do k = 1, km1
1083 do i=1,im
1084 if(cnvflg(i) .and. k < kmax(i)) then
1085 xlamud(i,k) = xlamx(i)
1086 endif
1087 enddo
1088 enddo
1089 else
1090 do k = 1, km1
1091 do i=1,im
1092 if(cnvflg(i) .and. k < kmax(i)) then
1093! xlamud(i,k) = crtlamd
1094 xlamud(i,k) = 0.001 * clamt(i)
1095 endif
1096 enddo
1097 enddo
1098 endif
1099c
1100c entrainment functions decreasing with height(fent),
1101c mimicking a cloud ensemble
1102c(bechtold et al., 2008)
1103c
1104 do k = 2, km1
1105 do i=1,im
1106 if(cnvflg(i).and.
1107 & (k > kbcon(i) .and. k < kmax(i))) then
1108 tem = qeso(i,k)/qeso(i,kbcon(i))
1109 fent1(i,k) = tem**2
1110 fent2(i,k) = tem**3
1111 endif
1112 enddo
1113 enddo
1114c
1115c final entrainment and detrainment rates as the sum of turbulent part and
1116c organized one depending on the environmental relative humidity
1117c(bechtold et al., 2008; derbyshire et al., 2011)
1118c
1119 if (hwrf_samfdeep) then
1120 do k = 2, km1
1121 do i=1,im
1122 if(cnvflg(i) .and.
1123 & (k > kbcon(i) .and. k < kmax(i))) then
1124 tem = cxlamet(i) * frh(i,k) * fent2(i,k)
1125 xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
1126 endif
1127 enddo
1128 enddo
1129 else
1130 do k = 2, km1
1131 do i=1,im
1132 if(cnvflg(i) .and.
1133 & (k > kbcon(i) .and. k < kmax(i))) then
1134 tem = cxlamet(i) * frh(i,k) * fent2(i,k)
1135 xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
1136 tem1 = cxlamdt(i) * frh(i,k)
1137 xlamud(i,k) = xlamud(i,k) + tem1
1138 endif
1139 enddo
1140 enddo
1141 endif
1142!
1143!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1144c
1145c determine updraft mass flux for the subcloud layers
1146c
1152 do k = km1, 1, -1
1153 do i = 1, im
1154 if (cnvflg(i)) then
1155 if(k < kbcon(i) .and. k >= kb(i)) then
1156 dz = zi(i,k+1) - zi(i,k)
1157 tem = 0.5*(xlamud(i,k)+xlamud(i,k+1))
1158 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-tem
1159 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
1160 endif
1161 endif
1162 enddo
1163 enddo
1164c
1165c compute mass flux above cloud base
1166c
1167 do i = 1, im
1168 flg(i) = cnvflg(i)
1169 enddo
1170 do k = 2, km1
1171 do i = 1, im
1172 if(flg(i))then
1173 if(k > kbcon(i) .and. k < kmax(i)) then
1174 dz = zi(i,k) - zi(i,k-1)
1175 tem = 0.5*(xlamud(i,k)+xlamud(i,k-1))
1176 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-tem
1177 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
1178 if(eta(i,k) <= 0.) then
1179 kmax(i) = k
1180 ktconn(i) = k
1181 flg(i) = .false.
1182 endif
1183 endif
1184 endif
1185 enddo
1186 enddo
1187c
1188c compute updraft cloud properties
1189c
1191 do i = 1, im
1192 if(cnvflg(i)) then
1193 indx = kb(i)
1194 hcko(i,indx) = heo(i,indx)
1195 ucko(i,indx) = uo(i,indx)
1196 vcko(i,indx) = vo(i,indx)
1197 pwavo(i) = 0.
1198 endif
1199 enddo
1200 if (.not.hwrf_samfdeep) then
1201! for tracers
1202 do n = 1, ntr
1203 do i = 1, im
1204 if(cnvflg(i)) then
1205 indx = kb(i)
1206 ecko(i,indx,n) = ctro(i,indx,n)
1207 ercko(i,indx,n) = ctro(i,indx,n)
1208 endif
1209 enddo
1210 enddo
1211 endif
1212c
1213c cloud property is modified by the entrainment process
1214c
1215! cm is an enhancement factor in entrainment rates for momentum
1216!
1218 do k = 2, km1
1219 do i = 1, im
1220 if (cnvflg(i)) then
1221 if(k > kb(i) .and. k < kmax(i)) then
1222 dz = zi(i,k) - zi(i,k-1)
1223 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1224 tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
1225 factor = 1. + tem - tem1
1226 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
1227 & (heo(i,k)+heo(i,k-1)))/factor
1228 dbyo(i,k) = hcko(i,k) - heso(i,k)
1229!
1230 tem = 0.5 * cm * tem
1231 factor = 1. + tem
1232 ptem = tem + pgcon
1233 ptem1= tem - pgcon
1234 ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*uo(i,k)
1235 & +ptem1*uo(i,k-1))/factor
1236 vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*vo(i,k)
1237 & +ptem1*vo(i,k-1))/factor
1238 endif
1239 endif
1240 enddo
1241 enddo
1242 if (.not.hwrf_samfdeep) then
1243 if (do_aerosols) then
1244 kk = itc -3
1245 else
1246 kk = ntr
1247 endif
1248 do n = 1, kk
1249 do k = 2, km1
1250 do i = 1, im
1251 if (cnvflg(i)) then
1252 if(k > kb(i) .and. k < kmax(i)) then
1253 dz = zi(i,k) - zi(i,k-1)
1254 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1255 tem = cq * tem
1256 factor = 1. + tem
1257 ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem*
1258 & (ctro(i,k,n)+ctro(i,k-1,n)))/factor
1259 ercko(i,k,n) = ecko(i,k,n)
1260 endif
1261 endif
1262 enddo
1263 enddo
1264 enddo
1265 if (do_aerosols) then
1266 do n = 1, ntc
1267 kk = n + itc -3
1268 do k = 2, km1
1269 do i = 1, im
1270 if (cnvflg(i)) then
1271 if(k > kb(i) .and. k < kmax(i)) then
1272 dz = zi(i,k) - zi(i,k-1)
1273 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1274 tem = cq * tem
1275 factor = 1. + tem
1276 ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *
1277 & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
1278 ercko(i,k,kk) = ecko(i,k,kk)
1279 chem_c(i,k,n) = fscav(n) * ecko(i,k,kk)
1280 tem = chem_c(i,k,n) / (1. + c0t(i,k) * dz)
1281 chem_pw(i,k,n) = c0t(i,k) * dz * tem * eta(i,k-1)
1282 ecko(i,k,kk) = tem + ecko(i,k,kk) - chem_c(i,k,n)
1283 endif
1284 endif
1285 enddo
1286 enddo
1287 enddo
1288 endif
1289 endif
1290c
1291c taking account into convection inhibition due to existence of
1292c dry layers below cloud base
1293c
1295 do i=1,im
1296 flg(i) = cnvflg(i)
1297 kbcon1(i) = kmax(i)
1298 enddo
1299 do k = 2, km1
1300 do i=1,im
1301 if (flg(i) .and. k < kmax(i)) then
1302 if(k >= kbcon(i) .and. dbyo(i,k) > 0.) then
1303 kbcon1(i) = k
1304 flg(i) = .false.
1305 endif
1306 endif
1307 enddo
1308 enddo
1309 do i=1,im
1310 if(cnvflg(i)) then
1311 if(kbcon1(i) == kmax(i)) cnvflg(i) = .false.
1312 endif
1313 enddo
1314 do i=1,im
1315 if(cnvflg(i)) then
1316 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
1317 if(tem > dthk) then
1318 cnvflg(i) = .false.
1319 endif
1320 endif
1321 enddo
1322!!
1323
1324 if(do_ca .and. ca_trigger)then
1325 do i=1,im
1326 if(ca_deep(i) > nthresh) cnvflg(i) = .true.
1327 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
1328 enddo
1329 endif
1330
1331 totflg = .true.
1332 do i = 1, im
1333 totflg = totflg .and. (.not. cnvflg(i))
1334 enddo
1335 if(totflg) return
1336!!
1337c
1338c calculate convective inhibition
1339c
1341 do k = 2, km1
1342 do i = 1, im
1343 if (cnvflg(i)) then
1344 if(k > kb(i) .and. k < kbcon1(i)) then
1345 dz1 = zo(i,k+1) - zo(i,k)
1346 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1347 rfact = 1. + fv * cp * gamma
1348 & * to(i,k) / hvap
1349 cina(i) = cina(i) +
1350! & dz1 * eta(i,k) * (grav / (cp * to(i,k)))
1351 & dz1 * (grav / (cp * to(i,k)))
1352 & * dbyo(i,k) / (1. + gamma)
1353 & * rfact
1354 val = 0.
1355 cina(i) = cina(i) +
1356! & dz1 * eta(i,k) * grav * fv *
1357 & dz1 * grav * fv *
1358 & max(val,(qeso(i,k) - qo(i,k)))
1359 endif
1360 endif
1361 enddo
1362 enddo
1364 if(hwrf_samfdeep) then
1365 do i = 1, im
1366 if(cnvflg(i)) then
1367 cinacr = cinacrmx
1368 if(cina(i) < cinacr) cnvflg(i) = .false.
1369 endif
1370 enddo
1371 else !gfs_samfdeep
1372 do i = 1, im
1373 if(cnvflg(i)) then
1374 if(islimsk(i) == 1) then
1375 w1 = w1l
1376 w2 = w2l
1377 w3 = w3l
1378 w4 = w4l
1379 else
1380 w1 = w1s
1381 w2 = w2s
1382 w3 = w3s
1383 w4 = w4s
1384 endif
1385 if(pdot(i) <= w4) then
1386 tem = (pdot(i) - w4) / (w3 - w4)
1387 elseif(pdot(i) >= -w4) then
1388 tem = - (pdot(i) + w4) / (w4 - w3)
1389 else
1390 tem = 0.
1391 endif
1392
1393 val1 = -1.
1394 tem = max(tem,val1)
1395 val2 = 1.
1396 tem = min(tem,val2)
1397 tem = 1. - tem
1398 tem1= .5*(cinacrmx-cinacrmn)
1399 cinacr = cinacrmx - tem * tem1
1400 if(cina(i) < cinacr) cnvflg(i) = .false.
1401 endif
1402 enddo
1403 endif !hwrf_samfdeep
1404!!
1405 if(do_ca .and. ca_trigger)then
1406 do i=1,im
1407 if(ca_deep(i) > nthresh) cnvflg(i) = .true.
1408 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
1409 enddo
1410 endif
1411
1412 totflg = .true.
1413 do i=1,im
1414 totflg = totflg .and. (.not. cnvflg(i))
1415 enddo
1416 if(totflg) return
1417!!
1418c
1419c determine first guess cloud top as the level of zero buoyancy
1420c
1422 do i = 1, im
1423 flg(i) = cnvflg(i)
1424 ktcon(i) = 1
1425 enddo
1426 do k = 2, km1
1427 do i = 1, im
1428 if (flg(i) .and. k < kmax(i)) then
1429 if(k > kbcon1(i) .and. dbyo(i,k) < 0.) then
1430 ktcon(i) = k
1431 flg(i) = .false.
1432 endif
1433 endif
1434 enddo
1435 enddo
1436c
1437 do i = 1, im
1438 if(cnvflg(i)) then
1439 if(ktcon(i) == 1 .and. ktconn(i) > 1) then
1440 ktcon(i) = ktconn(i)
1441 endif
1442 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
1443 if(tem < cthk) cnvflg(i) = .false.
1444 endif
1445 enddo
1446
1447
1448 if(do_ca .and. ca_trigger)then
1449 do i=1,im
1450 if(ca_deep(i) > nthresh) cnvflg(i) = .true.
1451 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
1452 enddo
1453 endif
1454
1455 totflg = .true.
1456 do i=1,im
1457 totflg = totflg .and. (.not. cnvflg(i))
1458 enddo
1459 if(totflg) return
1460!!
1461
1462c
1463c search for downdraft originating level above theta-e minimum
1464c
1466 do i = 1, im
1467 if(cnvflg(i)) then
1468 hmin(i) = heo(i,kbcon1(i))
1469 lmin(i) = kbmax(i)
1470 jmin(i) = kbmax(i)
1471 endif
1472 enddo
1473 do k = 2, km1
1474 do i = 1, im
1475 if (cnvflg(i) .and. k <= kbmax(i)) then
1476 if(k > kbcon1(i) .and. heo(i,k) < hmin(i)) then
1477 lmin(i) = k + 1
1478 hmin(i) = heo(i,k)
1479 endif
1480 endif
1481 enddo
1482 enddo
1483c
1484c make sure that jmin is within the cloud
1485c
1486 do i = 1, im
1487 if(cnvflg(i)) then
1488 jmin(i) = min(lmin(i),ktcon(i)-1)
1489 jmin(i) = max(jmin(i),kbcon1(i)+1)
1490 if(jmin(i) >= ktcon(i)) cnvflg(i) = .false.
1491 endif
1492 enddo
1493c
1494c specify upper limit of mass flux at cloud base
1495c
1497 do i = 1, im
1498 if(cnvflg(i)) then
1499 k = kbcon(i)
1500 dp = 1000. * del(i,k)
1501 xmbmax(i) = dp / (grav * dt2)
1502 endif
1503 enddo
1504c
1505c compute cloud moisture property and precipitation
1506c
1508 do i = 1, im
1509 if (cnvflg(i)) then
1510! aa1(i) = 0.
1511 qcko(i,kb(i)) = qo(i,kb(i))
1512 qrcko(i,kb(i)) = qo(i,kb(i))
1513! rhbar(i) = 0.
1514 endif
1515 enddo
1517 do k = 2, km1
1518 do i = 1, im
1519 if (cnvflg(i)) then
1520 if(k > kb(i) .and. k < ktcon(i)) then
1521 dz = zi(i,k) - zi(i,k-1)
1522 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1523 qrch = qeso(i,k)
1524 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1525cj
1526 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1527 tem = cq * tem
1528 factor = 1. + tem
1529 qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem*
1530 & (qo(i,k)+qo(i,k-1)))/factor
1531 qrcko(i,k) = qcko(i,k)
1532cj
1533 dq = eta(i,k) * (qcko(i,k) - qrch)
1534c
1535! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
1536c
1537c check if there is excess moisture to release latent heat
1538c
1539 if(k >= kbcon(i) .and. dq > 0.) then
1540 etah = .5 * (eta(i,k) + eta(i,k-1))
1541 dp = 1000. * del(i,k)
1542 if(ncloud > 0 .and. k > jmin(i)) then
1543 ptem = c0t(i,k) + c1
1544 qlk = dq / (eta(i,k) + etah * ptem * dz)
1545 dellal(i,k) = etah * c1 * dz * qlk * grav / dp
1546 else
1547 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1548 endif
1549! aa1(i) = aa1(i) - dz * grav * qlk * etah
1550! aa1(i) = aa1(i) - dz * grav * qlk
1551 buo(i,k) = buo(i,k) - grav * qlk
1552 qcko(i,k) = qlk + qrch
1553 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1554 pwavo(i) = pwavo(i) + pwo(i,k)
1555! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp
1556 cnvwt(i,k) = etah * qlk * grav / dp
1557 zdqca(i,k)=dq/eta(i,k)
1558 endif
1559!
1560! compute buoyancy and drag for updraft velocity
1561!
1562 if(k >= kbcon(i)) then
1563 rfact = 1. + fv * cp * gamma
1564 & * to(i,k) / hvap
1565 buo(i,k) = buo(i,k) + (grav / (cp * to(i,k)))
1566 & * dbyo(i,k) / (1. + gamma)
1567 & * rfact
1568 val = 0.
1569 buo(i,k) = buo(i,k) + grav * fv *
1570 & max(val,(qeso(i,k) - qo(i,k)))
1571 drag(i,k) = max(xlamue(i,k),xlamud(i,k))
1572!
1573 tem = ((uo(i,k)-uo(i,k-1))/dz)**2
1574 tem = tem+((vo(i,k)-vo(i,k-1))/dz)**2
1575 wush(i,k) = csmf * sqrt(tem)
1576!
1577 endif
1578!
1579 endif
1580 endif
1581 enddo
1582 enddo
1583c
1584! do i = 1, im
1585! if(cnvflg(i)) then
1586! indx = ktcon(i) - kb(i) - 1
1587! rhbar(i) = rhbar(i) / float(indx)
1588! endif
1589! enddo
1590c
1591c calculate cloud work function
1592c
1593! do k = 2, km1
1594! do i = 1, im
1595! if (cnvflg(i)) then
1596! if(k >= kbcon(i) .and. k < ktcon(i)) then
1597! dz1 = zo(i,k+1) - zo(i,k)
1598! gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1599! rfact = 1. + fv * cp * gamma
1600! & * to(i,k) / hvap
1601! aa1(i) = aa1(i) +
1602!! & dz1 * eta(i,k) * (grav / (cp * to(i,k)))
1603! & dz1 * (grav / (cp * to(i,k)))
1604! & * dbyo(i,k) / (1. + gamma)
1605! & * rfact
1606! val = 0.
1607! aa1(i) = aa1(i) +
1608!! & dz1 * eta(i,k) * grav * fv *
1609! & dz1 * grav * fv *
1610! & max(val,(qeso(i,k) - qo(i,k)))
1611! endif
1612! endif
1613! enddo
1614! enddo
1615!
1616! calculate cloud work function
1617!
1623 do i = 1, im
1624 if (cnvflg(i)) then
1625 aa1(i) = 0.
1626 endif
1627 enddo
1628 do k = 2, km1
1629 do i = 1, im
1630 if (cnvflg(i)) then
1631 if(k >= kbcon(i) .and. k < ktcon(i)) then
1632 dz1 = zo(i,k+1) - zo(i,k)
1633! aa1(i) = aa1(i) + buo(i,k) * dz1 * eta(i,k)
1634 aa1(i) = aa1(i) + buo(i,k) * dz1
1635 dbyo1(i,k) = hcko(i,k) - heso(i,k)
1636 endif
1637 endif
1638 enddo
1639 enddo
1640!
1642 do i = 1, im
1643 if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
1644 enddo
1645!!
1646 totflg = .true.
1647 do i=1,im
1648 totflg = totflg .and. (.not. cnvflg(i))
1649 enddo
1650 if(totflg) return
1651!!
1652c
1653c estimate the onvective overshooting as the level
1654c where the [aafac * cloud work function] becomes zero,
1655c which is the final cloud top
1656c
1658 do i = 1, im
1659 if (cnvflg(i)) then
1660 aa2(i) = aafac * aa1(i)
1661 endif
1662 enddo
1663c
1664 do i = 1, im
1665 flg(i) = cnvflg(i)
1666 ktcon1(i) = ktcon(i)
1667 enddo
1668 do k = 2, km1
1669 do i = 1, im
1670 if (flg(i)) then
1671 if(k >= ktcon(i) .and. k < kmax(i)) then
1672 dz1 = zo(i,k+1) - zo(i,k)
1673 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1674 rfact = 1. + fv * cp * gamma
1675 & * to(i,k) / hvap
1676 aa2(i) = aa2(i) +
1677! & dz1 * eta(i,k) * (grav / (cp * to(i,k)))
1678 & dz1 * (grav / (cp * to(i,k)))
1679 & * dbyo(i,k) / (1. + gamma)
1680 & * rfact
1681! val = 0.
1682! aa2(i) = aa2(i) +
1683!! & dz1 * eta(i,k) * grav * fv *
1684! & dz1 * grav * fv *
1685! & max(val,(qeso(i,k) - qo(i,k)))
1686!NRL MNM: Limit overshooting not to be deeper than half the actual cloud
1687 tem = 0.5 * (zi(i,ktcon(i))-zi(i,kbcon(i)))
1688 tem1 = zi(i,k)-zi(i,ktcon(i))
1689 if(aa2(i) < 0. .or. tem1 >= tem) then
1690 ktcon1(i) = k
1691 flg(i) = .false.
1692 endif
1693 endif
1694 endif
1695 enddo
1696 enddo
1697c
1698c compute cloud moisture property, detraining cloud water
1699c and precipitation in overshooting layers
1700c
1702 do k = 2, km1
1703 do i = 1, im
1704 if (cnvflg(i)) then
1705 if(k >= ktcon(i) .and. k < ktcon1(i)) then
1706 dz = zi(i,k) - zi(i,k-1)
1707 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1708 qrch = qeso(i,k)
1709 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1710cj
1711 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1712 tem = cq * tem
1713 factor = 1. + tem
1714 qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem*
1715 & (qo(i,k)+qo(i,k-1)))/factor
1716 qrcko(i,k) = qcko(i,k)
1717cj
1718 dq = eta(i,k) * (qcko(i,k) - qrch)
1719c
1720c check if there is excess moisture to release latent heat
1721c
1722 if(dq > 0.) then
1723 etah = .5 * (eta(i,k) + eta(i,k-1))
1724 dp = 1000. * del(i,k)
1725 if(ncloud > 0) then
1726 ptem = c0t(i,k) + c1
1727 qlk = dq / (eta(i,k) + etah * ptem * dz)
1728 dellal(i,k) = etah * c1 * dz * qlk * grav / dp
1729 else
1730 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1731 endif
1732 qcko(i,k) = qlk + qrch
1733 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1734 pwavo(i) = pwavo(i) + pwo(i,k)
1735! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp
1736 cnvwt(i,k) = etah * qlk * grav / dp
1737 zdqca(i,k)=dq/eta(i,k)
1738 endif
1739 endif
1740 endif
1741 enddo
1742 enddo
1743!
1744! compute updraft velocity square(wu2)
1746!
1747 if (hwrf_samfdeep) then
1748 do i = 1, im
1749 if (cnvflg(i)) then
1750 k = kbcon1(i)
1751 tem = po(i,k) / (rd * to(i,k))
1752 wucb = -0.01 * dot(i,k) / (tem * grav)
1753 if(wucb.gt.0.) then
1754 wu2(i,k) = wucb * wucb
1755 else
1756 wu2(i,k) = 0.
1757 endif
1758 endif
1759 enddo
1760 endif
1761!
1762 do k = 2, km1
1763 do i = 1, im
1764 if (cnvflg(i)) then
1765 if(k > kbcon1(i) .and. k < ktcon(i)) then
1766 dz = zi(i,k) - zi(i,k-1)
1767 tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
1768 tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
1769 tem2 = wush(i,k) * sqrt(wu2(i,k-1))
1770 tem2 = (tem1 - tem2) * dz
1771 ptem = (1. - tem) * wu2(i,k-1)
1772 ptem1 = 1. + tem
1773 wu2(i,k) = (ptem + tem2) / ptem1
1774 wu2(i,k) = max(wu2(i,k), 0.)
1775 endif
1776 endif
1777 enddo
1778 enddo
1779
1780 if(progsigma)then
1781 do k = 2, km1
1782 do i = 1, im
1783 if (cnvflg(i)) then
1784 if(k > kbcon1(i) .and. k < ktcon(i)) then
1785 rho = po(i,k)*100. / (rd * to(i,k))
1786 omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav
1787 omega_u(i,k)=max(omega_u(i,k),-80.)
1788 endif
1789 endif
1790 enddo
1791 enddo
1792 endif
1793!
1794! compute updraft velocity average over the whole cumulus
1795!
1797 do i = 1, im
1798 wc(i) = 0.
1799 sumx(i) = 0.
1800 enddo
1801 do k = 2, km1
1802 do i = 1, im
1803 if (cnvflg(i)) then
1804 if(k > kbcon1(i) .and. k < ktcon(i)) then
1805 dz = zi(i,k) - zi(i,k-1)
1806 tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1807 wc(i) = wc(i) + tem * dz
1808 sumx(i) = sumx(i) + dz
1809 endif
1810 endif
1811 enddo
1812 enddo
1813 do i = 1, im
1814 if(cnvflg(i)) then
1815 if(sumx(i) == 0.) then
1816 cnvflg(i)=.false.
1817 else
1818 wc(i) = wc(i) / sumx(i)
1819 endif
1820 val = 1.e-4
1821 if (wc(i) < val) cnvflg(i)=.false.
1822 endif
1823 enddo
1824c
1825
1827 if(progsigma)then
1828 do i = 1, im
1829 omegac(i) = 0.
1830 sumx(i) = 0.
1831 enddo
1832 do k = 2, km1
1833 do i = 1, im
1834 if (cnvflg(i)) then
1835 if(k > kbcon1(i) .and. k < ktcon(i)) then
1836 dp = 1000. * del(i,k)
1837 tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
1838 omegac(i) = omegac(i) + tem * dp
1839 sumx(i) = sumx(i) + dp
1840 endif
1841 endif
1842 enddo
1843 enddo
1844 do i = 1, im
1845 if(cnvflg(i)) then
1846 if(sumx(i) == 0.) then
1847 cnvflg(i)=.false.
1848 else
1849 omegac(i) = omegac(i) / sumx(i)
1850 endif
1851 val = -1.2
1852 if (omegac(i) > val) cnvflg(i)=.false.
1853 endif
1854 enddo
1855
1857 do k = 2, km1
1858 do i = 1, im
1859 if (cnvflg(i)) then
1860 if(k >= kbcon1(i) .and. k < ktcon(i)) then
1861 if(omega_u(i,k) .ne. 0.)then
1862 zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k))
1863 else
1864 zeta(i,k)=0.
1865 endif
1866 zeta(i,k)=max(0.,zeta(i,k))
1867 zeta(i,k)=min(1.,zeta(i,k))
1868 endif
1869 endif
1870 enddo
1871 enddo
1872
1873
1874 endif !if progsigma
1875
1876c exchange ktcon with ktcon1
1877c
1879 do i = 1, im
1880 if(cnvflg(i)) then
1881 kk = ktcon(i)
1882 ktcon(i) = ktcon1(i)
1883 ktcon1(i) = kk
1884 endif
1885 enddo
1886c
1887c this section is ready for cloud water
1888c
1890 if(ncloud > 0) then
1891c
1892c compute liquid and vapor separation at cloud top
1893c
1894 do i = 1, im
1895 if(cnvflg(i)) then
1896 k = ktcon(i) - 1
1897 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1898 qrch = qeso(i,k)
1899 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1900 dq = qcko(i,k) - qrch
1901c
1902c check if there is excess moisture to release latent heat
1903c
1904 if(dq > 0.) then
1905 qlko_ktcon(i) = dq
1906 qcko(i,k) = qrch
1907 zdqca(i,k) = dq
1908 endif
1909 endif
1910 enddo
1911 endif
1912c
1913
1914ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then
1915ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
1916ccccc endif
1917c
1918c------- downdraft calculations
1919c
1920c--- compute precipitation efficiency in terms of windshear
1921c
1928 do i = 1, im
1929 if(cnvflg(i)) then
1930 vshear(i) = 0.
1931 endif
1932 enddo
1933 do k = 2, km
1934 do i = 1, im
1935 if (cnvflg(i)) then
1936 if(k > kb(i) .and. k <= ktcon(i)) then
1937 shear = sqrt((uo(i,k)-uo(i,k-1)) ** 2
1938 & + (vo(i,k)-vo(i,k-1)) ** 2)
1939 vshear(i) = vshear(i) + shear
1940 endif
1941 endif
1942 enddo
1943 enddo
1944 do i = 1, im
1945 if(cnvflg(i)) then
1946 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
1947 e1=1.591-.639*vshear(i)
1948 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1949 edt(i)=1.-e1
1950 val = .9
1951 edt(i) = min(edt(i),val)
1952 val = .0
1953 edt(i) = max(edt(i),val)
1954 edto(i)=edt(i)
1955 edtx(i)=edt(i)
1956 endif
1957 enddo
1958c
1959c determine detrainment rate between 1 and kbcon
1960c
1966 do i = 1, im
1967 if(cnvflg(i)) then
1968 sumx(i) = 0.
1969 endif
1970 enddo
1971 do k = 1, km1
1972 do i = 1, im
1973 if(cnvflg(i)) then
1974 if(k >= 1 .and. k < kd94(i)) then
1975 dz = zi(i,k+1) - zi(i,k)
1976 sumx(i) = sumx(i) + dz
1977 endif
1978 endif
1979 enddo
1980 enddo
1981 do i = 1, im
1982 beta = betas
1983 if(islimsk(i) == 1) beta = betal
1984 if(cnvflg(i)) then
1985 dz = (sumx(i)+zi(i,1))/float(kd94(i))
1986 tem = 1./float(kd94(i))
1987 xlamd(i) = (1.-beta**tem)/dz
1988 endif
1989 enddo
1990c
1991c determine downdraft mass flux
1992c
1994 do k = km1, 1, -1
1995 do i = 1, im
1996 if (cnvflg(i) .and. k <= kmax(i)-1) then
1997 if(k < jmin(i) .and. k >= kd94(i)) then
1998 dz = zi(i,k+1) - zi(i,k)
1999 ptem = xlamddt(i) - xlamdet(i)
2000 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
2001 else if(k < kd94(i)) then
2002 dz = zi(i,k+1) - zi(i,k)
2003 ptem = xlamd(i) + xlamddt(i) - xlamdet(i)
2004 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
2005 endif
2006 endif
2007 enddo
2008 enddo
2009c
2010c--- downdraft moisture properties
2011c
2013 do i = 1, im
2014 if(cnvflg(i)) then
2015 jmn = jmin(i)
2016 hcdo(i,jmn) = heo(i,jmn)
2017 qcdo(i,jmn) = qo(i,jmn)
2018 qrcdo(i,jmn)= qo(i,jmn)
2019 ucdo(i,jmn) = uo(i,jmn)
2020 vcdo(i,jmn) = vo(i,jmn)
2021 pwevo(i) = 0.
2022 endif
2023 enddo
2024! for tracers
2025 if (.not.hwrf_samfdeep) then
2026 do n = 1, ntr
2027 do i = 1, im
2028 if(cnvflg(i)) then
2029 jmn = jmin(i)
2030 ecdo(i,jmn,n) = ctro(i,jmn,n)
2031 endif
2032 enddo
2033 enddo
2034 endif
2035cj
2037 do k = km1, 1, -1
2038 do i = 1, im
2039 if (cnvflg(i) .and. k < jmin(i)) then
2040 dz = zi(i,k+1) - zi(i,k)
2041 if(k >= kd94(i)) then
2042 tem = xlamdet(i) * dz
2043 tem1 = 0.5 * xlamddt(i) * dz
2044 else
2045 tem = xlamdet(i) * dz
2046 tem1 = 0.5 * (xlamd(i)+xlamddt(i)) * dz
2047 endif
2048 factor = 1. + tem - tem1
2049 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
2050 & (heo(i,k)+heo(i,k+1)))/factor
2051 dbyo(i,k) = hcdo(i,k) - heso(i,k)
2052!
2053 tem = 0.5 * cm * tem
2054 factor = 1. + tem
2055 ptem = tem - pgcon
2056 ptem1= tem + pgcon
2057 ucdo(i,k) = ((1.-tem)*ucdo(i,k+1)+ptem*uo(i,k+1)
2058 & +ptem1*uo(i,k))/factor
2059 vcdo(i,k) = ((1.-tem)*vcdo(i,k+1)+ptem*vo(i,k+1)
2060 & +ptem1*vo(i,k))/factor
2061 endif
2062 enddo
2063 enddo
2064 if(.not.hwrf_samfdeep) then
2065 do n = 1, ntr
2066 do k = km1, 1, -1
2067 do i = 1, im
2068 if (cnvflg(i) .and. k < jmin(i)) then
2069 dz = zi(i,k+1) - zi(i,k)
2070 tem = 0.5 * xlamdet(i) * dz
2071 tem = cq * tem
2072 factor = 1. + tem
2073 ecdo(i,k,n) = ((1.-tem)*ecdo(i,k+1,n)+tem*
2074 & (ctro(i,k,n)+ctro(i,k+1,n)))/factor
2075 endif
2076 enddo
2077 enddo
2078 enddo
2079 endif
2080c
2082 do k = km1, 1, -1
2083 do i = 1, im
2084 if (cnvflg(i) .and. k < jmin(i)) then
2085 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2086 qrcdo(i,k) = qeso(i,k)+
2087 & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
2088! detad = etad(i,k+1) - etad(i,k)
2089cj
2090 dz = zi(i,k+1) - zi(i,k)
2091 tem = 0.5 * xlamdet(i) * dz
2092 tem = cq * tem
2093 factor = 1. + tem
2094 qcdo(i,k) = ((1.-tem)*qrcdo(i,k+1)+tem*
2095 & (qo(i,k)+qo(i,k+1)))/factor
2096cj
2097! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
2098! & etad(i,k) * qrcdo(i,k)
2099! pwdo(i,k) = pwdo(i,k) - detad *
2100! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
2101cj
2102 pwdo(i,k) = etad(i,k) * (qcdo(i,k) - qrcdo(i,k))
2103 pwevo(i) = pwevo(i) + pwdo(i,k)
2104 endif
2105 enddo
2106 enddo
2107c
2108c--- final downdraft strength dependent on precip
2109c--- efficiency(edt), normalized condensate(pwav), and
2110c--- evaporate(pwev)
2111c
2113 do i = 1, im
2114 edtmax = edtmaxl
2115 if(islimsk(i) == 0) edtmax = edtmaxs
2116 if(cnvflg(i)) then
2117 if(pwevo(i) < 0.) then
2118 edto(i) = -edto(i) * pwavo(i) / pwevo(i)
2119 edto(i) = min(edto(i),edtmax)
2120 else
2121 edto(i) = 0.
2122 endif
2123 endif
2124 enddo
2125c
2126c--- downdraft cloudwork functions
2127c
2129 do k = km1, 1, -1
2130 do i = 1, im
2131 if (cnvflg(i) .and. k < jmin(i)) then
2132 gamma = el2orc * qeso(i,k) / to(i,k)**2
2133 dhh = hcdo(i,k)
2134 dt = to(i,k)
2135 dg = gamma
2136 dh = heso(i,k)
2137 dz = zo(i,k) - zo(i,k+1)
2138! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k)
2139 aa1(i) = aa1(i)+edto(i)*dz
2140 & *(grav/(cp*dt))*((dhh-dh)/(1.+dg))
2141 & *(1.+fv*cp*dg*dt/hvap)
2142 val = 0.
2143! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k)
2144 aa1(i) = aa1(i)+edto(i)*dz
2145 & *grav*fv*max(val,(qeso(i,k)-qo(i,k)))
2146 endif
2147 enddo
2148 enddo
2150 do i = 1, im
2151 if(cnvflg(i) .and. aa1(i) <= 0.) then
2152 cnvflg(i) = .false.
2153 endif
2154 enddo
2155!!
2156 totflg = .true.
2157 do i=1,im
2158 totflg = totflg .and. (.not. cnvflg(i))
2159 enddo
2160 if(totflg) return
2161!!
2162c
2163c--- what would the change be, that a cloud with unit mass
2164c--- will do to the environment?
2165c
2167 do k = 1, km
2168 do i = 1, im
2169 if(cnvflg(i) .and. k <= kmax(i)) then
2170 dellah(i,k) = 0.
2171 dellaq(i,k) = 0.
2172 dellau(i,k) = 0.
2173 dellav(i,k) = 0.
2174 endif
2175 enddo
2176 enddo
2177 if (.not.hwrf_samfdeep) then
2178 do n = 1, ntr
2179 do k = 1, km
2180 do i = 1, im
2181 if(cnvflg(i) .and. k <= kmax(i)) then
2182 dellae(i,k,n) = 0.
2183 endif
2184 enddo
2185 enddo
2186 enddo
2187 endif
2188 do i = 1, im
2189 if(cnvflg(i)) then
2190 dp = 1000. * del(i,1)
2191 tem = edto(i) * etad(i,1) * grav / dp
2192 dellah(i,1) = tem * (hcdo(i,1) - heo(i,1))
2193 dellaq(i,1) = tem * qrcdo(i,1)
2194 dellau(i,1) = tem * (ucdo(i,1) - uo(i,1))
2195 dellav(i,1) = tem * (vcdo(i,1) - vo(i,1))
2196 endif
2197 enddo
2198 if (.not.hwrf_samfdeep) then
2199 do n = 1, ntr
2200 do i = 1, im
2201 if(cnvflg(i)) then
2202 dp = 1000. * del(i,1)
2203 dellae(i,1,n) = edto(i) * etad(i,1) * ecdo(i,1,n)
2204 & * grav / dp
2205 endif
2206 enddo
2207 enddo
2208 endif
2209c
2210c--- changed due to subsidence and entrainment
2211c
2212 do k = 2, km1
2213 do i = 1, im
2214 if (cnvflg(i) .and. k < ktcon(i)) then
2215 aup = 1.
2216 if(k <= kb(i)) aup = 0.
2217 adw = 1.
2218 if(k > jmin(i)) adw = 0.
2219 dp = 1000. * del(i,k)
2220 dz = zi(i,k) - zi(i,k-1)
2221c
2222 dv1h = heo(i,k)
2223 dv2h = .5 * (heo(i,k) + heo(i,k-1))
2224 dv3h = heo(i,k-1)
2225c
2226 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
2227 tem1 = 0.5 * (xlamud(i,k)+xlamud(i,k-1))
2228c
2229 if(k <= kd94(i)) then
2230 ptem = xlamdet(i)
2231 ptem1 = xlamd(i)+xlamddt(i)
2232 else
2233 ptem = xlamdet(i)
2234 ptem1 = xlamddt(i)
2235 endif
2236
2237 factor = grav / dp
2238cj
2239 dellah(i,k) = dellah(i,k) +
2240 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h
2241 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h
2242 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz
2243 & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
2244 & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz
2245 & ) * factor
2246cj
2247 tem1 = -eta(i,k) * qrcko(i,k)
2248 tem2 = -eta(i,k-1) * qcko(i,k-1)
2249 ptem1 = -etad(i,k) * qrcdo(i,k)
2250 ptem2 = -etad(i,k-1) * qcdo(i,k-1)
2251 dellaq(i,k) = dellaq(i,k) +
2252 & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*factor
2253cj
2254 tem1=eta(i,k)*(uo(i,k)-ucko(i,k))
2255 tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1))
2256 ptem1=etad(i,k)*(uo(i,k)-ucdo(i,k))
2257 ptem2=etad(i,k-1)*(uo(i,k-1)-ucdo(i,k-1))
2258 dellau(i,k) = dellau(i,k) +
2259 & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*factor
2260cj
2261 tem1=eta(i,k)*(vo(i,k)-vcko(i,k))
2262 tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1))
2263 ptem1=etad(i,k)*(vo(i,k)-vcdo(i,k))
2264 ptem2=etad(i,k-1)*(vo(i,k-1)-vcdo(i,k-1))
2265 dellav(i,k) = dellav(i,k) +
2266 & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*factor
2267cj
2268 endif
2269 enddo
2270 enddo
2271 if (.not.hwrf_samfdeep) then
2272 do n = 1, ntr
2273 do k = 2, km1
2274 do i = 1, im
2275 if (cnvflg(i) .and. k < ktcon(i)) then
2276 aup = 1.
2277 if(k <= kb(i)) aup = 0.
2278 adw = 1.
2279 if(k > jmin(i)) adw = 0.
2280 dp = 1000. * del(i,k)
2281cj
2282 tem1 = -eta(i,k) * ercko(i,k,n)
2283 tem2 = -eta(i,k-1) * ecko(i,k-1,n)
2284 ptem1 = -etad(i,k) * ecdo(i,k,n)
2285 ptem2 = -etad(i,k-1) * ecdo(i,k-1,n)
2286 dellae(i,k,n) = dellae(i,k,n) +
2287 & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*grav/dp
2288cj
2289 endif
2290 enddo
2291 enddo
2292 enddo
2293 endif
2294c
2295c------- cloud top
2296c
2297 do i = 1, im
2298 if(cnvflg(i)) then
2299 indx = ktcon(i)
2300 dp = 1000. * del(i,indx)
2301 tem = eta(i,indx-1) * grav / dp
2302 dellah(i,indx) = tem * (hcko(i,indx-1) - heo(i,indx-1))
2303 dellaq(i,indx) = tem * qcko(i,indx-1)
2304 dellau(i,indx) = tem * (ucko(i,indx-1) - uo(i,indx-1))
2305 dellav(i,indx) = tem * (vcko(i,indx-1) - vo(i,indx-1))
2306c
2307c cloud water
2308c
2309 dellal(i,indx) = tem * qlko_ktcon(i)
2310 endif
2311 enddo
2312 if (.not.hwrf_samfdeep) then
2313 do n = 1, ntr
2314 do i = 1, im
2315 if(cnvflg(i)) then
2316 indx = ktcon(i)
2317 dp = 1000. * del(i,indx)
2318 dellae(i,indx,n) = eta(i,indx-1) *
2319 & ecko(i,indx-1,n) * grav / dp
2320 endif
2321 enddo
2322 enddo
2323 endif
2324!
2325! compute change rates due to environmental subsidence & uplift
2326! using a positive definite TVD flux-limiter scheme
2327!
2328! for moisture
2329!
2330 do k=1,km1
2331 do i=1,im
2332 if(cnvflg(i) .and. k <= ktcon(i)) then
2333 q_diff(i,k) = q1(i,k) - q1(i,k+1)
2334 endif
2335 enddo
2336 enddo
2337 do i=1,im
2338 if(cnvflg(i)) then
2339 if(q1(i,1) >= 0.) then
2340 q_diff(i,0) = max(0.,2.*q1(i,1)-q1(i,2))-
2341 & q1(i,1)
2342 else
2343 q_diff(i,0) = min(0.,2.*q1(i,1)-q1(i,2))-
2344 & q1(i,1)
2345 endif
2346 endif
2347 enddo
2348!
2349 flxtvd = 0.
2350 do k = 1, km1
2351 do i = 1, im
2352 if(cnvflg(i) .and. k < ktcon(i)) then
2353 tem = 0.
2354 if(k >= kb(i)) tem = eta(i,k)
2355 if(k <= jmin(i)) then
2356 tem = tem - edto(i) * etad(i,k)
2357 endif
2358 if(tem > 0.) then
2359 rrkp = 0.
2360 if(abs(q_diff(i,k)) > 1.e-22)
2361 & rrkp = q_diff(i,k+1) / q_diff(i,k)
2362 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
2363 tem1 = q1(i,k+1) +
2364 & phkp*(qo(i,k)-q1(i,k+1))
2365 flxtvd(i,k) = tem * tem1
2366 elseif(tem < 0.) then
2367 rrkp = 0.
2368 if(abs(q_diff(i,k)) > 1.e-22)
2369 & rrkp = q_diff(i,k-1) / q_diff(i,k)
2370 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
2371 tem1 = q1(i,k) +
2372 & phkp*(qo(i,k)-q1(i,k))
2373 flxtvd(i,k) = tem * tem1
2374 else
2375 tem1 = qo(i,k)
2376 flxtvd(i,k) = 0.
2377 endif
2378!
2379! subtract the double counting change rates at jmin+1 & kb beforehand
2380!
2381 if(k == jmin(i)) then
2382 dp = 1000. * del(i,k+1)
2383 dellaq(i,k+1) = dellaq(i,k+1) -
2384 & edto(i) * etad(i,k) * tem1 * grav/dp
2385 endif
2386 if(k == kb(i)) then
2387 dp = 1000. * del(i,k)
2388 dellaq(i,k) = dellaq(i,k) -
2389 & eta(i,k) * tem1 * grav/dp
2390 endif
2391!
2392 endif
2393 enddo
2394 enddo
2395!
2396 do k=1,km1
2397 do i=1,im
2398 if(cnvflg(i) .and. k <= ktcon(i)) then
2399 dp = 1000. * del(i,k)
2400 dellaq(i,k) = dellaq(i,k) +
2401 & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp
2402 endif
2403 enddo
2404 enddo
2405!
2406! for tracers including TKE & ozone
2407!
2408 if (.not.hwrf_samfdeep) then
2409!
2410 do n=1,ntr
2411 do k=1,km1
2412 do i=1,im
2413 if(cnvflg(i) .and. k <= ktcon(i)) then
2414 e_diff(i,k,n) = ctr(i,k,n) - ctr(i,k+1,n)
2415 endif
2416 enddo
2417 enddo
2418 do i=1,im
2419 if(cnvflg(i)) then
2420 if(ctr(i,1,n) >= 0.) then
2421 e_diff(i,0,n) = max(0.,2.*ctr(i,1,n)-ctr(i,2,n))-
2422 & ctr(i,1,n)
2423 else
2424 e_diff(i,0,n) = min(0.,2.*ctr(i,1,n)-ctr(i,2,n))-
2425 & ctr(i,1,n)
2426 endif
2427 endif
2428 enddo
2429 enddo
2430!
2431 do n=1,ntr
2432!
2433 flxtvd = 0.
2434 do k = 1, km1
2435 do i = 1, im
2436 if(cnvflg(i) .and. k < ktcon(i)) then
2437 tem = 0.
2438 if(k >= kb(i)) tem = eta(i,k)
2439 if(k <= jmin(i)) then
2440 tem = tem - edto(i) * etad(i,k)
2441 endif
2442 if(tem > 0.) then
2443 rrkp = 0.
2444 if(abs(e_diff(i,k,n)) > 1.e-22)
2445 & rrkp = e_diff(i,k+1,n) / e_diff(i,k,n)
2446 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
2447 tem1 = ctr(i,k+1,n) +
2448 & phkp*(ctro(i,k,n)-ctr(i,k+1,n))
2449 flxtvd(i,k) = tem * tem1
2450 elseif(tem < 0.) then
2451 rrkp = 0.
2452 if(abs(e_diff(i,k,n)) > 1.e-22)
2453 & rrkp = e_diff(i,k-1,n) / e_diff(i,k,n)
2454 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
2455 tem1 = ctr(i,k,n) +
2456 & phkp*(ctro(i,k,n)-ctr(i,k,n))
2457 flxtvd(i,k) = tem * tem1
2458 else
2459 tem1 = ctro(i,k,n)
2460 flxtvd(i,k) = 0.
2461 endif
2462!
2463! subtract the double counting change rates at jmin+1 & kb beforehand
2464!
2465 if(k == jmin(i)) then
2466 dp = 1000. * del(i,k+1)
2467 dellae(i,k+1,n) = dellae(i,k+1,n) -
2468 & edto(i)*etad(i,k) * tem1 * grav/dp
2469 endif
2470 if(k == kb(i)) then
2471 dp = 1000. * del(i,k)
2472 dellae(i,k,n) = dellae(i,k,n) -
2473 & eta(i,k) * tem1 * grav/dp
2474 endif
2475!
2476 endif
2477 enddo
2478 enddo
2479!
2480 do k=1,km1
2481 do i=1,im
2482 if(cnvflg(i) .and. k <= ktcon(i)) then
2483 dp = 1000. * del(i,k)
2484 dellae(i,k,n) = dellae(i,k,n) +
2485 & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp
2486 endif
2487 enddo
2488 enddo
2489!
2490 enddo
2491!
2492 endif
2493c
2494c------- final changed variable per unit mass flux
2495c
2497!
2498 if(progsigma)then
2499 dxcrtas=30.e3
2500 dxcrtuf=10.e3
2501 else
2502 dxcrtas=8.e3
2503 dxcrtuf=15.e3
2504 endif
2505
2506
2507 do i = 1, im
2508 asqecflg(i) = cnvflg(i)
2509 if(asqecflg(i) .and. gdx(i) < dxcrtas) then
2510 asqecflg(i) = .false.
2511 endif
2512 enddo
2513
2514!
2516 do k = 1, km
2517 do i = 1, im
2518 if (asqecflg(i) .and. k <= kmax(i)) then
2519 if(k > ktcon(i)) then
2520 qo(i,k) = q1(i,k)
2521 to(i,k) = t1(i,k)
2522 endif
2523 if(k <= ktcon(i)) then
2524 qo(i,k) = dellaq(i,k) * mbdt(i) + q1(i,k)
2525 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
2526 to(i,k) = dellat * mbdt(i) + t1(i,k)
2527 val = 1.e-10
2528 qo(i,k) = max(qo(i,k), val )
2529 endif
2530 endif
2531 enddo
2532 enddo
2533c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2534c
2535c--- the above changed environment is now used to calulate the
2536c--- effect the arbitrary cloud(with unit mass flux)
2537c--- would have on the stability,
2538c--- which then is used to calculate the real mass flux,
2539c--- necessary to keep this change in balance with the large-scale
2540c--- destabilization.
2541c
2542c--- environmental conditions again, first heights
2543c
2547 do k = 1, km
2548 do i = 1, im
2549 if(asqecflg(i) .and. k <= kmax(i)) then
2550 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
2551 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
2552 val = 1.e-8
2553 qeso(i,k) = max(qeso(i,k), val )
2554! tvo(i,k) = to(i,k) + fv * to(i,k) * qo(i,k)
2555 endif
2556 enddo
2557 enddo
2558c
2559c--- moist static energy
2560c
2561!! - Recalculate moist static energy and saturation moist static energy.
2562 do k = 1, km1
2563 do i = 1, im
2564 if(asqecflg(i) .and. k <= kmax(i)-1) then
2565 dz = .5 * (zo(i,k+1) - zo(i,k))
2566 dp = .5 * (pfld(i,k+1) - pfld(i,k))
2567 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
2568 pprime = pfld(i,k+1) + epsm1 * es
2569 qs = eps * es / pprime
2570 dqsdp = - qs / pprime
2571 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
2572 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
2573 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
2574 dt = (grav * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
2575 dq = dqsdt * dt + dqsdp * dp
2576 to(i,k) = to(i,k+1) + dt
2577 qo(i,k) = qo(i,k+1) + dq
2578 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
2579 endif
2580 enddo
2581 enddo
2582 do k = 1, km1
2583 do i = 1, im
2584 if(asqecflg(i) .and. k <= kmax(i)-1) then
2585 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
2586 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
2587 val1 = 1.e-8
2588 qeso(i,k) = max(qeso(i,k), val1)
2589 val2 = 1.e-10
2590 qo(i,k) = max(qo(i,k), val2 )
2591! qo(i,k) = min(qo(i,k),qeso(i,k))
2592 heo(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
2593 & cp * to(i,k) + hvap * qo(i,k)
2594 heso(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
2595 & cp * to(i,k) + hvap * qeso(i,k)
2596 endif
2597 enddo
2598 enddo
2599 do i = 1, im
2600 if(asqecflg(i)) then
2601 k = kmax(i)
2602 heo(i,k) = grav * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
2603 heso(i,k) = grav * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
2604c heo(i,k) = min(heo(i,k),heso(i,k))
2605 endif
2606 enddo
2607c
2608c**************************** static control
2609c
2610c------- moisture and cloud work functions
2611c
2613 do i = 1, im
2614 if(asqecflg(i)) then
2615 xaa0(i) = 0.
2616 xpwav(i) = 0.
2617 endif
2618 enddo
2619c
2620 do i = 1, im
2621 if(asqecflg(i)) then
2622 indx = kb(i)
2623 hcko(i,indx) = heo(i,indx)
2624 qcko(i,indx) = qo(i,indx)
2625 endif
2626 enddo
2627 do k = 2, km1
2628 do i = 1, im
2629 if (asqecflg(i)) then
2630 if(k > kb(i) .and. k <= ktcon(i)) then
2631 dz = zi(i,k) - zi(i,k-1)
2632 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2633 tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
2634 factor = 1. + tem - tem1
2635 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
2636 & (heo(i,k)+heo(i,k-1)))/factor
2637 endif
2638 endif
2639 enddo
2640 enddo
2641 do k = 2, km1
2642 do i = 1, im
2643 if (asqecflg(i)) then
2644 if(k > kb(i) .and. k < ktcon(i)) then
2645 dz = zi(i,k) - zi(i,k-1)
2646 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2647 xdby = hcko(i,k) - heso(i,k)
2648 xqrch = qeso(i,k)
2649 & + gamma * xdby / (hvap * (1. + gamma))
2650cj
2651 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2652 tem = cq * tem
2653 factor = 1. + tem
2654 qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem*
2655 & (qo(i,k)+qo(i,k-1)))/factor
2656cj
2657 dq = eta(i,k) * (qcko(i,k) - xqrch)
2658c
2659 if(k >= kbcon(i) .and. dq > 0.) then
2660 etah = .5 * (eta(i,k) + eta(i,k-1))
2661 if(ncloud > 0 .and. k > jmin(i)) then
2662 ptem = c0t(i,k) + c1
2663 qlk = dq / (eta(i,k) + etah * ptem * dz)
2664 else
2665 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
2666 endif
2667 if(k < ktcon1(i)) then
2668! xaa0(i) = xaa0(i) - dz * grav * qlk * etah
2669 xaa0(i) = xaa0(i) - dz * grav * qlk
2670 endif
2671 qcko(i,k) = qlk + xqrch
2672 xpw = etah * c0t(i,k) * dz * qlk
2673 xpwav(i) = xpwav(i) + xpw
2674 endif
2675 endif
2676 if(k >= kbcon(i) .and. k < ktcon1(i)) then
2677 dz1 = zo(i,k+1) - zo(i,k)
2678 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2679 rfact = 1. + fv * cp * gamma
2680 & * to(i,k) / hvap
2681 xaa0(i) = xaa0(i)
2682! & + dz1 * eta(i,k) * (grav / (cp * to(i,k)))
2683 & + dz1 * (grav / (cp * to(i,k)))
2684 & * xdby / (1. + gamma)
2685 & * rfact
2686 val=0.
2687 xaa0(i) = xaa0(i) +
2688! & dz1 * eta(i,k) * grav * fv *
2689 & dz1 * grav * fv *
2690 & max(val,(qeso(i,k) - qo(i,k)))
2691 endif
2692 endif
2693 enddo
2694 enddo
2695c
2696c------- downdraft calculations
2697c
2698c--- downdraft moisture properties
2699c
2701 do i = 1, im
2702 if(asqecflg(i)) then
2703 jmn = jmin(i)
2704 hcdo(i,jmn) = heo(i,jmn)
2705 qcdo(i,jmn) = qo(i,jmn)
2706 qrcd(i,jmn) = qo(i,jmn)
2707 xpwev(i) = 0.
2708 endif
2709 enddo
2710cj
2711 do k = km1, 1, -1
2712 do i = 1, im
2713 if (asqecflg(i) .and. k < jmin(i)) then
2714 dz = zi(i,k+1) - zi(i,k)
2715 if(k >= kd94(i)) then
2716 tem = xlamdet(i) * dz
2717 tem1 = 0.5 * xlamddt(i) * dz
2718 else
2719 tem = xlamdet(i) * dz
2720 tem1 = 0.5 * (xlamd(i)+xlamddt(i)) * dz
2721 endif
2722 factor = 1. + tem - tem1
2723 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
2724 & (heo(i,k)+heo(i,k+1)))/factor
2725 endif
2726 enddo
2727 enddo
2728cj
2729 do k = km1, 1, -1
2730 do i = 1, im
2731 if (asqecflg(i) .and. k < jmin(i)) then
2732 dq = qeso(i,k)
2733 dt = to(i,k)
2734 gamma = el2orc * dq / dt**2
2735 dh = hcdo(i,k) - heso(i,k)
2736 qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
2737! detad = etad(i,k+1) - etad(i,k)
2738cj
2739 dz = zi(i,k+1) - zi(i,k)
2740 tem = 0.5 * xlamdet(i) * dz
2741 tem = cq * tem
2742 factor = 1. + tem
2743 qcdo(i,k) = ((1.-tem)*qrcd(i,k+1)+tem*
2744 & (qo(i,k)+qo(i,k+1)))/factor
2745cj
2746! xpwd = etad(i,k+1) * qcdo(i,k+1) -
2747! & etad(i,k) * qrcd(i,k)
2748! xpwd = xpwd - detad *
2749! & .5 * (qrcd(i,k) + qrcd(i,k+1))
2750cj
2751 xpwd = etad(i,k) * (qcdo(i,k) - qrcd(i,k))
2752 xpwev(i) = xpwev(i) + xpwd
2753 endif
2754 enddo
2755 enddo
2756c
2757 do i = 1, im
2758 edtmax = edtmaxl
2759 if(islimsk(i) == 0) edtmax = edtmaxs
2760 if(asqecflg(i)) then
2761 if(xpwev(i) >= 0.) then
2762 edtx(i) = 0.
2763 else
2764 edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
2765 edtx(i) = min(edtx(i),edtmax)
2766 endif
2767 endif
2768 enddo
2769c
2770c
2771c--- downdraft cloudwork functions
2772c
2773c
2774 do k = km1, 1, -1
2775 do i = 1, im
2776 if (asqecflg(i) .and. k < jmin(i)) then
2777 gamma = el2orc * qeso(i,k) / to(i,k)**2
2778 dhh = hcdo(i,k)
2779 dt = to(i,k)
2780 dg = gamma
2781 dh = heso(i,k)
2782 dz = zo(i,k) - zo(i,k+1)
2783! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k)
2784 xaa0(i) = xaa0(i)+edtx(i)*dz
2785 & *(grav/(cp*dt))*((dhh-dh)/(1.+dg))
2786 & *(1.+fv*cp*dg*dt/hvap)
2787 val = 0.
2788! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k)
2789 xaa0(i) = xaa0(i)+edtx(i)*dz
2790 & *grav*fv*max(val,(qeso(i,k)-qo(i,k)))
2791 endif
2792 enddo
2793 enddo
2794c
2795c calculate critical cloud work function
2796c
2797! do i = 1, im
2798! if(cnvflg(i)) then
2799! if(pfld(i,ktcon(i)) < pcrit(15))then
2800! acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i)))
2801! & /(975.-pcrit(15))
2802! else if(pfld(i,ktcon(i)) > pcrit(1))then
2803! acrt(i)=acrit(1)
2804! else
2805! k = int((850. - pfld(i,ktcon(i)))/50.) + 2
2806! k = min(k,15)
2807! k = max(k,2)
2808! acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*
2809! & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
2810! endif
2811! endif
2812! enddo
2813! do i = 1, im
2814! if(cnvflg(i)) then
2815! if(islimsk(i) == 1) then
2816! w1 = w1l
2817! w2 = w2l
2818! w3 = w3l
2819! w4 = w4l
2820! else
2821! w1 = w1s
2822! w2 = w2s
2823! w3 = w3s
2824! w4 = w4s
2825! endif
2826c
2827c modify critical cloud workfunction by cloud base vertical velocity
2828c
2829! if(pdot(i) <= w4) then
2830! acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
2831! elseif(pdot(i) >= -w4) then
2832! acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
2833! else
2834! acrtfct(i) = 0.
2835! endif
2836! val1 = -1.
2837! acrtfct(i) = max(acrtfct(i),val1)
2838! val2 = 1.
2839! acrtfct(i) = min(acrtfct(i),val2)
2840! acrtfct(i) = 1. - acrtfct(i)
2841c
2842c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
2843c
2844c if(rhbar(i) >= .8) then
2845c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
2846c endif
2847c
2848c modify adjustment time scale by cloud base vertical velocity
2849c
2850! dtconv(i) = dt2 + max((1800. - dt2),0.) *
2851! & (pdot(i) - w2) / (w1 - w2)
2852c dtconv(i) = max(dtconv(i), dt2)
2853c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
2854!
2855! dtconv(i) = max(dtconv(i),dtmin)
2856! dtconv(i) = min(dtconv(i),dtmax)
2857c
2858! endif
2859! enddo
2860!
2861! compute convective turn-over time
2862!
2864 if(hwrf_samfdeep) then
2865 do i= 1, im
2866 if(cnvflg(i)) then
2867 tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
2868 dtconv(i) = tem / wc(i)
2869 dtconv(i) = max(dtconv(i),dtmin)
2870 dtconv(i) = min(dtconv(i),dtmax)
2871 endif
2872 enddo
2873 else
2874 do i= 1, im
2875 if(cnvflg(i)) then
2876 tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
2877 dtconv(i) = tem / wc(i)
2878 tfac = 1. + gdx(i) / 75000.
2879 dtconv(i) = tfac * dtconv(i)
2880 dtconv(i) = max(dtconv(i),dtmin)
2881 dtconv(i) = min(dtconv(i),dtmax)
2882 endif
2883 enddo
2884 endif
2885!
2887 do i= 1, im
2888 if(cnvflg(i)) then
2889 sumx(i) = 0.
2890 umean(i) = 0.
2891 endif
2892 enddo
2893 do k = 2, km1
2894 do i = 1, im
2895 if(cnvflg(i)) then
2896 if(k >= kbcon1(i) .and. k < ktcon1(i)) then
2897 dz = zi(i,k) - zi(i,k-1)
2898 tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
2899 umean(i) = umean(i) + tem * dz
2900 sumx(i) = sumx(i) + dz
2901 endif
2902 endif
2903 enddo
2904 enddo
2905 do i= 1, im
2906 if(cnvflg(i)) then
2907 umean(i) = umean(i) / sumx(i)
2908 umean(i) = max(umean(i), 1.)
2909 tauadv = gdx(i) / umean(i)
2910 advfac(i) = tauadv / dtconv(i)
2911 advfac(i) = min(advfac(i), 1.)
2912 endif
2913 enddo
2914
2916 if(progsigma)then
2917
2918!Initial computations, dynamic q-tendency
2919 if(first_time_step .and. .not.restart)then
2920 do k = 1,km
2921 do i = 1,im
2922 qadv(i,k)=0.
2923 enddo
2924 enddo
2925 else
2926 do k = 1,km
2927 do i = 1,im
2928 qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
2929 enddo
2930 enddo
2931 endif
2932
2933 do k = 1,km
2934 do i = 1,im
2935 tmfq(i,k)=tmf(i,k,1)
2936 enddo
2937 enddo
2938
2939 flag_shallow = .false.
2940 flag_mid = .false.
2941 call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
2942 & flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,
2943 & delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
2944 & sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
2945 endif
2946
2949
2950 do i= 1, im
2951 if(cnvflg(i) .and. .not.asqecflg(i)) then
2952 k = kbcon(i)
2953 rho = po(i,k)*100. / (rd*to(i,k))
2954 if(progsigma)then
2955 xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv)
2956 else
2957 xmb(i) = advfac(i)*betaw*rho*wc(i)
2958 endif
2959 endif
2960 enddo
2966 do i= 1, im
2967 if(asqecflg(i)) then
2968! fld(i)=(aa1(i)-acrt(i)*acrtfct(i))/dtconv(i)
2969 fld(i)=aa1(i)/dtconv(i)
2970 if(fld(i) <= 0.) then
2971 asqecflg(i) = .false.
2972 cnvflg(i) = .false.
2973 endif
2974 endif
2980 if(asqecflg(i)) then
2981c xaa0(i) = max(xaa0(i),0.)
2982 xk(i) = (xaa0(i) - aa1(i)) / mbdt(i)
2983 if(xk(i) >= 0.) then
2984 asqecflg(i) = .false.
2985 cnvflg(i) = .false.
2986 endif
2987 endif
2988c
2989c--- kernel, cloud base mass flux
2990c
2997 if(asqecflg(i)) then
2998 xmb(i) = -advfac(i) * fld(i) / xk(i)
2999 endif
3000 enddo
3001!!
3002
3004 totflg = .true.
3005 do i=1,im
3006 totflg = totflg .and. (.not. cnvflg(i))
3007 enddo
3008 if(totflg) return
3009!!
3010!
3012 do i = 1, im
3013 if(cnvflg(i)) then
3014 tem = min(max(xlamue(i,kbcon(i)), 7.e-5), 3.e-4)
3015 tem = 0.2 / tem
3016 tem1 = 3.14 * tem * tem
3017 sigmagfm(i) = tem1 / garea(i)
3018 sigmagfm(i) = max(sigmagfm(i), 0.001)
3019 sigmagfm(i) = min(sigmagfm(i), 0.999)
3020 endif
3021 enddo
3022!
3024 do i = 1, im
3025 if(cnvflg(i)) then
3026 if (gdx(i) < dxcrtuf) then
3027 if(progsigma)then
3028 scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i))
3029 else
3030 scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
3031 endif
3032 scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
3033 else
3034 scaldfunc(i) = 1.0
3035 endif
3036 xmb(i) = xmb(i) * scaldfunc(i)
3037 xmb(i) = min(xmb(i),xmbmax(i))
3038 endif
3039 enddo
3040!
3042!
3043! if (do_aerosols)
3044! & call samfdeepcnv_aerosols(im, im, km, itc, ntc, ntr, delt,
3045! & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kd94, ktcon, fscav,
3046! & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp,
3047! & qtr, qaero)
3048!
3049c
3050c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
3051c
3052 do k = 1, km
3053 do i = 1, im
3054 if (cnvflg(i) .and. k <= kmax(i)) then
3055 to(i,k) = t1(i,k)
3056 qo(i,k) = q1(i,k)
3057 uo(i,k) = u1(i,k)
3058 vo(i,k) = v1(i,k)
3059 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
3060 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
3061 val = 1.e-8
3062 qeso(i,k) = max(qeso(i,k), val )
3063 endif
3064 enddo
3065 enddo
3066 if (.not.hwrf_samfdeep) then
3067 do n = 1, ntr
3068 kk = n+2
3069 do k = 1, km
3070 do i = 1, im
3071 if (cnvflg(i) .and. k <= kmax(i)) then
3072 ctro(i,k,n) = qtr(i,k,kk)
3073 endif
3074 enddo
3075 enddo
3076 enddo
3077 endif
3078c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3079c
3080c--- feedback: simply the changes from the cloud with unit mass flux
3081c--- multiplied by the mass flux necessary to keep the
3082c--- equilibrium with the larger-scale.
3083c
3088 do i = 1, im
3089 delhbar(i) = 0.
3090 delqbar(i) = 0.
3091 deltbar(i) = 0.
3092 delubar(i) = 0.
3093 delvbar(i) = 0.
3094 qcond(i) = 0.
3095 enddo
3096 if (.not.hwrf_samfdeep) then
3097 do n = 1, ntr
3098 do i = 1, im
3099 delebar(i,n) = 0.
3100 enddo
3101 enddo
3102 endif
3103 do k = 1, km
3104 do i = 1, im
3105 if (cnvflg(i) .and. k <= kmax(i)) then
3106 if(k <= ktcon(i)) then
3107 tem2 = xmb(i) * dt2
3108 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
3109 t1(i,k) = t1(i,k) + tem2 * dellat
3110 q1(i,k) = q1(i,k) + tem2 * dellaq(i,k)
3111! tem = tem2 / rcs(i)
3112! u1(i,k) = u1(i,k) + dellau(i,k) * tem
3113! v1(i,k) = v1(i,k) + dellav(i,k) * tem
3114 u1(i,k) = u1(i,k) + tem2 * dellau(i,k)
3115 v1(i,k) = v1(i,k) + tem2 * dellav(i,k)
3116 dp = 1000. * del(i,k)
3117 tem = xmb(i) * dp / grav
3118 delhbar(i) = delhbar(i) + tem * dellah(i,k)
3119 delqbar(i) = delqbar(i) + tem * dellaq(i,k)
3120 deltbar(i) = deltbar(i) + tem * dellat
3121 delubar(i) = delubar(i) + tem * dellau(i,k)
3122 delvbar(i) = delvbar(i) + tem * dellav(i,k)
3123 endif
3124 endif
3125 enddo
3126 enddo
3127!
3128! Negative moisture is set to zero after borrowing it from
3129! positive values within the mass-flux transport layers
3130!
3131 do i = 1,im
3132 tsumn(i) = 0.
3133 tsump(i) = 0.
3134 rtnp(i) = 1.
3135 enddo
3136 do k = 1,km1
3137 do i = 1,im
3138 if(cnvflg(i) .and. k <= ktcon(i)) then
3139 tem = q1(i,k) * delp(i,k) / grav
3140 if(q1(i,k) < 0.) tsumn(i) = tsumn(i) + tem
3141 if(q1(i,k) > 0.) tsump(i) = tsump(i) + tem
3142 endif
3143 enddo
3144 enddo
3145 do i = 1,im
3146 if(cnvflg(i)) then
3147 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
3148 if(tsump(i) > abs(tsumn(i))) then
3149 rtnp(i) = tsumn(i) / tsump(i)
3150 else
3151 rtnp(i) = tsump(i) / tsumn(i)
3152 endif
3153 endif
3154 endif
3155 enddo
3156 do k = 1,km1
3157 do i = 1,im
3158 if(cnvflg(i) .and. k <= ktcon(i)) then
3159 if(rtnp(i) < 0.) then
3160 if(tsump(i) > abs(tsumn(i))) then
3161 if(q1(i,k) < 0.) q1(i,k) = 0.
3162 if(q1(i,k) > 0.) q1(i,k) = (1.+rtnp(i))*q1(i,k)
3163 else
3164 if(q1(i,k) < 0.) q1(i,k) = (1.+rtnp(i))*q1(i,k)
3165 if(q1(i,k) > 0.) q1(i,k) = 0.
3166 endif
3167 endif
3168 endif
3169 enddo
3170 enddo
3171!
3172 if (.not.hwrf_samfdeep) then
3173 indx = ntk - 2
3174 do n = 1, ntr
3175!
3176 do k = 1, km
3177 do i = 1, im
3178 if (cnvflg(i) .and. k <= kmax(i)) then
3179 if(k <= ktcon(i)) then
3180 ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2
3181 dp = 1000. * del(i,k)
3182 delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/grav
3183 endif
3184 endif
3185 enddo
3186 enddo
3187!
3188! Negative TKE, ozone, and aerosols are set to zero after borrowing them
3189! from positive values within the mass-flux transport layers
3190!
3191 do i = 1,im
3192 tsumn(i) = 0.
3193 tsump(i) = 0.
3194 rtnp(i) = 1.
3195 enddo
3196 do k = 1,km1
3197 do i = 1,im
3198 if(cnvflg(i) .and. k <= ktcon(i)) then
3199 if(n == indx) then
3200 if(k > 1) then
3201 dz = zi(i,k) - zi(i,k-1)
3202 else
3203 dz = zi(i,k)
3204 endif
3205 tem = ctr(i,k,n) * dz
3206 else
3207 tem = ctr(i,k,n) * delp(i,k) / grav
3208 endif
3209 if(ctr(i,k,n) < 0.) tsumn(i) = tsumn(i) + tem
3210 if(ctr(i,k,n) > 0.) tsump(i) = tsump(i) + tem
3211 endif
3212 enddo
3213 enddo
3214 do i = 1,im
3215 if(cnvflg(i)) then
3216 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
3217 if(tsump(i) > abs(tsumn(i))) then
3218 rtnp(i) = tsumn(i) / tsump(i)
3219 else
3220 rtnp(i) = tsump(i) / tsumn(i)
3221 endif
3222 endif
3223 endif
3224 enddo
3225 do k = 1,km1
3226 do i = 1,im
3227 if(cnvflg(i) .and. k <= ktcon(i)) then
3228 if(rtnp(i) < 0.) then
3229 if(tsump(i) > abs(tsumn(i))) then
3230 if(ctr(i,k,n)<0.) ctr(i,k,n)=0.
3231 if(ctr(i,k,n)>0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n)
3232 else
3233 if(ctr(i,k,n)<0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n)
3234 if(ctr(i,k,n)>0.) ctr(i,k,n)=0.
3235 endif
3236 endif
3237 endif
3238 enddo
3239 enddo
3240!
3241 kk = n+2
3242 do k = 1, km
3243 do i = 1, im
3244 if(cnvflg(i) .and. k <= ktcon(i)) then
3245 qtr(i,k,kk) = ctr(i,k,n)
3246 endif
3247 enddo
3248 enddo
3249!
3250 enddo
3251!
3252 if (do_aerosols) then
3253!
3254 do n = 1, ntc
3255!
3256! convert wet deposition to total mass deposited over dt2 and dp
3257 do k = 2, km1
3258 do i = 1, im
3259 if (cnvflg(i)) then
3260 if(k > kb(i) .and. k < ktcon(i)) then
3261 dp = 1000. * del(i,k)
3262 wet_dep(i,k,n) = chem_pw(i,k,n)*grav/dp
3263 wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*dt2*dp
3264 endif
3265 endif
3266 enddo
3267 enddo
3268!
3269 kk = n + itc - 1
3270 do k = 2, km1
3271 do i = 1, im
3272 if (cnvflg(i)) then
3273 if(k > kb(i) .and. k < ktcon(i)) then
3274 dp = 1000. * del(i,k)
3275 if (qtr(i,k,kk) < 0.) then
3276! borrow negative mass from wet deposition
3277 tem = -qtr(i,k,kk)*dp
3278 if(wet_dep(i,k,n) >= tem) then
3279 wet_dep(i,k,n) = wet_dep(i,k,n) - tem
3280 qtr(i,k,kk) = 0.
3281 else
3282 wet_dep(i,k,n) = 0.
3283 qtr(i,k,kk) = qtr(i,k,kk)+wet_dep(i,k,n)/dp
3284 endif
3285 endif
3286 endif
3287 endif
3288 enddo
3289 enddo
3290!
3291 enddo
3292!
3293 endif
3294!
3295 endif
3296!
3298 do k = 1, km
3299 do i = 1, im
3300 if (cnvflg(i) .and. k <= kmax(i)) then
3301 if(k <= ktcon(i)) then
3302 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
3303 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
3304 val = 1.e-8
3305 qeso(i,k) = max(qeso(i,k), val )
3306 endif
3307 endif
3308 enddo
3309 enddo
3310c
3312 do i = 1, im
3313 rntot(i) = 0.
3314 delqev(i) = 0.
3315 delq2(i) = 0.
3316 flg(i) = cnvflg(i)
3317 enddo
3318 do k = km, 1, -1
3319 do i = 1, im
3320 if (cnvflg(i) .and. k <= kmax(i)) then
3321 if(k < ktcon(i)) then
3322 aup = 1.
3323 if(k <= kb(i)) aup = 0.
3324 adw = 1.
3325 if(k >= jmin(i)) adw = 0.
3326 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
3327 rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
3328 endif
3329 endif
3330 enddo
3331 enddo
3335 do k = km, 1, -1
3336 do i = 1, im
3337 if (k <= kmax(i)) then
3338 deltv(i) = 0.
3339 delq(i) = 0.
3340 qevap(i) = 0.
3341 if(cnvflg(i) .and. k < ktcon(i)) then
3342 aup = 1.
3343 if(k <= kb(i)) aup = 0.
3344 adw = 1.
3345 if(k >= jmin(i)) adw = 0.
3346 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
3347 rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
3348 endif
3349 if(flg(i) .and. k < ktcon(i)) then
3350! evef = edt(i) * evfact
3351! if(islimsk(i) == 1) evef=edt(i) * evfactl
3352! if(islimsk(i) == 1) evef=.07
3353 qcond(i) = evef * (q1(i,k) - qeso(i,k))
3354 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
3355 dp = 1000. * del(i,k)
3356 tem = grav / dp
3357 tem1 = dp / grav
3358 if(rn(i) > 0. .and. qcond(i) < 0.) then
3359 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
3360 qevap(i) = min(qevap(i), rn(i)*1000.*tem)
3361 delq2(i) = delqev(i) + .001 * qevap(i) * tem1
3362 endif
3363 if(rn(i) > 0. .and. qcond(i) < 0. .and.
3364 & delq2(i) > rntot(i)) then
3365 qevap(i) = 1000.* tem * (rntot(i) - delqev(i))
3366 flg(i) = .false.
3367 endif
3368 if(rn(i) > 0. .and. qevap(i) > 0.) then
3369 q1(i,k) = q1(i,k) + qevap(i)
3370 t1(i,k) = t1(i,k) - elocp * qevap(i)
3371 rn(i) = rn(i) - .001 * qevap(i) * tem1
3372 deltv(i) = - elocp*qevap(i)/dt2
3373 delq(i) = + qevap(i)/dt2
3374 delqev(i) = delqev(i) + .001 * qevap(i) * tem1
3375 endif
3376 delqbar(i) = delqbar(i) + delq(i) * tem1
3377 deltbar(i) = deltbar(i) + deltv(i) * tem1
3378 endif
3379 endif
3380 enddo
3381 enddo
3382
3383!LB:
3384 if(do_ca)then
3385 do i = 1,im
3386 rainevap(i)=delqev(i)
3387 enddo
3388 endif
3389cj
3390! do i = 1, 4
3391! if(me == 31 .and. cnvflg(i)) then
3392! if(cnvflg(i)) then
3393! if(i==1) print*,'ntr:ntk= ',ntr,ntk
3394! print *, ' deep delhbar, delqbar, deltbar = ',
3395! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
3396! print *, ' deep delebar ozone = ',delebar(i,1)
3397! print *, ' deep delebar tke = ',delebar(i,2)
3398! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
3399! print *, ' precip =', hvap*rn(i)*1000./dt2
3400! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
3401! endif
3402! enddo
3403c
3404c precipitation rate converted to actual precip
3405c in unit of m instead of kg
3406c
3407 do i = 1, im
3408 if(cnvflg(i)) then
3409c
3410c in the event of upper level rain evaporation and lower level downdraft
3411c moistening, rn can become negative, in this case, we back out of the
3412c heating and the moistening
3413c
3414 if(rn(i) < 0. .and. .not.flg(i)) rn(i) = 0.
3415 if(rn(i) <= 0.) then
3416 rn(i) = 0.
3417 else
3418 ktop(i) = ktcon(i)
3419 kbot(i) = kbcon(i)
3420 kcnv(i) = 1
3421 cldwrk(i) = aa1(i)
3422 endif
3423 endif
3424 enddo
3425c
3426c convective cloud water
3427c
3429 do k = 1, km
3430 do i = 1, im
3431 if (cnvflg(i) .and. rn(i) > 0.) then
3432 if (k >= kbcon(i) .and. k < ktcon(i)) then
3433 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
3434 endif
3435 endif
3436 enddo
3437 enddo
3438c
3439c convective cloud cover
3440c
3442 do k = 1, km
3443 do i = 1, im
3444 if (cnvflg(i) .and. rn(i) > 0.) then
3445 if (k >= kbcon(i) .and. k < ktcon(i)) then
3446 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
3447 cnvc(i,k) = min(cnvc(i,k), 0.6)
3448 cnvc(i,k) = max(cnvc(i,k), 0.0)
3449 endif
3450 endif
3451 enddo
3452 enddo
3453c
3454c cloud water
3455c
3457 if (ncloud > 0) then
3458!
3459 do k = 1, km
3460 do i = 1, im
3461 if (cnvflg(i) .and. rn(i) > 0.) then
3462! if (k > kb(i) .and. k <= ktcon(i)) then
3463 if (k >= kbcon(i) .and. k <= ktcon(i)) then
3464 tem = dellal(i,k) * xmb(i) * dt2
3465 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
3466 if (qtr(i,k,2) > -999.0) then
3467 qtr(i,k,1) = qtr(i,k,1) + tem * tem1 ! ice
3468 qtr(i,k,2) = qtr(i,k,2) + tem *(1.0-tem1) ! water
3469 else
3470 qtr(i,k,1) = qtr(i,k,1) + tem
3471 endif
3472 endif
3473 endif
3474 enddo
3475 enddo
3476!
3477 endif
3478c
3480 do k = 1, km
3481 do i = 1, im
3482 if(cnvflg(i) .and. rn(i) <= 0.) then
3483 if (k <= kmax(i)) then
3484 t1(i,k) = to(i,k)
3485 q1(i,k) = qo(i,k)
3486 u1(i,k) = uo(i,k)
3487 v1(i,k) = vo(i,k)
3488 endif
3489 endif
3490 enddo
3491 enddo
3492 if (.not.hwrf_samfdeep) then
3493 do n = 1, ntr
3494 kk = n+2
3495 do k = 1, km
3496 do i = 1, im
3497 if(cnvflg(i) .and. rn(i) <= 0.) then
3498 if (k <= kmax(i)) then
3499 qtr(i,k,kk)= ctro(i,k,n)
3500 endif
3501 endif
3502 enddo
3503 enddo
3504 enddo
3505 if (do_aerosols) then
3506 do n = 1, ntc
3507 do k = 2, km1
3508 do i = 1, im
3509 if(cnvflg(i) .and. rn(i) <= 0.) then
3510 if (k <= ktcon(i)) then
3511 wet_dep(i,k,n) = 0.
3512 endif
3513 endif
3514 enddo
3515 enddo
3516 enddo
3517 endif
3519! if (do_aerosols) then
3520! do n = 1, ntc
3521! kk = n + itc - 1
3522! do k = 1, km
3523! do i = 1, im
3524! if(cnvflg(i) .and. rn(i) > 0.) then
3525! if (k <= kmax(i)) qtr(i,k,kk) = qaero(i,k,n)
3526! endif
3527! enddo
3528! enddo
3529! enddo
3530! endif
3531 endif
3532!
3533! hchuang code change
3534!
3536!
3538 do k = 1, km
3539 do i = 1, im
3540 if(cnvflg(i) .and. rn(i) > 0.) then
3541 if(k >= kb(i) .and. k < ktop(i)) then
3542 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
3543 endif
3544 endif
3545 enddo
3546 enddo
3548 do i = 1, im
3549 if(cnvflg(i) .and. rn(i) > 0.) then
3550 k = ktop(i)-1
3551 dt_mf(i,k) = ud_mf(i,k)
3552 endif
3553 enddo
3555 do k = 1, km
3556 do i = 1, im
3557 if(cnvflg(i) .and. rn(i) > 0.) then
3558 if(k >= 1 .and. k <= jmin(i)) then
3559 dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
3560 endif
3561 endif
3562 enddo
3563 enddo
3564!
3565! include TKE contribution from deep convection
3566!
3567 if (.not.hwrf_samfdeep) then
3568 if (ntk > 0) then
3569!
3570 do k = 2, km1
3571 do i = 1, im
3572 if(cnvflg(i) .and. rn(i) > 0.) then
3573 if(k > kb(i) .and. k < ktop(i)) then
3574 tem = 0.5 * (eta(i,k-1) + eta(i,k)) * xmb(i)
3575 tem1 = pfld(i,k) * 100. / (rd * t1(i,k))
3576 if(progsigma)then
3577 tem2 = sigmab(i)
3578 else
3579 tem2 = max(sigmagfm(i), betaw)
3580 endif
3581 ptem = tem / (tem2 * tem1)
3582 qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem
3583 endif
3584 endif
3585 enddo
3586 enddo
3587!
3588 do k = 2, km1
3589 do i = 1, im
3590 if(cnvflg(i) .and. rn(i) > 0.) then
3591 if(k > 1 .and. k <= jmin(i)) then
3592 tem = 0.5*edto(i)*(etad(i,k-1)+etad(i,k))*xmb(i)
3593 tem1 = pfld(i,k) * 100. / (rd * t1(i,k))
3594 if(progsigma)then
3595 tem2 = sigmab(i)
3596 else
3597 tem2 = max(sigmagfm(i), betaw)
3598 endif
3599 ptem = tem / (tem2 * tem1)
3600 qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem
3601 endif
3602 endif
3603 enddo
3604 enddo
3605!
3606 endif
3607!!
3608 if(mp_phys == mp_phys_mg) then
3609 do k=1,km
3610 do i=1,im
3611 qlcn(i,k) = qtr(i,k,2) - qlcn(i,k)
3612 qicn(i,k) = qtr(i,k,1) - qicn(i,k)
3613 cf_upi(i,k) = cnvc(i,k)
3614 w_upi(i,k) = ud_mf(i,k)*t1(i,k)*rd /
3615 & (dt2*max(sigmagfm(i),1.e-12)*prslp(i,k))
3616 cnv_mfd(i,k) = ud_mf(i,k)/dt2
3617 clcn(i,k) = cnvc(i,k)
3618 cnv_fice(i,k) = qicn(i,k)
3619 & / max(1.e-10,qlcn(i,k)+qicn(i,k))
3620 enddo
3621 enddo
3622 endif
3623 endif ! (.not.hwrf_samfdeep)
3624 return
3625 end subroutine samfdeepcnv_run
3626
3628
3629 end module samfdeepcnv
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
samfdeepcnv::samfdeepcnv_run
subroutine samfdeepcnv_run(im, km, first_time_step, restart, tmf, qmicro, itc, ntc, cliq, cp, cvap, eps, epsm1, fv, grav, hvap, rd, rv, t0c, delt, ntk, ntr, delp, prslp, psp, phil, qtr, prevsq, q, q1, t1, u1, v1, fscav, hwrf_samfdeep, progsigma, cldwrk, rn, kbot, ktop, kcnv, islimsk, garea, dot, ncloud, hpbl, ud_mf, dd_mf, dt_mf, cnvw, cnvc, qlcn, qicn, w_upi, cf_upi, cnv_mfd, cnv_dqldt, clcn, cnv_fice, cnv_ndrop, cnv_nice, mp_phys, mp_phys_mg, clam, c0s, c1, betal, betas, evef, pgcon, asolfac, do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, rainevap, sigmain, sigmaout, betadcu, betamcu, betascu, maxmf, do_mynnedmf, errmsg, errflg)
Definition samfdeepcnv.f:88
progsigma
This module contains the subroutine that calculates the prognostic updraft area fraction that is used...
Definition progsigma_calc.f90:7
samfdeepcnv
This module contains the CCPP-compliant scale-aware mass-flux deep convection scheme.
Definition samfdeepcnv.f:7