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