CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
sascnvn.F
Go to the documentation of this file.
1
11
14 module sascnvn
15
16 implicit none
17
18 private
19
20 public :: sascnvn_init, sascnvn_run
21
22 contains
23
24!!
25!! \section arg_table_sascnvn_init Argument Table
26!! \htmlinclude sascnvn_init.html
27!!
28 subroutine sascnvn_init(imfdeepcnv,imfdeepcnv_sas,errmsg,errflg)
29!
30 integer, intent(in) :: imfdeepcnv, imfdeepcnv_sas
31 character(len=*), intent(out) :: errmsg
32 integer, intent(out) :: errflg
33!
34! initialize CCPP error handling variables
35 errmsg = ''
36 errflg = 0
37!
38 if (imfdeepcnv/=imfdeepcnv_sas) then
39 write(errmsg,'(*(a))') 'Logic error: sascnvn incompatible with',&
40 & ' value of imfdeepcnv'
41 errflg = 1
42 return
43 endif
44!
45 end subroutine sascnvn_init
46
92 subroutine sascnvn_run(
93 & grav,cp,hvap,rv,fv,t0c,rgas,cvap,cliq,eps,epsm1, &
94 & im,km,jcap,delt,delp,prslp,psp,phil,qlc,qli, &
95 & q1,t1,u1,v1,cldwrk,rn,kbot,ktop,kcnv,islimsk, &
96 & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc, &
97 & qlcn,qicn,w_upi,cf_upi,cnv_mfd, &
98 & cnv_dqldt,clcn,cnv_fice,cnv_ndrop,cnv_nice,mp_phys, &
99 & mp_phys_mg,clam,c0,c1,betal,betas,evfact,evfactl,pgcon, &
100 & errmsg,errflg)
101!
102 use machine , only : kind_phys
103 use funcphys , only : fpvs
104! use physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
105! &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
106! &, cvap => con_cvap, cliq => con_cliq &
107! &, eps => con_eps, epsm1 => con_epsm1,rgas => con_rd
108 implicit none
109!
110! Interface variables
111!
112 real(kind=kind_phys), intent(in) :: grav, cp, hvap, rv, fv, t0c, &
113 & rgas, cvap, cliq, eps, epsm1
114 integer, intent(in) :: im, km, jcap, ncloud, &
115 & mp_phys, mp_phys_mg
116 integer, intent(inout) :: kbot(:), ktop(:), kcnv(:)
117 integer, intent(in) :: islimsk(:)
118 real(kind=kind_phys), intent(in) :: delt, clam, c0, c1, pgcon
119 real(kind=kind_phys), intent(in) :: betal, betas, evfact, evfactl
120 real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), &
121 & prslp(:,:), dot(:,:), &
122 & phil(:,:)
123 real(kind=kind_phys), intent(inout) :: &
124 & qlc(:,:), qli(:,:), &
125 & q1(:,:), t1(:,:), &
126 & u1(:,:), v1(:,:), &
127 & cnvw(:,:), cnvc(:,:)
128 real(kind=kind_phys), intent(out) :: cldwrk(:), rn(:), &
129 & dt_mf(:,:), dd_mf(:,:)
130 real(kind=kind_phys), intent(out), optional :: ud_mf(:,:)
131 real(kind=kind_phys), intent(inout), optional :: &
132 & qlcn(:,:), qicn(:,:), &
133 & w_upi(:,:), cnv_mfd(:,:), &
134 & cnv_dqldt(:,:), clcn(:,:), &
135 & cnv_fice(:,:), cnv_ndrop(:,:),&
136 & cnv_nice(:,:), cf_upi(:,:)
137 character(len=*), intent(out) :: errmsg
138 integer, intent(out) :: errflg
139!
140! Local variables
141!
142 integer i, indx, jmn, k, kk, km1
143! integer latd,lond
144!
145 real(kind=kind_phys) cxlamu, xlamde, xlamdd
146!
147! real(kind=kind_phys) detad
148 real(kind=kind_phys) adw, aup, aafac,
149 & beta,
150 & dellat, delta,
151 & desdt, dg,
152 & dh, dhh, dp,
153 & dq, dqsdp, dqsdt, dt,
154 & dt2, dtmax, dtmin, dv1h,
155 & dv1q, dv2h, dv2q, dv1u,
156 & dv1v, dv2u, dv2v, dv3q,
157 & dv3h, dv3u, dv3v,
158 & dz, dz1, e1, edtmax,
159 & edtmaxl, edtmaxs, el2orc, elocp,
160 & es, etah, cthk, dthk,
161 & evef, fact1,
162 & fact2, factor, fjcap, fkm,
163 & g, gamma, pprime,
164 & qlk, qrch, qs,
165 & rain, rfact, shear, tem1,
166 & val, val1,
167 & val2, w1, w1l, w1s,
168 & w2, w2l, w2s, w3,
169 & w3l, w3s, w4, w4l,
170 & w4s, xdby, xpw, xpwd,
171 & xqrch, mbdt, tem,
172 & ptem, ptem1
173!
174 integer kb(im), kbcon(im), kbcon1(im),
175 & ktcon(im), ktcon1(im),
176 & jmin(im), lmin(im), kbmax(im),
177 & kbm(im), kmax(im)
178!
179 real(kind=kind_phys) ps(im), del(im,km), prsl(im,km)
180!
181 real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im),
182 & delhbar(im), delq(im), delq2(im),
183 & delqbar(im), delqev(im), deltbar(im),
184 & deltv(im), dtconv(im), edt(im),
185 & edto(im), edtx(im), fld(im),
186 & hcdo(im,km), hmax(im), hmin(im),
187 & ucdo(im,km), vcdo(im,km),aa2(im),
188 & pbcdif(im), pdot(im), po(im,km),
189 & pwavo(im), pwevo(im), xlamud(im),
190 & qcdo(im,km), qcond(im), qevap(im),
191 & rntot(im), vshear(im), xaa0(im),
192 & xk(im), xlamd(im),
193 & xmb(im), xmbmax(im), xpwav(im),
194 & xpwev(im), delubar(im),delvbar(im)
195!
196 real(kind=kind_phys) cincr, cincrmax, cincrmin
197!
198! physical parameters
199! parameter(g=grav)
200! parameter(elocp=hvap/cp,
201! & el2orc=hvap*hvap/(rv*cp))
202! parameter(c0=.002,c1=.002,delta=fv)
203! parameter(delta=fv)
204! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
205 parameter(cthk=150.,cincrmax=180.,cincrmin=120.,dthk=25.)
206!
207 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
208 & uo(im,km), vo(im,km), qeso(im,km)
209! cloud water
210! real(kind=kind_phys) tvo(im,km)
211 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
212 & dbyo(im,km), zo(im,km), xlamue(im,km),
213 & fent1(im,km), fent2(im,km), frh(im,km),
214 & heo(im,km), heso(im,km),
215 & qrcd(im,km), dellah(im,km), dellaq(im,km),
216 & dellau(im,km), dellav(im,km), hcko(im,km),
217 & ucko(im,km), vcko(im,km), qcko(im,km),
218 & eta(im,km), etad(im,km), zi(im,km),
219 & qrcko(im,km), qrcdo(im,km),
220 & pwo(im,km), pwdo(im,km),
221 & tx1(im), sumx(im), cnvwt(im,km)
222! &, rhbar(im)
223!
224 logical totflg, cnvflg(im), flg(im)
225!
226 real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
227! save pcrit, acritt
228 data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,
229 & 350.,300.,250.,200.,150./
230 data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,
231 & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
232! gdas derived acrit
233! data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
234! & .743,.813,.886,.947,1.138,1.377,1.896/
235 real(kind=kind_phys) tf, tcr, tcrf
236 parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
237!
238!-----------------------------------------------------------------------
239!************************************************************************
240! replace (derived) constants above with regular variables
241 g = grav
242 elocp = hvap/cp
243 el2orc = hvap*hvap/(rv*cp)
244 delta = fv
245 fact1 = (cvap-cliq)/rv
246 fact2 = hvap/rv-fact1*t0c
247!************************************************************************
248! initialize CCPP error handling variables
249 errmsg = ''
250 errflg = 0
251!************************************************************************
254!************************************************************************
255! convert input pa terms to cb terms -- moorthi
256 ps = psp * 0.001
257 prsl = prslp * 0.001
258 del = delp * 0.001
259!************************************************************************
260!
261!
262 km1 = km - 1
264!
265! initialize arrays
266!
267 do i=1,im
268 cnvflg(i) = .true.
269 rn(i)=0.
270 kbot(i)=km+1
271 ktop(i)=0
272 kbcon(i)=km
273 ktcon(i)=1
274 dtconv(i) = 3600.
275 cldwrk(i) = 0.
276 pdot(i) = 0.
277 pbcdif(i)= 0.
278 lmin(i) = 1
279 jmin(i) = 1
280 qlko_ktcon(i) = 0.
281 edt(i) = 0.
282 edto(i) = 0.
283 edtx(i) = 0.
284 acrt(i) = 0.
285 acrtfct(i) = 1.
286 aa1(i) = 0.
287 aa2(i) = 0.
288 xaa0(i) = 0.
289 pwavo(i)= 0.
290 pwevo(i)= 0.
291 xpwav(i)= 0.
292 xpwev(i)= 0.
293 vshear(i) = 0.
294 enddo
296 do k = 1, km
297 do i = 1, im
298 cnvw(i,k) = 0.
299 cnvc(i,k) = 0.
300 enddo
301 enddo
303! hchuang code change
304 do k = 1, km
305 do i = 1, im
306 ud_mf(i,k) = 0.
307 dd_mf(i,k) = 0.
308 dt_mf(i,k) = 0.
309 if(mp_phys == mp_phys_mg) then
310 qlcn(i,k) = 0.0
311 qicn(i,k) = 0.0
312 w_upi(i,k) = 0.0
313 cf_upi(i,k) = 0.0
314 cnv_mfd(i,k) = 0.0
315! cnv_prc3(i,k) = 0.0
316 cnv_dqldt(i,k) = 0.0
317 clcn(i,k) = 0.0
318 cnv_fice(i,k) = 0.0
319 cnv_ndrop(i,k) = 0.0
320 cnv_nice(i,k) = 0.0
321 end if
322 enddo
323 enddo
325!
326 do k = 1, 15
327 acrit(k) = acritt(k) * (975. - pcrit(k))
328 enddo
329 dt2 = delt
330 val = 1200.
331 dtmin = max(dt2, val )
332 val = 3600.
333 dtmax = max(dt2, val )
334! model tunable parameters are all here
335 mbdt = 10.
336 edtmaxl = .3
337 edtmaxs = .3
338! clam = .1
339 aafac = .1
340! betal = .15
341! betas = .15
342! betal = .05
343! betas = .05
344! evef = 0.07
345! evfact = 0.3
346! evfactl = 0.3
347!
348 cxlamu = 1.0e-4
349 xlamde = 1.0e-4
350 xlamdd = 1.0e-4
351!
352! pgcon = 0.7 ! gregory et al. (1997, qjrms)
353! pgcon = 0.55 ! zhang & wu (2003,jas)
354 fjcap = (float(jcap) / 126.) ** 2
355 val = 1.
356 fjcap = max(fjcap,val)
357 fkm = (float(km) / 28.) ** 2
358 fkm = max(fkm,val)
359 w1l = -8.e-3
360 w2l = -4.e-2
361 w3l = -5.e-3
362 w4l = -5.e-4
363 w1s = -2.e-4
364 w2s = -2.e-3
365 w3s = -1.e-3
366 w4s = -2.e-5
368!
369! define top layer for search of the downdraft originating layer
370! and the maximum thetae for updraft
371!
372 do i=1,im
373 kbmax(i) = km
374 kbm(i) = km
375 kmax(i) = km
376 tx1(i) = 1.0 / ps(i)
377 enddo
378!
379 do k = 1, km
380 do i=1,im
381 if (prsl(i,k)*tx1(i) .gt. 0.04) kmax(i) = k + 1
382 if (prsl(i,k)*tx1(i) .gt. 0.45) kbmax(i) = k + 1
383 if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
384 enddo
385 enddo
386 do i=1,im
387 kmax(i) = min(km,kmax(i))
388 kbmax(i) = min(kbmax(i),kmax(i))
389 kbm(i) = min(kbm(i),kmax(i))
390 enddo
391!
392! hydrostatic height assume zero terr and initially assume
393! updraft entrainment rate as an inverse function of height
394!
396 do k = 1, km
397 do i=1,im
398 zo(i,k) = phil(i,k) / g
399 enddo
400 enddo
402 do k = 1, km1
403 do i=1,im
404 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
405 xlamue(i,k) = clam / zi(i,k)
406 enddo
407 enddo
409!
410!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
411! convert surface pressure to mb from cb
412!
413 do k = 1, km
414 do i = 1, im
415 if (k .le. kmax(i)) then
416 pfld(i,k) = prsl(i,k) * 10.0
417 eta(i,k) = 1.
418 fent1(i,k)= 1.
419 fent2(i,k)= 1.
420 frh(i,k) = 0.
421 hcko(i,k) = 0.
422 qcko(i,k) = 0.
423 qrcko(i,k)= 0.
424 ucko(i,k) = 0.
425 vcko(i,k) = 0.
426 etad(i,k) = 1.
427 hcdo(i,k) = 0.
428 qcdo(i,k) = 0.
429 ucdo(i,k) = 0.
430 vcdo(i,k) = 0.
431 qrcd(i,k) = 0.
432 qrcdo(i,k)= 0.
433 dbyo(i,k) = 0.
434 pwo(i,k) = 0.
435 pwdo(i,k) = 0.
436 dellal(i,k) = 0.
437 to(i,k) = t1(i,k)
438 qo(i,k) = q1(i,k)
439 uo(i,k) = u1(i,k)
440 vo(i,k) = v1(i,k)
441! uo(i,k) = u1(i,k) * rcs(i)
442! vo(i,k) = v1(i,k) * rcs(i)
443 cnvwt(i,k)= 0.
444 endif
445 enddo
446 enddo
447!
448! column variables
449! p is pressure of the layer (mb)
450! t is temperature at t-dt (k)..tn
451! q is mixing ratio at t-dt (kg/kg)..qn
452! to is temperature at t+dt (k)... this is after advection and turbulan
453! qo is mixing ratio at t+dt (kg/kg)..q1
454!
456 do k = 1, km
457 do i=1,im
458 if (k .le. kmax(i)) then
459 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
460 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
461 val1 = 1.e-8
462 qeso(i,k) = max(qeso(i,k), val1)
463 val2 = 1.e-10
464 qo(i,k) = max(qo(i,k), val2 )
465! qo(i,k) = min(qo(i,k),qeso(i,k))
466! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
467 endif
468 enddo
469 enddo
470!
471! compute moist static energy
472!
474 do k = 1, km
475 do i=1,im
476 if (k .le. kmax(i)) then
477! tem = g * zo(i,k) + cp * to(i,k)
478 tem = phil(i,k) + cp * to(i,k)
479 heo(i,k) = tem + hvap * qo(i,k)
480 heso(i,k) = tem + hvap * qeso(i,k)
481! heo(i,k) = min(heo(i,k),heso(i,k))
482 endif
483 enddo
484 enddo
485
487!
488! determine level with largest moist static energy
489! this is the level where updraft starts
490!
492 do i=1,im
493 hmax(i) = heo(i,1)
494 kb(i) = 1
495 enddo
496 do k = 2, km
497 do i=1,im
498 if (k .le. kbm(i)) then
499 if(heo(i,k).gt.hmax(i)) then
500 kb(i) = k
501 hmax(i) = heo(i,k)
502 endif
503 endif
504 enddo
505 enddo
507!
508 do k = 1, km1
509 do i=1,im
510 if (k .le. kmax(i)-1) then
511 dz = .5 * (zo(i,k+1) - zo(i,k))
512 dp = .5 * (pfld(i,k+1) - pfld(i,k))
513 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
514 pprime = pfld(i,k+1) + epsm1 * es
515 qs = eps * es / pprime
516 dqsdp = - qs / pprime
517 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
518 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
519 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
520 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
521 dq = dqsdt * dt + dqsdp * dp
522 to(i,k) = to(i,k+1) + dt
523 qo(i,k) = qo(i,k+1) + dq
524 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
525 endif
526 enddo
527 enddo
529!
530 do k = 1, km1
531 do i=1,im
532 if (k .le. kmax(i)-1) then
533 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
534 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
535 val1 = 1.e-8
536 qeso(i,k) = max(qeso(i,k), val1)
537 val2 = 1.e-10
538 qo(i,k) = max(qo(i,k), val2 )
539! qo(i,k) = min(qo(i,k),qeso(i,k))
540 frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), 1.)
541 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
542 & cp * to(i,k) + hvap * qo(i,k)
543 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
544 & cp * to(i,k) + hvap * qeso(i,k)
545 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
546 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
547 endif
548 enddo
549 enddo
550!
551! look for the level of free convection as cloud base
552!
554 do i=1,im
555 flg(i) = .true.
556 kbcon(i) = kmax(i)
557 enddo
558 do k = 1, km1
559 do i=1,im
560 if (flg(i).and.k.le.kbmax(i)) then
561 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
562 kbcon(i) = k
563 flg(i) = .false.
564 endif
565 endif
566 enddo
567 enddo
569!
570 do i=1,im
571 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
572 enddo
573!!
574 totflg = .true.
575 do i=1,im
576 totflg = totflg .and. (.not. cnvflg(i))
577 enddo
578 if(totflg) return
579!!
580!
581! determine critical convective inhibition
582! as a function of vertical velocity at cloud base.
583!
585 do i=1,im
586 if(cnvflg(i)) then
587! pdot(i) = 10.* dot(i,kbcon(i))
588 pdot(i) = 0.01 * dot(i,kbcon(i)) ! now dot is in pa/s
589 endif
590 enddo
591 do i=1,im
592 if(cnvflg(i)) then
593 if(islimsk(i) == 1) then
594 w1 = w1l
595 w2 = w2l
596 w3 = w3l
597 w4 = w4l
598 else
599 w1 = w1s
600 w2 = w2s
601 w3 = w3s
602 w4 = w4s
603 endif
604 if(pdot(i).le.w4) then
605 tem = (pdot(i) - w4) / (w3 - w4)
606 elseif(pdot(i).ge.-w4) then
607 tem = - (pdot(i) + w4) / (w4 - w3)
608 else
609 tem = 0.
610 endif
611 val1 = -1.
612 tem = max(tem,val1)
613 val2 = 1.
614 tem = min(tem,val2)
615 tem = 1. - tem
616 tem1= .5*(cincrmax-cincrmin)
617 cincr = cincrmax - tem * tem1
618 pbcdif(i) = pfld(i,kb(i)) - pfld(i,kbcon(i))
619 if(pbcdif(i).gt.cincr) then
620 cnvflg(i) = .false.
621 endif
622 endif
623 enddo
624!!
625 totflg = .true.
626 do i=1,im
627 totflg = totflg .and. (.not. cnvflg(i))
628 enddo
629 if(totflg) return
630!!
631!
632! assume that updraft entrainment rate above cloud base is
633! same as that at cloud base
634!
640 do k = 2, km1
641 do i=1,im
642 if(cnvflg(i).and.
643 & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
644 xlamue(i,k) = xlamue(i,kbcon(i))
645 endif
646 enddo
647 enddo
648!
649! assume the detrainment rate for the updrafts to be same as
650! the entrainment rate at cloud base
651!
653 do i = 1, im
654 if(cnvflg(i)) then
655 xlamud(i) = xlamue(i,kbcon(i))
656 endif
657 enddo
658!
659! functions rapidly decreasing with height, mimicking a cloud ensemble
660! (bechtold et al., 2008)
661!
662 do k = 2, km1
663 do i=1,im
664 if(cnvflg(i).and.
665 & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
666 tem = qeso(i,k)/qeso(i,kbcon(i))
667 fent1(i,k) = tem**2
668 fent2(i,k) = tem**3
669 endif
670 enddo
671 enddo
672!
673! final entrainment rate as the sum of turbulent part and organized entrainment
674! depending on the environmental relative humidity
675! (bechtold et al., 2008)
676!
677 do k = 2, km1
678 do i=1,im
679 if(cnvflg(i).and.
680 & (k.ge.kbcon(i).and.k.lt.kmax(i))) then
681 tem = cxlamu * frh(i,k) * fent2(i,k)
682 xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
683 endif
684 enddo
685 enddo
686!
687! determine updraft mass flux for the subcloud layers
688!
694 do k = km1, 1, -1
695 do i = 1, im
696 if (cnvflg(i)) then
697 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
698 dz = zi(i,k+1) - zi(i,k)
699 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
700 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
701 endif
702 endif
703 enddo
704 enddo
705!
706! compute mass flux above cloud base
707!
708 do k = 2, km1
709 do i = 1, im
710 if(cnvflg(i))then
711 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
712 dz = zi(i,k) - zi(i,k-1)
713 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
714 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
715 endif
716 endif
717 enddo
718 enddo
719!
720! compute updraft cloud properties
721!
723 do i = 1, im
724 if(cnvflg(i)) then
725 indx = kb(i)
726 hcko(i,indx) = heo(i,indx)
727 ucko(i,indx) = uo(i,indx)
728 vcko(i,indx) = vo(i,indx)
729 pwavo(i) = 0.
730 endif
731 enddo
732!
733! cloud property is modified by the entrainment process
734!
736 do k = 2, km1
737 do i = 1, im
738 if (cnvflg(i)) then
739 if(k.gt.kb(i).and.k.lt.kmax(i)) then
740 dz = zi(i,k) - zi(i,k-1)
741 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
742 tem1 = 0.5 * xlamud(i) * dz
743 factor = 1. + tem - tem1
744 ptem = 0.5 * tem + pgcon
745 ptem1= 0.5 * tem - pgcon
746 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
747 & (heo(i,k)+heo(i,k-1)))/factor
748 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)
749 & +ptem1*uo(i,k-1))/factor
750 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)
751 & +ptem1*vo(i,k-1))/factor
752 dbyo(i,k) = hcko(i,k) - heso(i,k)
753 endif
754 endif
755 enddo
756 enddo
757!
758! taking account into convection inhibition due to existence of
759! dry layers below cloud base
760!
762 do i=1,im
763 flg(i) = cnvflg(i)
764 kbcon1(i) = kmax(i)
765 enddo
766 do k = 2, km1
767 do i=1,im
768 if (flg(i).and.k.lt.kmax(i)) then
769 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
770 kbcon1(i) = k
771 flg(i) = .false.
772 endif
773 endif
774 enddo
775 enddo
776 do i=1,im
777 if(cnvflg(i)) then
778 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
779 endif
780 enddo
781 do i=1,im
782 if(cnvflg(i)) then
783 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
784 if(tem.gt.dthk) then
785 cnvflg(i) = .false.
786 endif
787 endif
788 enddo
789!!
790 totflg = .true.
791 do i = 1, im
792 totflg = totflg .and. (.not. cnvflg(i))
793 enddo
794 if(totflg) return
795!!
796!
797! determine first guess cloud top as the level of zero buoyancy
798!
800 do i = 1, im
801 flg(i) = cnvflg(i)
802 ktcon(i) = 1
803 enddo
804 do k = 2, km1
805 do i = 1, im
806 if (flg(i).and.k .lt. kmax(i)) then
807 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
808 ktcon(i) = k
809 flg(i) = .false.
810 endif
811 endif
812 enddo
813 enddo
814!
815 do i = 1, im
816 if(cnvflg(i)) then
817 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
818 if(tem.lt.cthk) cnvflg(i) = .false.
819 endif
820 enddo
821!!
822 totflg = .true.
823 do i = 1, im
824 totflg = totflg .and. (.not. cnvflg(i))
825 enddo
826 if(totflg) return
827!!
828!
829! search for downdraft originating level above theta-e minimum
830!
832 do i = 1, im
833 if(cnvflg(i)) then
834 hmin(i) = heo(i,kbcon1(i))
835 lmin(i) = kbmax(i)
836 jmin(i) = kbmax(i)
837 endif
838 enddo
839 do k = 2, km1
840 do i = 1, im
841 if (cnvflg(i) .and. k .le. kbmax(i)) then
842 if(k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
843 lmin(i) = k + 1
844 hmin(i) = heo(i,k)
845 endif
846 endif
847 enddo
848 enddo
849!
850! make sure that jmin(i) is within the cloud
851!
852 do i = 1, im
853 if(cnvflg(i)) then
854 jmin(i) = min(lmin(i),ktcon(i)-1)
855 jmin(i) = max(jmin(i),kbcon1(i)+1)
856 if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
857 endif
858 enddo
859!
860! specify upper limit of mass flux at cloud base
861!
863 do i = 1, im
864 if(cnvflg(i)) then
865! xmbmax(i) = .1
866!
867 k = kbcon(i)
868 dp = 1000. * del(i,k)
869 xmbmax(i) = dp / (g * dt2)
870!
871! tem = dp / (g * dt2)
872! xmbmax(i) = min(tem, xmbmax(i))
873 endif
874 enddo
875!
876! compute cloud moisture property and precipitation
877!
879 do i = 1, im
880 if (cnvflg(i)) then
881 aa1(i) = 0.
882 qcko(i,kb(i)) = qo(i,kb(i))
883 qrcko(i,kb(i)) = qo(i,kb(i))
884! rhbar(i) = 0.
885 endif
886 enddo
888 do k = 2, km1
889 do i = 1, im
890 if (cnvflg(i)) then
891 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
892 dz = zi(i,k) - zi(i,k-1)
893 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
894 qrch = qeso(i,k)
895 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
896!
897 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
898 tem1 = 0.5 * xlamud(i) * dz
899 factor = 1. + tem - tem1
900 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
901 & (qo(i,k)+qo(i,k-1)))/factor
902 qrcko(i,k) = qcko(i,k)
903!
904 dq = eta(i,k) * (qcko(i,k) - qrch)
905!
906! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
907!
908! check if there is excess moisture to release latent heat
909!
910 if(k.ge.kbcon(i).and.dq.gt.0.) then
911 etah = .5 * (eta(i,k) + eta(i,k-1))
912 dp = 1000. * del(i,k)
913 if(ncloud.gt.0..and.k.gt.jmin(i)) then
914 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
915 dellal(i,k) = etah * c1 * dz * qlk * g / dp
916 else
917 qlk = dq / (eta(i,k) + etah * c0 * dz)
918 endif
919 aa1(i) = aa1(i) - dz * g * qlk
920 qcko(i,k) = qlk + qrch
921 pwo(i,k) = etah * c0 * dz * qlk
922 pwavo(i) = pwavo(i) + pwo(i,k)
923! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
924 cnvwt(i,k) = etah * qlk * g / dp
925 endif
926 endif
927 endif
928 enddo
929 enddo
930!
931! do i = 1, im
932! if(cnvflg(i)) then
933! indx = ktcon(i) - kb(i) - 1
934! rhbar(i) = rhbar(i) / float(indx)
935! endif
936! enddo
937!
938! calculate cloud work function
939!
945 do k = 2, km1
946 do i = 1, im
947 if (cnvflg(i)) then
948 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
949 dz1 = zo(i,k+1) - zo(i,k)
950 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
951 rfact = 1. + delta * cp * gamma
952 & * to(i,k) / hvap
953 aa1(i) = aa1(i) +
954 & dz1 * (g / (cp * to(i,k)))
955 & * dbyo(i,k) / (1. + gamma)
956 & * rfact
957 val = 0.
958 aa1(i)=aa1(i)+
959 & dz1 * g * delta *
960 & max(val,(qeso(i,k) - qo(i,k)))
961 endif
962 endif
963 enddo
964 enddo
966 do i = 1, im
967 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
968 enddo
969!!
970 totflg = .true.
971 do i=1,im
972 totflg = totflg .and. (.not. cnvflg(i))
973 enddo
974 if(totflg) return
975!!
976!
977! estimate the onvective overshooting as the level
978! where the [aafac * cloud work function] becomes zero,
979! which is the final cloud top
980!
982 do i = 1, im
983 if (cnvflg(i)) then
984 aa2(i) = aafac * aa1(i)
985 endif
986 enddo
987!
988 do i = 1, im
989 flg(i) = cnvflg(i)
990 ktcon1(i) = ktcon(i)
991 enddo
992 do k = 2, km1
993 do i = 1, im
994 if (flg(i)) then
995 if(k.ge.ktcon(i).and.k.lt.kmax(i)) then
996 dz1 = zo(i,k+1) - zo(i,k)
997 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
998 rfact = 1. + delta * cp * gamma
999 & * to(i,k) / hvap
1000 aa2(i) = aa2(i) +
1001 & dz1 * (g / (cp * to(i,k)))
1002 & * dbyo(i,k) / (1. + gamma)
1003 & * rfact
1004!NRL MNM: Limit overshooting not to be deeper than the actual cloud
1005 tem = 0.5 * (zi(i,ktcon(i))-zi(i,kbcon(i)))
1006 tem1 = zi(i,k)-zi(i,ktcon(i))
1007 if(aa2(i) < 0. .or. tem1 >= tem) then
1008 ktcon1(i) = k
1009 flg(i) = .false.
1010 endif
1011 endif
1012 endif
1013 enddo
1014 enddo
1015!
1016! compute cloud moisture property, detraining cloud water
1017! and precipitation in overshooting layers
1018!
1020 do k = 2, km1
1021 do i = 1, im
1022 if (cnvflg(i)) then
1023 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
1024 dz = zi(i,k) - zi(i,k-1)
1025 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1026 qrch = qeso(i,k)
1027 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1028!
1029 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1030 tem1 = 0.5 * xlamud(i) * dz
1031 factor = 1. + tem - tem1
1032 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
1033 & (qo(i,k)+qo(i,k-1)))/factor
1034 qrcko(i,k) = qcko(i,k)
1035!
1036 dq = eta(i,k) * (qcko(i,k) - qrch)
1037!
1038! check if there is excess moisture to release latent heat
1039!
1040 if(dq.gt.0.) then
1041 etah = .5 * (eta(i,k) + eta(i,k-1))
1042 dp = 1000. * del(i,k)
1043 if(ncloud.gt.0.) then
1044 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
1045 dellal(i,k) = etah * c1 * dz * qlk * g / dp
1046 else
1047 qlk = dq / (eta(i,k) + etah * c0 * dz)
1048 endif
1049 qcko(i,k) = qlk + qrch
1050 pwo(i,k) = etah * c0 * dz * qlk
1051 pwavo(i) = pwavo(i) + pwo(i,k)
1052! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
1053 cnvwt(i,k) = etah * qlk * g / dp
1054 endif
1055 endif
1056 endif
1057 enddo
1058 enddo
1059!
1060! exchange ktcon with ktcon1
1061!
1063 do i = 1, im
1064 if(cnvflg(i)) then
1065 kk = ktcon(i)
1066 ktcon(i) = ktcon1(i)
1067 ktcon1(i) = kk
1068 endif
1069 enddo
1070!
1071! this section is ready for cloud water
1072!
1074 if(ncloud.gt.0) then
1075!
1076! compute liquid and vapor separation at cloud top
1077!
1078 do i = 1, im
1079 if(cnvflg(i)) then
1080 k = ktcon(i) - 1
1081 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1082 qrch = qeso(i,k)
1083 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1084 dq = qcko(i,k) - qrch
1085!
1086! check if there is excess moisture to release latent heat
1087!
1088 if(dq.gt.0.) then
1089 qlko_ktcon(i) = dq
1090 qcko(i,k) = qrch
1091 endif
1092 endif
1093 enddo
1094 endif
1095!
1096! if(lat.eq.latd.and.lon.eq.lond.and.cnvflg(i)) then
1097! print *, ' aa1(i) before dwndrft =', aa1(i)
1098! endif
1099!
1100!------- downdraft calculations
1101!
1102!--- compute precipitation efficiency in terms of windshear
1103!
1110 do i = 1, im
1111 if(cnvflg(i)) then
1112 vshear(i) = 0.
1113 endif
1114 enddo
1115 do k = 2, km
1116 do i = 1, im
1117 if (cnvflg(i)) then
1118 if(k.gt.kb(i).and.k.le.ktcon(i)) then
1119 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
1120 & + (vo(i,k)-vo(i,k-1)) ** 2)
1121 vshear(i) = vshear(i) + shear
1122 endif
1123 endif
1124 enddo
1125 enddo
1126 do i = 1, im
1127 if(cnvflg(i)) then
1128 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
1129 e1=1.591-.639*vshear(i)
1130 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1131 edt(i)=1.-e1
1132 val = .9
1133 edt(i) = min(edt(i),val)
1134 val = .0
1135 edt(i) = max(edt(i),val)
1136 edto(i)=edt(i)
1137 edtx(i)=edt(i)
1138 endif
1139 enddo
1140!
1141! determine detrainment rate between 1 and kbcon
1142!
1148 do i = 1, im
1149 if(cnvflg(i)) then
1150 sumx(i) = 0.
1151 endif
1152 enddo
1153 do k = 1, km1
1154 do i = 1, im
1155 if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
1156 dz = zi(i,k+1) - zi(i,k)
1157 sumx(i) = sumx(i) + dz
1158 endif
1159 enddo
1160 enddo
1161 do i = 1, im
1162 beta = betas
1163 if(islimsk(i) == 1) beta = betal
1164 if(cnvflg(i)) then
1165 dz = (sumx(i)+zi(i,1))/float(kbcon(i))
1166 tem = 1./float(kbcon(i))
1167 xlamd(i) = (1.-beta**tem)/dz
1168 endif
1169 enddo
1170!
1171! determine downdraft mass flux
1172!
1174 do k = km1, 1, -1
1175 do i = 1, im
1176 if (cnvflg(i) .and. k .le. kmax(i)-1) then
1177 if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
1178 dz = zi(i,k+1) - zi(i,k)
1179 ptem = xlamdd - xlamde
1180 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1181 else if(k.lt.kbcon(i)) then
1182 dz = zi(i,k+1) - zi(i,k)
1183 ptem = xlamd(i) + xlamdd - xlamde
1184 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1185 endif
1186 endif
1187 enddo
1188 enddo
1189!
1190!--- downdraft moisture properties
1191!
1193 do i = 1, im
1194 if(cnvflg(i)) then
1195 jmn = jmin(i)
1196 hcdo(i,jmn) = heo(i,jmn)
1197 qcdo(i,jmn) = qo(i,jmn)
1198 qrcdo(i,jmn)= qo(i,jmn)
1199 ucdo(i,jmn) = uo(i,jmn)
1200 vcdo(i,jmn) = vo(i,jmn)
1201 pwevo(i) = 0.
1202 endif
1203 enddo
1204!j
1206 do k = km1, 1, -1
1207 do i = 1, im
1208 if (cnvflg(i) .and. k.lt.jmin(i)) then
1209 dz = zi(i,k+1) - zi(i,k)
1210 if(k.ge.kbcon(i)) then
1211 tem = xlamde * dz
1212 tem1 = 0.5 * xlamdd * dz
1213 else
1214 tem = xlamde * dz
1215 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1216 endif
1217 factor = 1. + tem - tem1
1218 ptem = 0.5 * tem - pgcon
1219 ptem1= 0.5 * tem + pgcon
1220 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
1221 & (heo(i,k)+heo(i,k+1)))/factor
1222 ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1)
1223 & +ptem1*uo(i,k))/factor
1224 vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1)
1225 & +ptem1*vo(i,k))/factor
1226 dbyo(i,k) = hcdo(i,k) - heso(i,k)
1227 endif
1228 enddo
1229 enddo
1230!
1232 do k = km1, 1, -1
1233 do i = 1, im
1234 if (cnvflg(i).and.k.lt.jmin(i)) then
1235 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1236 qrcdo(i,k) = qeso(i,k)+
1237 & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
1238! detad = etad(i,k+1) - etad(i,k)
1239!
1240 dz = zi(i,k+1) - zi(i,k)
1241 if(k.ge.kbcon(i)) then
1242 tem = xlamde * dz
1243 tem1 = 0.5 * xlamdd * dz
1244 else
1245 tem = xlamde * dz
1246 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1247 endif
1248 factor = 1. + tem - tem1
1249 qcdo(i,k) = ((1.-tem1)*qrcdo(i,k+1)+tem*0.5*
1250 & (qo(i,k)+qo(i,k+1)))/factor
1251!
1252! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
1253! & etad(i,k) * qrcdo(i,k)
1254! pwdo(i,k) = pwdo(i,k) - detad *
1255! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
1256!
1257 pwdo(i,k) = etad(i,k) * (qcdo(i,k) - qrcdo(i,k))
1258 pwevo(i) = pwevo(i) + pwdo(i,k)
1259 endif
1260 enddo
1261 enddo
1262!
1263!--- final downdraft strength dependent on precip
1264!--- efficiency (edt), normalized condensate (pwav), and
1265!--- evaporate (pwev)
1266!
1268 do i = 1, im
1269 edtmax = edtmaxl
1270 if(islimsk(i) == 0) edtmax = edtmaxs
1271 if(cnvflg(i)) then
1272 if(pwevo(i).lt.0.) then
1273 edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1274 edto(i) = min(edto(i),edtmax)
1275 else
1276 edto(i) = 0.
1277 endif
1278 endif
1279 enddo
1280!
1281!--- downdraft cloudwork functions
1282!
1284 do k = km1, 1, -1
1285 do i = 1, im
1286 if (cnvflg(i) .and. k .lt. jmin(i)) then
1287 gamma = el2orc * qeso(i,k) / to(i,k)**2
1288 dhh=hcdo(i,k)
1289 dt=to(i,k)
1290 dg=gamma
1291 dh=heso(i,k)
1292 dz=-1.*(zo(i,k+1)-zo(i,k))
1293 aa1(i)=aa1(i)+edto(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg))
1294 & *(1.+delta*cp*dg*dt/hvap)
1295 val=0.
1296 aa1(i)=aa1(i)+edto(i)*
1297 & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
1298 endif
1299 enddo
1300 enddo
1302 do i = 1, im
1303 if(cnvflg(i).and.aa1(i).le.0.) then
1304 cnvflg(i) = .false.
1305 endif
1306 enddo
1307!!
1308 totflg = .true.
1309 do i=1,im
1310 totflg = totflg .and. (.not. cnvflg(i))
1311 enddo
1312 if(totflg) return
1313!!
1314!
1315!--- what would the change be, that a cloud with unit mass
1316!--- will do to the environment?
1317!
1319 do k = 1, km
1320 do i = 1, im
1321 if(cnvflg(i) .and. k .le. kmax(i)) then
1322 dellah(i,k) = 0.
1323 dellaq(i,k) = 0.
1324 dellau(i,k) = 0.
1325 dellav(i,k) = 0.
1326 endif
1327 enddo
1328 enddo
1329 do i = 1, im
1330 if(cnvflg(i)) then
1331 dp = 1000. * del(i,1)
1332 dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1)
1333 & - heo(i,1)) * g / dp
1334 dellaq(i,1) = edto(i) * etad(i,1) * (qrcdo(i,1)
1335 & - qo(i,1)) * g / dp
1336 dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1)
1337 & - uo(i,1)) * g / dp
1338 dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1)
1339 & - vo(i,1)) * g / dp
1340 endif
1341 enddo
1342!
1343!--- changed due to subsidence and entrainment
1344!
1345 do k = 2, km1
1346 do i = 1, im
1347 if (cnvflg(i).and.k.lt.ktcon(i)) then
1348 aup = 1.
1349 if(k.le.kb(i)) aup = 0.
1350 adw = 1.
1351 if(k.gt.jmin(i)) adw = 0.
1352 dp = 1000. * del(i,k)
1353 dz = zi(i,k) - zi(i,k-1)
1354!
1355 dv1h = heo(i,k)
1356 dv2h = .5 * (heo(i,k) + heo(i,k-1))
1357 dv3h = heo(i,k-1)
1358 dv1q = qo(i,k)
1359 dv2q = .5 * (qo(i,k) + qo(i,k-1))
1360 dv3q = qo(i,k-1)
1361 dv1u = uo(i,k)
1362 dv2u = .5 * (uo(i,k) + uo(i,k-1))
1363 dv3u = uo(i,k-1)
1364 dv1v = vo(i,k)
1365 dv2v = .5 * (vo(i,k) + vo(i,k-1))
1366 dv3v = vo(i,k-1)
1367!
1368 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1369 tem1 = xlamud(i)
1370!
1371 if(k.le.kbcon(i)) then
1372 ptem = xlamde
1373 ptem1 = xlamd(i)+xlamdd
1374 else
1375 ptem = xlamde
1376 ptem1 = xlamdd
1377 endif
1378!
1379 dellah(i,k) = dellah(i,k) +
1380 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h
1381 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h
1382 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz
1383 & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1384 & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz
1385 & ) *g/dp
1386!
1387 dellaq(i,k) = dellaq(i,k) +
1388 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q
1389 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q
1390 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz
1391 & + aup*tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz
1392 & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qcdo(i,k-1))*dz
1393 & ) *g/dp
1394!
1395 dellau(i,k) = dellau(i,k) +
1396 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1u
1397 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3u
1398 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz
1399 & + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz
1400 & + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz
1401 & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u)
1402 & ) *g/dp
1403!
1404 dellav(i,k) = dellav(i,k) +
1405 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1v
1406 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3v
1407 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz
1408 & + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz
1409 & + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz
1410 & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v)
1411 & ) *g/dp
1412!
1413 endif
1414 enddo
1415 enddo
1416!
1417!------- cloud top
1418!
1419 do i = 1, im
1420 if(cnvflg(i)) then
1421 indx = ktcon(i)
1422 dp = 1000. * del(i,indx)
1423 dv1h = heo(i,indx-1)
1424 dellah(i,indx) = eta(i,indx-1) *
1425 & (hcko(i,indx-1) - dv1h) * g / dp
1426 dv1q = qo(i,indx-1)
1427 dellaq(i,indx) = eta(i,indx-1) *
1428 & (qcko(i,indx-1) - dv1q) * g / dp
1429 dv1u = uo(i,indx-1)
1430 dellau(i,indx) = eta(i,indx-1) *
1431 & (ucko(i,indx-1) - dv1u) * g / dp
1432 dv1v = vo(i,indx-1)
1433 dellav(i,indx) = eta(i,indx-1) *
1434 & (vcko(i,indx-1) - dv1v) * g / dp
1435!
1436! cloud water
1437!
1438 dellal(i,indx) = eta(i,indx-1) *
1439 & qlko_ktcon(i) * g / dp
1440 endif
1441 enddo
1442!
1443!------- final changed variable per unit mass flux
1444!
1446 do k = 1, km
1447 do i = 1, im
1448 if (cnvflg(i).and.k .le. kmax(i)) then
1449 if(k.gt.ktcon(i)) then
1450 qo(i,k) = q1(i,k)
1451 to(i,k) = t1(i,k)
1452 endif
1453 if(k.le.ktcon(i)) then
1454 qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
1455 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1456 to(i,k) = dellat * mbdt + t1(i,k)
1457 val = 1.e-10
1458 qo(i,k) = max(qo(i,k), val )
1459 endif
1460 endif
1461 enddo
1462 enddo
1463!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1464!
1465!--- the above changed environment is now used to calulate the
1466!--- effect the arbitrary cloud (with unit mass flux)
1467!--- would have on the stability,
1468!--- which then is used to calculate the real mass flux,
1469!--- necessary to keep this change in balance with the large-scale
1470!--- destabilization.
1471!
1472!--- environmental conditions again, first heights
1473!
1477 do k = 1, km
1478 do i = 1, im
1479 if(cnvflg(i) .and. k .le. kmax(i)) then
1480 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1481 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
1482 val = 1.e-8
1483 qeso(i,k) = max(qeso(i,k), val )
1484! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
1485 endif
1486 enddo
1487 enddo
1488!
1489!--- moist static energy
1490!
1491!! - Recalculate moist static energy and saturation moist static energy.
1492 do k = 1, km1
1493 do i = 1, im
1494 if(cnvflg(i) .and. k .le. kmax(i)-1) then
1495 dz = .5 * (zo(i,k+1) - zo(i,k))
1496 dp = .5 * (pfld(i,k+1) - pfld(i,k))
1497 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
1498 pprime = pfld(i,k+1) + epsm1 * es
1499 qs = eps * es / pprime
1500 dqsdp = - qs / pprime
1501 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1502 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
1503 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1504 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
1505 dq = dqsdt * dt + dqsdp * dp
1506 to(i,k) = to(i,k+1) + dt
1507 qo(i,k) = qo(i,k+1) + dq
1508 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
1509 endif
1510 enddo
1511 enddo
1512 do k = 1, km1
1513 do i = 1, im
1514 if(cnvflg(i) .and. k .le. kmax(i)-1) then
1515 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1516 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
1517 val1 = 1.e-8
1518 qeso(i,k) = max(qeso(i,k), val1)
1519 val2 = 1.e-10
1520 qo(i,k) = max(qo(i,k), val2 )
1521! qo(i,k) = min(qo(i,k),qeso(i,k))
1522 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
1523 & cp * to(i,k) + hvap * qo(i,k)
1524 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
1525 & cp * to(i,k) + hvap * qeso(i,k)
1526 endif
1527 enddo
1528 enddo
1529 do i = 1, im
1530 if(cnvflg(i)) then
1531 k = kmax(i)
1532 heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
1533 heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
1534! heo(i,k) = min(heo(i,k),heso(i,k))
1535 endif
1536 enddo
1537!
1538!**************************** static control
1539!
1540!------- moisture and cloud work functions
1541!
1543 do i = 1, im
1544 if(cnvflg(i)) then
1545 xaa0(i) = 0.
1546 xpwav(i) = 0.
1547 endif
1548 enddo
1549!
1550 do i = 1, im
1551 if(cnvflg(i)) then
1552 indx = kb(i)
1553 hcko(i,indx) = heo(i,indx)
1554 qcko(i,indx) = qo(i,indx)
1555 endif
1556 enddo
1557 do k = 2, km1
1558 do i = 1, im
1559 if (cnvflg(i)) then
1560 if(k.gt.kb(i).and.k.le.ktcon(i)) then
1561 dz = zi(i,k) - zi(i,k-1)
1562 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1563 tem1 = 0.5 * xlamud(i) * dz
1564 factor = 1. + tem - tem1
1565 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
1566 & (heo(i,k)+heo(i,k-1)))/factor
1567 endif
1568 endif
1569 enddo
1570 enddo
1571 do k = 2, km1
1572 do i = 1, im
1573 if (cnvflg(i)) then
1574 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
1575 dz = zi(i,k) - zi(i,k-1)
1576 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1577 xdby = hcko(i,k) - heso(i,k)
1578 xqrch = qeso(i,k)
1579 & + gamma * xdby / (hvap * (1. + gamma))
1580!
1581 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1582 tem1 = 0.5 * xlamud(i) * dz
1583 factor = 1. + tem - tem1
1584 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
1585 & (qo(i,k)+qo(i,k-1)))/factor
1586!
1587 dq = eta(i,k) * (qcko(i,k) - xqrch)
1588!
1589 if(k.ge.kbcon(i).and.dq.gt.0.) then
1590 etah = .5 * (eta(i,k) + eta(i,k-1))
1591 if(ncloud.gt.0..and.k.gt.jmin(i)) then
1592 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
1593 else
1594 qlk = dq / (eta(i,k) + etah * c0 * dz)
1595 endif
1596 if(k.lt.ktcon1(i)) then
1597 xaa0(i) = xaa0(i) - dz * g * qlk
1598 endif
1599 qcko(i,k) = qlk + xqrch
1600 xpw = etah * c0 * dz * qlk
1601 xpwav(i) = xpwav(i) + xpw
1602 endif
1603 endif
1604 if(k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
1605 dz1 = zo(i,k+1) - zo(i,k)
1606 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1607 rfact = 1. + delta * cp * gamma
1608 & * to(i,k) / hvap
1609 xaa0(i) = xaa0(i)
1610 & + dz1 * (g / (cp * to(i,k)))
1611 & * xdby / (1. + gamma)
1612 & * rfact
1613 val=0.
1614 xaa0(i)=xaa0(i)+
1615 & dz1 * g * delta *
1616 & max(val,(qeso(i,k) - qo(i,k)))
1617 endif
1618 endif
1619 enddo
1620 enddo
1621!
1622!------- downdraft calculations
1623!
1624!--- downdraft moisture properties
1625!
1627 do i = 1, im
1628 if(cnvflg(i)) then
1629 jmn = jmin(i)
1630 hcdo(i,jmn) = heo(i,jmn)
1631 qcdo(i,jmn) = qo(i,jmn)
1632 qrcd(i,jmn) = qo(i,jmn)
1633 xpwev(i) = 0.
1634 endif
1635 enddo
1636!
1637 do k = km1, 1, -1
1638 do i = 1, im
1639 if (cnvflg(i) .and. k.lt.jmin(i)) then
1640 dz = zi(i,k+1) - zi(i,k)
1641 if(k.ge.kbcon(i)) then
1642 tem = xlamde * dz
1643 tem1 = 0.5 * xlamdd * dz
1644 else
1645 tem = xlamde * dz
1646 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1647 endif
1648 factor = 1. + tem - tem1
1649 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
1650 & (heo(i,k)+heo(i,k+1)))/factor
1651 endif
1652 enddo
1653 enddo
1654!
1655 do k = km1, 1, -1
1656 do i = 1, im
1657 if (cnvflg(i) .and. k .lt. jmin(i)) then
1658 dq = qeso(i,k)
1659 dt = to(i,k)
1660 gamma = el2orc * dq / dt**2
1661 dh = hcdo(i,k) - heso(i,k)
1662 qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
1663! detad = etad(i,k+1) - etad(i,k)
1664!
1665 dz = zi(i,k+1) - zi(i,k)
1666 if(k.ge.kbcon(i)) then
1667 tem = xlamde * dz
1668 tem1 = 0.5 * xlamdd * dz
1669 else
1670 tem = xlamde * dz
1671 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1672 endif
1673 factor = 1. + tem - tem1
1674 qcdo(i,k) = ((1.-tem1)*qrcd(i,k+1)+tem*0.5*
1675 & (qo(i,k)+qo(i,k+1)))/factor
1676!
1677! xpwd = etad(i,k+1) * qcdo(i,k+1) -
1678! & etad(i,k) * qrcd(i,k)
1679! xpwd = xpwd - detad *
1680! & .5 * (qrcd(i,k) + qrcd(i,k+1))
1681!
1682 xpwd = etad(i,k) * (qcdo(i,k) - qrcd(i,k))
1683 xpwev(i) = xpwev(i) + xpwd
1684 endif
1685 enddo
1686 enddo
1687!
1688 do i = 1, im
1689 edtmax = edtmaxl
1690 if(islimsk(i) == 0) edtmax = edtmaxs
1691 if(cnvflg(i)) then
1692 if(xpwev(i).ge.0.) then
1693 edtx(i) = 0.
1694 else
1695 edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1696 edtx(i) = min(edtx(i),edtmax)
1697 endif
1698 endif
1699 enddo
1700!
1701!
1702!--- downdraft cloudwork functions
1703!
1704!
1705 do k = km1, 1, -1
1706 do i = 1, im
1707 if (cnvflg(i) .and. k.lt.jmin(i)) then
1708 gamma = el2orc * qeso(i,k) / to(i,k)**2
1709 dhh=hcdo(i,k)
1710 dt= to(i,k)
1711 dg= gamma
1712 dh= heso(i,k)
1713 dz=-1.*(zo(i,k+1)-zo(i,k))
1714 xaa0(i)=xaa0(i)+edtx(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg))
1715 & *(1.+delta*cp*dg*dt/hvap)
1716 val=0.
1717 xaa0(i)=xaa0(i)+edtx(i)*
1718 & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
1719 endif
1720 enddo
1721 enddo
1722!
1723! calculate critical cloud work function
1724!
1727 do i = 1, im
1728 if(cnvflg(i)) then
1729 if(pfld(i,ktcon(i)).lt.pcrit(15))then
1730 acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i)))
1731 & /(975.-pcrit(15))
1732 else if(pfld(i,ktcon(i)).gt.pcrit(1))then
1733 acrt(i)=acrit(1)
1734 else
1735 k = int((850. - pfld(i,ktcon(i)))/50.) + 2
1736 k = min(k,15)
1737 k = max(k,2)
1738 acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*
1739 & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
1740 endif
1741 endif
1742 enddo
1744 do i = 1, im
1745 if(cnvflg(i)) then
1746 if(islimsk(i) == 1) then
1747 w1 = w1l
1748 w2 = w2l
1749 w3 = w3l
1750 w4 = w4l
1751 else
1752 w1 = w1s
1753 w2 = w2s
1754 w3 = w3s
1755 w4 = w4s
1756 endif
1757!
1758! modify critical cloud workfunction by cloud base vertical velocity
1759!
1760 if(pdot(i).le.w4) then
1761 acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
1762 elseif(pdot(i).ge.-w4) then
1763 acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
1764 else
1765 acrtfct(i) = 0.
1766 endif
1767 val1 = -1.
1768 acrtfct(i) = max(acrtfct(i),val1)
1769 val2 = 1.
1770 acrtfct(i) = min(acrtfct(i),val2)
1771 acrtfct(i) = 1. - acrtfct(i)
1772!
1773! modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
1774!
1775! if(rhbar(i).ge..8) then
1776! acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
1777! endif
1778!
1779! modify adjustment time scale by cloud base vertical velocity
1780!
1782 dtconv(i) = dt2 + max((1800. - dt2),0.) *
1783 & (pdot(i) - w2) / (w1 - w2)
1784! dtconv(i) = max(dtconv(i), dt2)
1785! dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
1786 dtconv(i) = max(dtconv(i),dtmin)
1787 dtconv(i) = min(dtconv(i),dtmax)
1788!
1789 endif
1790 enddo
1791!
1792!--- large scale forcing
1793!
1799 do i= 1, im
1800 if(cnvflg(i)) then
1801 fld(i)=(aa1(i)-acrt(i)* acrtfct(i))/dtconv(i)
1802 if(fld(i).le.0.) cnvflg(i) = .false.
1803 endif
1809 if(cnvflg(i)) then
1810! xaa0(i) = max(xaa0(i),0.)
1811 xk(i) = (xaa0(i) - aa1(i)) / mbdt
1812 if(xk(i).ge.0.) cnvflg(i) = .false.
1813 endif
1814!
1815!--- kernel, cloud base mass flux
1816!
1821 if(cnvflg(i)) then
1822 xmb(i) = -fld(i) / xk(i)
1823 xmb(i) = min(xmb(i),xmbmax(i))
1824 endif
1825 enddo
1826!!
1828 totflg = .true.
1829 do i=1,im
1830 totflg = totflg .and. (.not. cnvflg(i))
1831 enddo
1832 if(totflg) return
1833!!
1834!
1835! restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
1836!
1837
1838 do k = 1, km
1839 do i = 1, im
1840 if (cnvflg(i) .and. k .le. kmax(i)) then
1841 to(i,k) = t1(i,k)
1842 qo(i,k) = q1(i,k)
1843 uo(i,k) = u1(i,k)
1844 vo(i,k) = v1(i,k)
1845 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1846 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
1847 val = 1.e-8
1848 qeso(i,k) = max(qeso(i,k), val )
1849 endif
1850 enddo
1851 enddo
1852!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1853!
1854!--- feedback: simply the changes from the cloud with unit mass flux
1855!--- multiplied by the mass flux necessary to keep the
1856!--- equilibrium with the larger-scale.
1857!
1862 do i = 1, im
1863 delhbar(i) = 0.
1864 delqbar(i) = 0.
1865 deltbar(i) = 0.
1866 delubar(i) = 0.
1867 delvbar(i) = 0.
1868 qcond(i) = 0.
1869 enddo
1870 do k = 1, km
1871 do i = 1, im
1872 if (cnvflg(i) .and. k .le. kmax(i)) then
1873 if(k.le.ktcon(i)) then
1874 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1875 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1876 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1877! tem = 1./rcs(i)
1878! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1879! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1880 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
1881 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
1882 dp = 1000. * del(i,k)
1883 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
1884 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
1885 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
1886 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
1887 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
1888 endif
1889 endif
1890 enddo
1891 enddo
1893 do k = 1, km
1894 do i = 1, im
1895 if (cnvflg(i) .and. k .le. kmax(i)) then
1896 if(k.le.ktcon(i)) then
1897 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1898 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
1899 val = 1.e-8
1900 qeso(i,k) = max(qeso(i,k), val )
1901 endif
1902 endif
1903 enddo
1904 enddo
1905!
1907 do i = 1, im
1908 rntot(i) = 0.
1909 delqev(i) = 0.
1910 delq2(i) = 0.
1911 flg(i) = cnvflg(i)
1912 enddo
1913 do k = km, 1, -1
1914 do i = 1, im
1915 if (cnvflg(i) .and. k .le. kmax(i)) then
1916 if(k.lt.ktcon(i)) then
1917 aup = 1.
1918 if(k.le.kb(i)) aup = 0.
1919 adw = 1.
1920 if(k.ge.jmin(i)) adw = 0.
1921 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
1922 rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
1923 endif
1924 endif
1925 enddo
1926 enddo
1930 do k = km, 1, -1
1931 do i = 1, im
1932 if (k .le. kmax(i)) then
1933 deltv(i) = 0.
1934 delq(i) = 0.
1935 qevap(i) = 0.
1936 if(cnvflg(i).and.k.lt.ktcon(i)) then
1937 aup = 1.
1938 if(k.le.kb(i)) aup = 0.
1939 adw = 1.
1940 if(k.ge.jmin(i)) adw = 0.
1941 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
1942 rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
1943 endif
1944 if(flg(i).and.k.lt.ktcon(i)) then
1945 evef = edt(i) * evfact
1946 if(islimsk(i) == 1) evef=edt(i) * evfactl
1947! if(islimsk(i) == 1) evef=.07
1948! if(islimsk(i) == 1) evef = 0.
1949 qcond(i) = evef * (q1(i,k) - qeso(i,k))
1950 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
1951 dp = 1000. * del(i,k)
1952 if(rn(i).gt.0..and.qcond(i).lt.0.) then
1953 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
1954 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
1955 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
1956 endif
1957 if(rn(i).gt.0..and.qcond(i).lt.0..and.
1958 & delq2(i).gt.rntot(i)) then
1959 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
1960 flg(i) = .false.
1961 endif
1962 if(rn(i).gt.0..and.qevap(i).gt.0.) then
1963 q1(i,k) = q1(i,k) + qevap(i)
1964 t1(i,k) = t1(i,k) - elocp * qevap(i)
1965 rn(i) = rn(i) - .001 * qevap(i) * dp / g
1966 deltv(i) = - elocp*qevap(i)/dt2
1967 delq(i) = + qevap(i)/dt2
1968 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
1969 endif
1970 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
1971 delqbar(i) = delqbar(i) + delq(i)*dp/g
1972 deltbar(i) = deltbar(i) + deltv(i)*dp/g
1973 endif
1974 endif
1975 enddo
1976 enddo
1977!
1978! do i = 1, im
1979! if(me.eq.31.and.cnvflg(i)) then
1980! if(cnvflg(i)) then
1981! print *, ' deep delhbar, delqbar, deltbar = ',
1982! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
1983! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
1984! print *, ' precip =', hvap*rn(i)*1000./dt2
1985! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
1986! endif
1987! enddo
1988!
1989! precipitation rate converted to actual precip
1990! in unit of m instead of kg
1991!
1992 do i = 1, im
1993 if(cnvflg(i)) then
1994!
1995! in the event of upper level rain evaporation and lower level downdraft
1996! moistening, rn can become negative, in this case, we back out of the
1997! heating and the moistening
1998!
1999
2000 if(rn(i).lt.0..and..not.flg(i)) rn(i) = 0.
2001 if(rn(i).le.0.) then
2002 rn(i) = 0.
2003 else
2004 ktop(i) = ktcon(i)
2005 kbot(i) = kbcon(i)
2006 kcnv(i) = 1
2007 cldwrk(i) = aa1(i)
2008 endif
2009 endif
2010 enddo
2011!
2012! convective cloud water
2013!
2015 do k = 1, km
2016 do i = 1, im
2017 if (cnvflg(i) .and. rn(i).gt.0.) then
2018 if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
2019 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
2020 endif
2021 endif
2022 enddo
2023 enddo
2024!
2025! convective cloud cover
2026!
2028 do k = 1, km
2029 do i = 1, im
2030 if (cnvflg(i) .and. rn(i).gt.0.) then
2031 if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
2032 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
2033 cnvc(i,k) = min(cnvc(i,k), 0.6)
2034 cnvc(i,k) = max(cnvc(i,k), 0.0)
2035 endif
2036 endif
2037 enddo
2038 enddo
2039
2040!
2041! cloud water
2042!
2044 if (ncloud.gt.0) then
2045!
2046 do k = 1, km
2047 do i = 1, im
2048 if (cnvflg(i) .and. rn(i).gt.0.) then
2049 if (k.gt.kb(i).and.k.le.ktcon(i)) then
2050 tem = dellal(i,k) * xmb(i) * dt2
2051 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2052 if (qlc(i,k) .gt. -999.0) then
2053 qli(i,k) = qli(i,k) + tem * tem1 ! ice
2054 qlc(i,k) = qlc(i,k) + tem *(1.0-tem1) ! water
2055 else
2056 qli(i,k) = qli(i,k) + tem
2057 endif
2058 endif
2059 endif
2060 enddo
2061 enddo
2062!
2063 endif
2064!
2066 do k = 1, km
2067 do i = 1, im
2068 if(cnvflg(i).and.rn(i).le.0.) then
2069 if (k .le. kmax(i)) then
2070 t1(i,k) = to(i,k)
2071 q1(i,k) = qo(i,k)
2072 u1(i,k) = uo(i,k)
2073 v1(i,k) = vo(i,k)
2074 endif
2075 endif
2076 enddo
2077 enddo
2078!
2079! hchuang code change
2080!
2082 do k = 1, km
2083 do i = 1, im
2084 if(cnvflg(i).and.rn(i).gt.0.) then
2085 if(k.ge.kb(i) .and. k.lt.ktop(i)) then
2086 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
2087 endif
2088 endif
2089 enddo
2090 enddo
2092 do i = 1, im
2093 if(cnvflg(i).and.rn(i).gt.0.) then
2094 k = ktop(i)-1
2095 dt_mf(i,k) = ud_mf(i,k)
2096 endif
2097 enddo
2099 do k = 1, km
2100 do i = 1, im
2101 if(cnvflg(i).and.rn(i).gt.0.) then
2102 if(k.ge.1 .and. k.le.jmin(i)) then
2103 dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
2104 endif
2105 endif
2106 enddo
2107 enddo
2108
2109 if(mp_phys == mp_phys_mg) then
2110 do k=1,km
2111 do i=1,im
2112 qlcn(i,k) = qlc(i,k)
2113 qicn(i,k) = qli(i,k)
2114 cf_upi(i,k) = cnvc(i,k)
2115 w_upi(i,k) = ud_mf(i,k)*t1(i,k)*rgas /
2116 & (dt2*max(cf_upi(i,k),1.e-12)*prslp(i,k))
2117 cnv_mfd(i,k) = ud_mf(i,k)/dt2
2118 clcn(i,k) = cnvc(i,k)
2119 cnv_fice(i,k) = qicn(i,k)
2120 & / max(1.e-10,qlcn(i,k)+qicn(i,k))
2121 enddo
2122 enddo
2123 endif
2124
2125!!
2126 return
2129 end subroutine sascnvn_run
2130
2131 end module sascnvn
2132! \section original Original Documentation
2133! Penetrative convection is simulated following Pan and Wu (1994), which is based on Arakawa and Schubert(1974) as simplified by Grell (1993) and with a saturated downdraft. Convection occurs when the cloud work function (CWF) exceeds a certain threshold. Mass flux of the cloud is determined using a quasi-equilibrium assumption based on this threshold CWF. The CWF is a function of temperature and moisture in each air column of the model gridpoint. The temperature and moisture profiles are adjusted towards the equilibrium CWF within a specified time scale using the deduced mass flux. A major simplification of the original Arakawa-Shubert scheme is to consider only the deepest cloud and not the spectrum of clouds. The cloud model incorporates a downdraft mechanism as well as the evaporation of precipitation. Entrainment of the updraft and detrainment of the downdraft in the sub-cloud layers are included. Downdraft strength is based on the vertical wind shear through the cloud. The critical CWF is a function of the cloud base vertical motion. As the large-scale rising motion becomes strong, the CWF [similar to convective available potential energy (CAPE)] is allowed to approach zero (therefore approaching neutral stability).
2134!
2135! Mass fluxes induced in the updraft and the downdraft are allowed to transport momentum. The momentum exchange is calculated through the mass flux formulation in a manner similar to that for heat and moisture. The effect of the convection-induced pressure gradient force on cumulus momentum transport is parameterized in terms of mass flux and vertical wind shear (Han and Pan, 2006). As a result, the cumulus momentum exchange is reduced by about 55 % compared to the full exchange.
2136!
2137! The entrainment rate in cloud layers is dependent upon environmental humidity (Han and Pan, 2010). A drier environment increases the entrainment, suppressing the convection. The entrainment rate in sub-cloud layers is given as inversely proportional to height. The detrainment rate is assumed to be a constant in all layers and equal to the entrainment rate value at cloud base, which is O(10-4). The liquid water in the updraft layer is assumed to be detrained from the layers above the level of the minimum moist static energy into the grid-scale cloud water with conversion parameter of 0.002 m-1, which is same as the rain conversion parameter.
2138!
2139! Following Han and Pan (2010), the trigger condition is that a parcel lifted from the convection starting level without entrainment must reach its level of free convection within 120-180 hPa of ascent, proportional to the large-scale vertical velocity. This is intended to produce more convection in large-scale convergent regions but less convection in large-scale subsidence regions. Another important trigger mechanism is to include the effect of environmental humidity in the sub-cloud layer, taking into account convection inhibition due to existence of dry layers below cloud base. On the other hand, the cloud parcel might overshoot beyond the level of neutral buoyancy due to its inertia, eventually stopping its overshoot at cloud top. The CWF is used to model the overshoot. The overshoot of the cloud top is stopped at the height where a parcel lifted from the neutral buoyancy level with energy equal to 10% of the CWF would first have zero energy.
2140!
2141! Deep convection parameterization (SAS) modifications include:
2142! - Detraining cloud water from every updraft layer
2143! - Starting convection from the level of maximum moist static energy within PBL
2144! - Random cloud top is eliminated and only deepest cloud is considered
2145! - Cloud water is detrained from every cloud layer
2146! - Finite entrainment and detrainment rates for heat, moisture, and momentum are specified
2147! - Similar to shallow convection scheme,
2148! - entrainment rate is given to be inversely proportional to height in sub-cloud layers
2149! - detrainment rate is set to be a constant as entrainment rate at the cloud base.
2150! -Above cloud base, an organized entrainment is added, which is a function of environmental relative humidity
subroutine, public sascnvn_run(
This subroutine contains the entirety of the SAS deep convection scheme.
Definition sascnvn.F:93
subroutine, public sascnvn_init(imfdeepcnv, imfdeepcnv_sas, errmsg, errflg)
Definition sascnvn.F:29
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26