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