CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
shalcnv.F
1
3
10 module shalcnv
11
12 implicit none
13
14 private
15
16 public :: shalcnv_init, shalcnv_run
17
18 contains
19
20!!
21!! \section arg_table_shalcnv_init Argument Table
22!! \htmlinclude shalcnv_init.html
23!!
24 subroutine shalcnv_init(do_shoc,shal_cnv,imfshalcnv, &
25 & imfshalcnv_sas,errmsg,errflg)
26!
27 logical, intent(in) :: do_shoc,shal_cnv
28 integer, intent(in) :: imfshalcnv, imfshalcnv_sas
29 character(len=*), intent(out) :: errmsg
30 integer, intent(out) :: errflg
31!
32! initialize CCPP error handling variables
33 errmsg = ''
34 errflg = 0
35!
36 if (do_shoc .or. .not.shal_cnv .or. &
37 & imfshalcnv/=imfshalcnv_sas) then
38 write(errmsg,'(*(a))') 'Logic error: shalcnv incompatible with',&
39 & ' control flags do_shoc, shal_cnv or imfshalcnv'
40 errflg = 1
41 return
42 endif
43!
44 end subroutine shalcnv_init
45
90 subroutine shalcnv_run( &
91 & grav,cp,hvap,rv,fv,t0c,rd,cvap,cliq,eps,epsm1, &
92 & im,km,jcap,delt,delp,prslp,psp,phil,qlc,qli, &
93 & q1,t1,u1,v1,rn,kbot,ktop,kcnv,islimsk, &
94 & dot,ncloud,hpbl,heat,evap,ud_mf,dt_mf,cnvw,cnvc, &
95 & clam,c0,c1,pgcon,errmsg,errflg)
96!
97 use machine , only : kind_phys
98 use funcphys , only : fpvs
99! use physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
100! &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
101! &, rd => con_rd, cvap => con_cvap, cliq => con_cliq &
102! &, eps => con_eps, epsm1 => con_epsm1
103 implicit none
104!
105! Interface variables
106!
107 real(kind=kind_phys), intent(in) :: grav, cp, hvap, rv, fv, t0c, &
108 & rd, cvap, cliq, eps, epsm1
109 integer, intent(in) :: im, km, jcap, ncloud
110 integer, intent(inout) :: kbot(:), ktop(:), kcnv(:)
111 integer, intent(in) :: islimsk(:)
112 real(kind=kind_phys), intent(in) :: delt, clam, c0, c1, pgcon
113 real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), &
114 & prslp(:,:), dot(:,:), &
115 & phil(:,:), hpbl(:), &
116 & heat(:), evap(:)
117 real(kind=kind_phys), intent(inout) :: &
118 & qlc(:,:), qli(:,:), &
119 & q1(:,:), t1(:,:), &
120 & u1(:,:), v1(:,:), &
121 & cnvw(:,:), cnvc(:,:)
122 real(kind=kind_phys), intent(out) :: rn(:), dt_mf(:,:)
123 real(kind=kind_phys), intent(out), optional :: ud_mf(:,:)
124 character(len=*), intent(out) :: errmsg
125 integer, intent(out) :: errflg
126!
127! Local variables
128!
129 integer i,j,indx, k, kk, km1
130 integer kpbl(im)
131!
132 real(kind=kind_phys) dellat, delta,
133 & desdt,
134 & dp,
135 & dq, dqsdp, dqsdt, dt,
136 & dt2, dtmax, dtmin, dv1h,
137 & dv1q, dv2h, dv2q, dv1u,
138 & dv1v, dv2u, dv2v, dv3q,
139 & dv3h, dv3u, dv3v,
140 & dz, dz1, e1,
141 & el2orc, elocp, aafac,
142 & es, etah, h1, dthk,
143 & evef, evfact, evfactl, fact1,
144 & fact2, factor, fjcap,
145 & g, gamma, pprime, betaw,
146 & qlk, qrch, qs,
147 & rfact, shear, tem1,
148 & val, val1,
149 & val2, w1, w1l, w1s,
150 & w2, w2l, w2s, w3,
151 & w3l, w3s, w4, w4l,
152 & w4s, tem, ptem, ptem1
153!
154 integer kb(im), kbcon(im), kbcon1(im),
155 & ktcon(im), ktcon1(im),
156 & kbm(im), kmax(im)
157!
158 real(kind=kind_phys) aa1(im),
159 & delhbar(im), delq(im), delq2(im),
160 & delqbar(im), delqev(im), deltbar(im),
161 & deltv(im), edt(im),
162 & wstar(im), sflx(im),
163 & pdot(im), po(im,km),
164 & qcond(im), qevap(im), hmax(im),
165 & rntot(im), vshear(im),
166 & xlamud(im), xmb(im), xmbmax(im),
167 & delubar(im), delvbar(im),
168 & ps(im), del(im,km), prsl(im,km)
169!
170 real(kind=kind_phys) cincr, cincrmax, cincrmin
171!
172! physical parameters
173! parameter(g=grav)
174! parameter(elocp=hvap/cp,
175! & el2orc=hvap*hvap/(rv*cp))
176! parameter(c0=.002,c1=5.e-4,delta=fv)
177! parameter(delta=fv)
178! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
179 parameter(cincrmax=180.,cincrmin=120.,dthk=25.)
180 parameter(h1=0.33333333)
181! local variables and arrays
182 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
183 & uo(im,km), vo(im,km), qeso(im,km)
184! cloud water
185! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
186 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km),
187 & dbyo(im,km), zo(im,km), xlamue(im,km),
188 & heo(im,km), heso(im,km),
189 & dellah(im,km), dellaq(im,km),
190 & dellau(im,km), dellav(im,km), hcko(im,km),
191 & ucko(im,km), vcko(im,km), qcko(im,km),
192 & qrcko(im,km), eta(im,km),
193 & zi(im,km), pwo(im,km),
194 & tx1(im), cnvwt(im,km)
195!
196 logical totflg, cnvflg(im), flg(im)
197!
198 real(kind=kind_phys) tf, tcr, tcrf
199 parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
200!
201!-----------------------------------------------------------------------
202!
203!************************************************************************
204! replace (derived) constants above with regular variables
205 g = grav
206 elocp = hvap/cp
207 el2orc = hvap*hvap/(rv*cp)
208 delta = fv
209 fact1 = (cvap-cliq)/rv
210 fact2 = hvap/rv-fact1*t0c
211!************************************************************************
212! initialize CCPP error handling variables
213 errmsg = ''
214 errflg = 0
215!************************************************************************
216! convert input pa terms to cb terms -- moorthi
219 ps = psp * 0.001
220 prsl = prslp * 0.001
221 del = delp * 0.001
222!************************************************************************
223!
224 km1 = km - 1
225!
226! compute surface buoyancy flux
227!
233 do i=1,im
234 sflx(i) = heat(i)+fv*t1(i,1)*evap(i)
235 enddo
236!
237! initialize arrays
238!
240 do i=1,im
241 cnvflg(i) = .true.
242 if(kcnv(i).eq.1) cnvflg(i) = .false.
243 if(sflx(i).le.0.) cnvflg(i) = .false.
244 if(cnvflg(i)) then
245 kbot(i)=km+1
246 ktop(i)=0
247 endif
248 rn(i)=0.
249 kbcon(i)=km
250 ktcon(i)=1
251 kb(i)=km
252 pdot(i) = 0.
253 qlko_ktcon(i) = 0.
254 edt(i) = 0.
255 aa1(i) = 0.
256 vshear(i) = 0.
257 enddo
259! hchuang code change
260 do k = 1, km
261 do i = 1, im
262 ud_mf(i,k) = 0.
263 dt_mf(i,k) = 0.
264 enddo
265 enddo
266!!
268 totflg = .true.
269 do i=1,im
270 totflg = totflg .and. (.not. cnvflg(i))
271 enddo
272 if(totflg) return
273!!
274!
276 dt2 = delt
277 val = 1200.
278 dtmin = max(dt2, val )
279 val = 3600.
280 dtmax = max(dt2, val )
281! model tunable parameters are all here
282! clam = .3
283 aafac = .1
284 betaw = .03
285! evef = 0.07
286 evfact = 0.3
287 evfactl = 0.3
288!
289! pgcon = 0.7 ! gregory et al. (1997, qjrms)
290! pgcon = 0.55 ! zhang & wu (2003,jas)
291 fjcap = (float(jcap) / 126.) ** 2
292 val = 1.
293 fjcap = max(fjcap,val)
294 w1l = -8.e-3
295 w2l = -4.e-2
296 w3l = -5.e-3
297 w4l = -5.e-4
298 w1s = -2.e-4
299 w2s = -2.e-3
300 w3s = -1.e-3
301 w4s = -2.e-5
302!
303! define top layer for search of the downdraft originating layer
304! and the maximum thetae for updraft
305!
307 do i=1,im
308 kbm(i) = km
309 kmax(i) = km
310 tx1(i) = 1.0 / ps(i)
311 enddo
312!
313 do k = 1, km
314 do i=1,im
315 if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
316 if (prsl(i,k)*tx1(i) .gt. 0.60) kmax(i) = k + 1
317 enddo
318 enddo
319 do i=1,im
320 kbm(i) = min(kbm(i),kmax(i))
321 enddo
322!
323! hydrostatic height assume zero terr and compute
324! updraft entrainment rate as an inverse function of height
325!
327 do k = 1, km
328 do i=1,im
329 zo(i,k) = phil(i,k) / g
330 enddo
331 enddo
333 do k = 1, km1
334 do i=1,im
335 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
336 xlamue(i,k) = clam / zi(i,k)
337 enddo
338 enddo
339 do i=1,im
340 xlamue(i,km) = xlamue(i,km1)
341 enddo
342!
343! pbl height
344!
346 do i=1,im
347 flg(i) = cnvflg(i)
348 kpbl(i)= 1
349 enddo
350 do k = 2, km1
351 do i=1,im
352 if (flg(i).and.zo(i,k).le.hpbl(i)) then
353 kpbl(i) = k
354 else
355 flg(i) = .false.
356 endif
357 enddo
358 enddo
359 do i=1,im
360 kpbl(i)= min(kpbl(i),kbm(i))
361 enddo
362!
363!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
364! convert surface pressure to mb from cb
365!
367 do k = 1, km
368 do i = 1, im
369 if (cnvflg(i) .and. k .le. kmax(i)) then
370 pfld(i,k) = prsl(i,k) * 10.0
371 eta(i,k) = 1.
372 hcko(i,k) = 0.
373 qcko(i,k) = 0.
374 qrcko(i,k)= 0.
375 ucko(i,k) = 0.
376 vcko(i,k) = 0.
377 dbyo(i,k) = 0.
378 pwo(i,k) = 0.
379 dellal(i,k) = 0.
380 to(i,k) = t1(i,k)
381 qo(i,k) = q1(i,k)
382 uo(i,k) = u1(i,k)
383 vo(i,k) = v1(i,k)
384! uo(i,k) = u1(i,k) * rcs(i)
385! vo(i,k) = v1(i,k) * rcs(i)
386 cnvwt(i,k) = 0.
387 endif
388 enddo
389 enddo
390!
391! column variables
392! p is pressure of the layer (mb)
393! t is temperature at t-dt (k)..tn
394! q is mixing ratio at t-dt (kg/kg)..qn
395! to is temperature at t+dt (k)... this is after advection and turbulan
396! qo is mixing ratio at t+dt (kg/kg)..q1
397!
399 do k = 1, km
400 do i=1,im
401 if (cnvflg(i) .and. k .le. kmax(i)) then
402 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
403 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
404 val1 = 1.e-8
405 qeso(i,k) = max(qeso(i,k), val1)
406 val2 = 1.e-10
407 qo(i,k) = max(qo(i,k), val2 )
408! qo(i,k) = min(qo(i,k),qeso(i,k))
409! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
410 endif
411 enddo
412 enddo
413!
414! compute moist static energy
415!
417 do k = 1, km
418 do i=1,im
419 if (cnvflg(i) .and. k .le. kmax(i)) then
420! tem = g * zo(i,k) + cp * to(i,k)
421 tem = phil(i,k) + cp * to(i,k)
422 heo(i,k) = tem + hvap * qo(i,k)
423 heso(i,k) = tem + hvap * qeso(i,k)
424! heo(i,k) = min(heo(i,k),heso(i,k))
425 endif
426 enddo
427 enddo
428!
429! determine level with largest moist static energy within pbl
430! this is the level where updraft starts
431!
434 do i=1,im
435 if (cnvflg(i)) then
436 hmax(i) = heo(i,1)
437 kb(i) = 1
438 endif
439 enddo
440 do k = 2, km
441 do i=1,im
442 if (cnvflg(i).and.k.le.kpbl(i)) then
443 if(heo(i,k).gt.hmax(i)) then
444 kb(i) = k
445 hmax(i) = heo(i,k)
446 endif
447 endif
448 enddo
449 enddo
450!
452 do k = 1, km1
453 do i=1,im
454 if (cnvflg(i) .and. k .le. kmax(i)-1) then
455 dz = .5 * (zo(i,k+1) - zo(i,k))
456 dp = .5 * (pfld(i,k+1) - pfld(i,k))
457 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
458 pprime = pfld(i,k+1) + epsm1 * es
459 qs = eps * es / pprime
460 dqsdp = - qs / pprime
461 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
462 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
463 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
464 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
465 dq = dqsdt * dt + dqsdp * dp
466 to(i,k) = to(i,k+1) + dt
467 qo(i,k) = qo(i,k+1) + dq
468 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
469 endif
470 enddo
471 enddo
472!
474 do k = 1, km1
475 do i=1,im
476 if (cnvflg(i) .and. k .le. kmax(i)-1) then
477 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
478 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
479 val1 = 1.e-8
480 qeso(i,k) = max(qeso(i,k), val1)
481 val2 = 1.e-10
482 qo(i,k) = max(qo(i,k), val2 )
483! qo(i,k) = min(qo(i,k),qeso(i,k))
484 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
485 & cp * to(i,k) + hvap * qo(i,k)
486 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
487 & cp * to(i,k) + hvap * qeso(i,k)
488 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
489 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
490 endif
491 enddo
492 enddo
493!
494! look for the level of free convection as cloud base
495!!> - Search below the index "kbm" for the level of free convection (LFC) where the condition \f$h_b > h^*\f$ is first met, where \f$h_b, h^*\f$ are the state moist static energy at the parcel's starting level and saturation moist static energy, respectively. Set "kbcon" to the index of the LFC.
496 do i=1,im
497 flg(i) = cnvflg(i)
498 if(flg(i)) kbcon(i) = kmax(i)
499 enddo
500 do k = 2, km1
501 do i=1,im
502 if (flg(i).and.k.lt.kbm(i)) then
503 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
504 kbcon(i) = k
505 flg(i) = .false.
506 endif
507 endif
508 enddo
509 enddo
510!
512 do i=1,im
513 if(cnvflg(i)) then
514 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
515 endif
516 enddo
517!!
518 totflg = .true.
519 do i=1,im
520 totflg = totflg .and. (.not. cnvflg(i))
521 enddo
522 if(totflg) return
523!!
524!
525! determine critical convective inhibition
526! as a function of vertical velocity at cloud base.
527!
529 do i=1,im
530 if(cnvflg(i)) then
531! pdot(i) = 10.* dot(i,kbcon(i))
532 pdot(i) = 0.01 * dot(i,kbcon(i)) ! now dot is in pa/s
533 endif
534 enddo
535 do i=1,im
536 if(cnvflg(i)) then
537 if(islimsk(i) == 1) then
538 w1 = w1l
539 w2 = w2l
540 w3 = w3l
541 w4 = w4l
542 else
543 w1 = w1s
544 w2 = w2s
545 w3 = w3s
546 w4 = w4s
547 endif
548 if(pdot(i).le.w4) then
549 ptem = (pdot(i) - w4) / (w3 - w4)
550 elseif(pdot(i).ge.-w4) then
551 ptem = - (pdot(i) + w4) / (w4 - w3)
552 else
553 ptem = 0.
554 endif
555 val1 = -1.
556 ptem = max(ptem,val1)
557 val2 = 1.
558 ptem = min(ptem,val2)
559 ptem = 1. - ptem
560 ptem1= .5*(cincrmax-cincrmin)
561 cincr = cincrmax - ptem * ptem1
562 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
563 if(tem1.gt.cincr) then
564 cnvflg(i) = .false.
565 endif
566 endif
567 enddo
568!!
569 totflg = .true.
570 do i=1,im
571 totflg = totflg .and. (.not. cnvflg(i))
572 enddo
573 if(totflg) return
574!!
575!
576! assume the detrainment rate for the updrafts to be same as
577! the entrainment rate at cloud base
578!
580 do i = 1, im
581 if(cnvflg(i)) then
582 xlamud(i) = xlamue(i,kbcon(i))
583 endif
584 enddo
585!
586! determine updraft mass flux for the subcloud layers
587!
593 do k = km1, 1, -1
594 do i = 1, im
595 if (cnvflg(i)) then
596 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
597 dz = zi(i,k+1) - zi(i,k)
598 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
599 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
600 endif
601 endif
602 enddo
603 enddo
604!
605! compute mass flux above cloud base
606!
607 do k = 2, km1
608 do i = 1, im
609 if(cnvflg(i))then
610 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
611 dz = zi(i,k) - zi(i,k-1)
612 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
613 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
614 endif
615 endif
616 enddo
617 enddo
618!
619! compute updraft cloud property
620!
622 do i = 1, im
623 if(cnvflg(i)) then
624 indx = kb(i)
625 hcko(i,indx) = heo(i,indx)
626 ucko(i,indx) = uo(i,indx)
627 vcko(i,indx) = vo(i,indx)
628 endif
629 enddo
630!
632 do k = 2, km1
633 do i = 1, im
634 if (cnvflg(i)) then
635 if(k.gt.kb(i).and.k.lt.kmax(i)) then
636 dz = zi(i,k) - zi(i,k-1)
637 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
638 tem1 = 0.5 * xlamud(i) * dz
639 factor = 1. + tem - tem1
640 ptem = 0.5 * tem + pgcon
641 ptem1= 0.5 * tem - pgcon
642 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
643 & (heo(i,k)+heo(i,k-1)))/factor
644 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)
645 & +ptem1*uo(i,k-1))/factor
646 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)
647 & +ptem1*vo(i,k-1))/factor
648 dbyo(i,k) = hcko(i,k) - heso(i,k)
649 endif
650 endif
651 enddo
652 enddo
653!
654! taking account into convection inhibition due to existence of
655! dry layers below cloud base
656!
658 do i=1,im
659 flg(i) = cnvflg(i)
660 kbcon1(i) = kmax(i)
661 enddo
662 do k = 2, km1
663 do i=1,im
664 if (flg(i).and.k.lt.kbm(i)) then
665 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
666 kbcon1(i) = k
667 flg(i) = .false.
668 endif
669 endif
670 enddo
671 enddo
672 do i=1,im
673 if(cnvflg(i)) then
674 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
675 endif
676 enddo
677 do i=1,im
678 if(cnvflg(i)) then
679 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
680 if(tem.gt.dthk) then
681 cnvflg(i) = .false.
682 endif
683 endif
684 enddo
685!!
686 totflg = .true.
687 do i = 1, im
688 totflg = totflg .and. (.not. cnvflg(i))
689 enddo
690 if(totflg) return
691!!
692!
693! determine first guess cloud top as the level of zero buoyancy
694! limited to the level of sigma=0.7
695!
697 do i = 1, im
698 flg(i) = cnvflg(i)
699 if(flg(i)) ktcon(i) = kbm(i)
700 enddo
701 do k = 2, km1
702 do i=1,im
703 if (flg(i).and.k .lt. kbm(i)) then
704 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
705 ktcon(i) = k
706 flg(i) = .false.
707 endif
708 endif
709 enddo
710 enddo
711!
712! turn off shallow convection if cloud top is less than pbl top
713!
714! do i=1,im
715! if(cnvflg(i)) then
716! kk = kpbl(i)+1
717! if(ktcon(i).le.kk) cnvflg(i) = .false.
718! endif
719! enddo
720!!
721! totflg = .true.
722! do i = 1, im
723! totflg = totflg .and. (.not. cnvflg(i))
724! enddo
725! if(totflg) return
726!!
727!
728! specify upper limit of mass flux at cloud base
729!
731 do i = 1, im
732 if(cnvflg(i)) then
733! xmbmax(i) = .1
734!
735 k = kbcon(i)
736 dp = 1000. * del(i,k)
737 xmbmax(i) = dp / (g * dt2)
738!
739! tem = dp / (g * dt2)
740! xmbmax(i) = min(tem, xmbmax(i))
741 endif
742 enddo
743!
744! compute cloud moisture property and precipitation
745!
747 do i = 1, im
748 if (cnvflg(i)) then
749 aa1(i) = 0.
750 qcko(i,kb(i)) = qo(i,kb(i))
751 qrcko(i,kb(i)) = qo(i,kb(i))
752 endif
753 enddo
755 do k = 2, km1
756 do i = 1, im
757 if (cnvflg(i)) then
758 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
759 dz = zi(i,k) - zi(i,k-1)
760 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
761 qrch = qeso(i,k)
762 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
763!j
764 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
765 tem1 = 0.5 * xlamud(i) * dz
766 factor = 1. + tem - tem1
767 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
768 & (qo(i,k)+qo(i,k-1)))/factor
769 qrcko(i,k) = qcko(i,k)
770!j
771 dq = eta(i,k) * (qcko(i,k) - qrch)
772!
773! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
774!
775! below lfc check if there is excess moisture to release latent heat
776!
777 if(k.ge.kbcon(i).and.dq.gt.0.) then
778 etah = .5 * (eta(i,k) + eta(i,k-1))
779 if(ncloud.gt.0.) then
780 dp = 1000. * del(i,k)
781 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
782 dellal(i,k) = etah * c1 * dz * qlk * g / dp
783 else
784 qlk = dq / (eta(i,k) + etah * c0 * dz)
785 endif
786 aa1(i) = aa1(i) - dz * g * qlk
787 qcko(i,k)= qlk + qrch
788 pwo(i,k) = etah * c0 * dz * qlk
789 cnvwt(i,k) = etah * qlk * g / dp
790 endif
791 endif
792 endif
793 enddo
794 enddo
795!
796! calculate cloud work function
797!
803 do k = 2, km1
804 do i = 1, im
805 if (cnvflg(i)) then
806 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
807 dz1 = zo(i,k+1) - zo(i,k)
808 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
809 rfact = 1. + delta * cp * gamma
810 & * to(i,k) / hvap
811 aa1(i) = aa1(i) +
812 & dz1 * (g / (cp * to(i,k)))
813 & * dbyo(i,k) / (1. + gamma)
814 & * rfact
815 val = 0.
816 aa1(i)=aa1(i)+
817 & dz1 * g * delta *
818 & max(val,(qeso(i,k) - qo(i,k)))
819 endif
820 endif
821 enddo
822 enddo
824 do i = 1, im
825 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
826 enddo
827!!
828 totflg = .true.
829 do i=1,im
830 totflg = totflg .and. (.not. cnvflg(i))
831 enddo
832 if(totflg) return
833!!
834!
835! estimate the onvective overshooting as the level
836! where the [aafac * cloud work function] becomes zero,
837! which is the final cloud top
838! limited to the level of sigma=0.7
839!
841 do i = 1, im
842 if (cnvflg(i)) then
843 aa1(i) = aafac * aa1(i)
844 endif
845 enddo
846!
847 do i = 1, im
848 flg(i) = cnvflg(i)
849 ktcon1(i) = kbm(i)
850 enddo
851 do k = 2, km1
852 do i = 1, im
853 if (flg(i)) then
854 if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
855 dz1 = zo(i,k+1) - zo(i,k)
856 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
857 rfact = 1. + delta * cp * gamma
858 & * to(i,k) / hvap
859 aa1(i) = aa1(i) +
860 & dz1 * (g / (cp * to(i,k)))
861 & * dbyo(i,k) / (1. + gamma)
862 & * rfact
863 if(aa1(i).lt.0.) then
864 ktcon1(i) = k
865 flg(i) = .false.
866 endif
867 endif
868 endif
869 enddo
870 enddo
871!
872! compute cloud moisture property, detraining cloud water
873! and precipitation in overshooting layers
874!
876 do k = 2, km1
877 do i = 1, im
878 if (cnvflg(i)) then
879 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
880 dz = zi(i,k) - zi(i,k-1)
881 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
882 qrch = qeso(i,k)
883 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
884!j
885 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
886 tem1 = 0.5 * xlamud(i) * dz
887 factor = 1. + tem - tem1
888 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
889 & (qo(i,k)+qo(i,k-1)))/factor
890 qrcko(i,k) = qcko(i,k)
891!j
892 dq = eta(i,k) * (qcko(i,k) - qrch)
893!
894! check if there is excess moisture to release latent heat
895!
896 if(dq.gt.0.) then
897 etah = .5 * (eta(i,k) + eta(i,k-1))
898 if(ncloud.gt.0.) then
899 dp = 1000. * del(i,k)
900 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
901 dellal(i,k) = etah * c1 * dz * qlk * g / dp
902 else
903 qlk = dq / (eta(i,k) + etah * c0 * dz)
904 endif
905 qcko(i,k) = qlk + qrch
906 pwo(i,k) = etah * c0 * dz * qlk
907 cnvwt(i,k) = etah * qlk * g / dp
908 endif
909 endif
910 endif
911 enddo
912 enddo
913!
914! exchange ktcon with ktcon1
915!
916 do i = 1, im
917 if(cnvflg(i)) then
918 kk = ktcon(i)
919 ktcon(i) = ktcon1(i)
920 ktcon1(i) = kk
921 endif
922 enddo
923!
924! this section is ready for cloud water
925!
926 if(ncloud.gt.0) then
927!
928! compute liquid and vapor separation at cloud top
929!
931 do i = 1, im
932 if(cnvflg(i)) then
933 k = ktcon(i) - 1
934 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
935 qrch = qeso(i,k)
936 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
937 dq = qcko(i,k) - qrch
938!
939! check if there is excess moisture to release latent heat
940!
941 if(dq.gt.0.) then
942 qlko_ktcon(i) = dq
943 qcko(i,k) = qrch
944 endif
945 endif
946 enddo
947 endif
948!
949!--- compute precipitation efficiency in terms of windshear
950!
951!! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 :
952!! \f[
953!! E = 1.591 - 0.639\frac{\Delta V}{\Delta z} + 0.0953\left(\frac{\Delta V}{\Delta z}\right)^2 - 0.00496\left(\frac{\Delta V}{\Delta z}\right)^3
954!! \f]
955!! where \f$\Delta V\f$ is the integrated horizontal shear over the cloud depth, \f$\Delta z\f$, (the ratio is converted to units of \f$10^{-3} s^{-1}\f$). The variable "edt" is \f$1-E\f$ and is constrained to the range \f$[0,0.9]\f$.
956 do i = 1, im
957 if(cnvflg(i)) then
958 vshear(i) = 0.
959 endif
960 enddo
961 do k = 2, km
962 do i = 1, im
963 if (cnvflg(i)) then
964 if(k.gt.kb(i).and.k.le.ktcon(i)) then
965 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
966 & + (vo(i,k)-vo(i,k-1)) ** 2)
967 vshear(i) = vshear(i) + shear
968 endif
969 endif
970 enddo
971 enddo
972 do i = 1, im
973 if(cnvflg(i)) then
974 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
975 e1=1.591-.639*vshear(i)
976 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
977 edt(i)=1.-e1
978 val = .9
979 edt(i) = min(edt(i),val)
980 val = .0
981 edt(i) = max(edt(i),val)
982 endif
983 enddo
984!
985!--- what would the change be, that a cloud with unit mass
986!--- will do to the environment?
987!
990 do k = 1, km
991 do i = 1, im
992 if(cnvflg(i) .and. k .le. kmax(i)) then
993 dellah(i,k) = 0.
994 dellaq(i,k) = 0.
995 dellau(i,k) = 0.
996 dellav(i,k) = 0.
997 endif
998 enddo
999 enddo
1000!
1001!--- changed due to subsidence and entrainment
1002!
1003 do k = 2, km1
1004 do i = 1, im
1005 if (cnvflg(i)) then
1006 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
1007 dp = 1000. * del(i,k)
1008 dz = zi(i,k) - zi(i,k-1)
1009!
1010 dv1h = heo(i,k)
1011 dv2h = .5 * (heo(i,k) + heo(i,k-1))
1012 dv3h = heo(i,k-1)
1013 dv1q = qo(i,k)
1014 dv2q = .5 * (qo(i,k) + qo(i,k-1))
1015 dv3q = qo(i,k-1)
1016 dv1u = uo(i,k)
1017 dv2u = .5 * (uo(i,k) + uo(i,k-1))
1018 dv3u = uo(i,k-1)
1019 dv1v = vo(i,k)
1020 dv2v = .5 * (vo(i,k) + vo(i,k-1))
1021 dv3v = vo(i,k-1)
1022!
1023 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1024 tem1 = xlamud(i)
1025!j
1026 dellah(i,k) = dellah(i,k) +
1027 & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h
1028 & - tem*eta(i,k-1)*dv2h*dz
1029 & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1030 & ) *g/dp
1031!j
1032 dellaq(i,k) = dellaq(i,k) +
1033 & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q
1034 & - tem*eta(i,k-1)*dv2q*dz
1035 & + tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz
1036 & ) *g/dp
1037!j
1038 dellau(i,k) = dellau(i,k) +
1039 & ( eta(i,k)*dv1u - eta(i,k-1)*dv3u
1040 & - tem*eta(i,k-1)*dv2u*dz
1041 & + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz
1042 & - pgcon*eta(i,k-1)*(dv1u-dv3u)
1043 & ) *g/dp
1044!j
1045 dellav(i,k) = dellav(i,k) +
1046 & ( eta(i,k)*dv1v - eta(i,k-1)*dv3v
1047 & - tem*eta(i,k-1)*dv2v*dz
1048 & + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz
1049 & - pgcon*eta(i,k-1)*(dv1v-dv3v)
1050 & ) *g/dp
1051!j
1052 endif
1053 endif
1054 enddo
1055 enddo
1056!
1057!------- cloud top
1058!
1059 do i = 1, im
1060 if(cnvflg(i)) then
1061 indx = ktcon(i)
1062 dp = 1000. * del(i,indx)
1063 dv1h = heo(i,indx-1)
1064 dellah(i,indx) = eta(i,indx-1) *
1065 & (hcko(i,indx-1) - dv1h) * g / dp
1066 dv1q = qo(i,indx-1)
1067 dellaq(i,indx) = eta(i,indx-1) *
1068 & (qcko(i,indx-1) - dv1q) * g / dp
1069 dv1u = uo(i,indx-1)
1070 dellau(i,indx) = eta(i,indx-1) *
1071 & (ucko(i,indx-1) - dv1u) * g / dp
1072 dv1v = vo(i,indx-1)
1073 dellav(i,indx) = eta(i,indx-1) *
1074 & (vcko(i,indx-1) - dv1v) * g / dp
1075!
1076! cloud water
1077!
1078 dellal(i,indx) = eta(i,indx-1) *
1079 & qlko_ktcon(i) * g / dp
1080 endif
1081 enddo
1082!
1083! mass flux at cloud base for shallow convection
1084! (grant, 2001)
1085!
1091 do i= 1, im
1092 if(cnvflg(i)) then
1093 k = kbcon(i)
1094! ptem = g*sflx(i)*zi(i,k)/t1(i,1)
1095 ptem = g*sflx(i)*hpbl(i)/t1(i,1)
1096 wstar(i) = ptem**h1
1097 tem = po(i,k)*100. / (rd*t1(i,k))
1098 xmb(i) = betaw*tem*wstar(i)
1099 xmb(i) = min(xmb(i),xmbmax(i))
1100 endif
1101 enddo
1104!
1105!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1106!
1107 do k = 1, km
1108 do i = 1, im
1109 if (cnvflg(i) .and. k .le. kmax(i)) then
1110 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1111 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
1112 val = 1.e-8
1113 qeso(i,k) = max(qeso(i,k), val )
1114 endif
1115 enddo
1116 enddo
1117!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1118!
1122 do i = 1, im
1123 delhbar(i) = 0.
1124 delqbar(i) = 0.
1125 deltbar(i) = 0.
1126 delubar(i) = 0.
1127 delvbar(i) = 0.
1128 qcond(i) = 0.
1129 enddo
1130 do k = 1, km
1131 do i = 1, im
1132 if (cnvflg(i)) then
1133 if(k.gt.kb(i).and.k.le.ktcon(i)) then
1134 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1135 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1136 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1137! tem = 1./rcs(i)
1138! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1139! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1140 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
1141 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
1142 dp = 1000. * del(i,k)
1143 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
1144 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
1145 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
1146 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
1147 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
1148 endif
1149 endif
1150 enddo
1151 enddo
1153 do k = 1, km
1154 do i = 1, im
1155 if (cnvflg(i)) then
1156 if(k.gt.kb(i).and.k.le.ktcon(i)) then
1157 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
1158 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
1159 val = 1.e-8
1160 qeso(i,k) = max(qeso(i,k), val )
1161 endif
1162 endif
1163 enddo
1164 enddo
1165!
1167 do i = 1, im
1168 rntot(i) = 0.
1169 delqev(i) = 0.
1170 delq2(i) = 0.
1171 flg(i) = cnvflg(i)
1172 enddo
1173 do k = km, 1, -1
1174 do i = 1, im
1175 if (cnvflg(i)) then
1176 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
1177 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
1178 endif
1179 endif
1180 enddo
1181 enddo
1182!
1183! evaporating rain
1184!
1188 do k = km, 1, -1
1189 do i = 1, im
1190 if (k .le. kmax(i)) then
1191 deltv(i) = 0.
1192 delq(i) = 0.
1193 qevap(i) = 0.
1194 if(cnvflg(i)) then
1195 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
1196 rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
1197 endif
1198 endif
1199 if(flg(i).and.k.lt.ktcon(i)) then
1200 evef = edt(i) * evfact
1201 if(islimsk(i) == 1) evef=edt(i) * evfactl
1202! if(islimsk(i) == 1) evef=.07
1203! if(islimsk(i) == 1) evef = 0.
1204 qcond(i) = evef * (q1(i,k) - qeso(i,k))
1205 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
1206 dp = 1000. * del(i,k)
1207 if(rn(i).gt.0..and.qcond(i).lt.0.) then
1208 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
1209 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
1210 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
1211 endif
1212 if(rn(i).gt.0..and.qcond(i).lt.0..and.
1213 & delq2(i).gt.rntot(i)) then
1214 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
1215 flg(i) = .false.
1216 endif
1217 if(rn(i).gt.0..and.qevap(i).gt.0.) then
1218 tem = .001 * dp / g
1219 tem1 = qevap(i) * tem
1220 if(tem1.gt.rn(i)) then
1221 qevap(i) = rn(i) / tem
1222 rn(i) = 0.
1223 else
1224 rn(i) = rn(i) - tem1
1225 endif
1226 q1(i,k) = q1(i,k) + qevap(i)
1227 t1(i,k) = t1(i,k) - elocp * qevap(i)
1228 deltv(i) = - elocp*qevap(i)/dt2
1229 delq(i) = + qevap(i)/dt2
1230 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
1231 endif
1232 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
1233 delqbar(i) = delqbar(i) + delq(i)*dp/g
1234 deltbar(i) = deltbar(i) + deltv(i)*dp/g
1235 endif
1236 endif
1237 enddo
1238 enddo
1239!j
1240! do i = 1, im
1241! if(me.eq.31.and.cnvflg(i)) then
1242! if(cnvflg(i)) then
1243! print *, ' shallow delhbar, delqbar, deltbar = ',
1244! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
1245! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
1246! print *, ' precip =', hvap*rn(i)*1000./dt2
1247! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
1248! endif
1249! enddo
1250!j
1251 do i = 1, im
1252 if(cnvflg(i)) then
1253 if(rn(i).lt.0..or..not.flg(i)) rn(i) = 0.
1254 ktop(i) = ktcon(i)
1255 kbot(i) = kbcon(i)
1256 kcnv(i) = 0
1257 endif
1258 enddo
1259!
1260! convective cloud water
1261!
1263 do k = 1, km
1264 do i = 1, im
1265 if (cnvflg(i) .and. rn(i).gt.0.) then
1266 if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1267 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
1268 endif
1269 endif
1270 enddo
1271 enddo
1272
1273!
1274! convective cloud cover
1275!
1277 do k = 1, km
1278 do i = 1, im
1279 if (cnvflg(i) .and. rn(i).gt.0.) then
1280 if (k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1281 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
1282 cnvc(i,k) = min(cnvc(i,k), 0.2)
1283 cnvc(i,k) = max(cnvc(i,k), 0.0)
1284 endif
1285 endif
1286 enddo
1287 enddo
1288
1289!
1290! cloud water
1291!
1293 if (ncloud.gt.0) then
1294!
1295 do k = 1, km1
1296 do i = 1, im
1297 if (cnvflg(i)) then
1298 if (k.gt.kb(i).and.k.le.ktcon(i)) then
1299 tem = dellal(i,k) * xmb(i) * dt2
1300 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
1301 if (qlc(i,k) .gt. -999.0) then
1302 qli(i,k) = qli(i,k) + tem * tem1 ! ice
1303 qlc(i,k) = qlc(i,k) + tem *(1.0-tem1) ! water
1304 else
1305 qli(i,k) = qli(i,k) + tem
1306 endif
1307 endif
1308 endif
1309 enddo
1310 enddo
1311!
1312 endif
1313!
1314! hchuang code change
1315!
1317 do k = 1, km
1318 do i = 1, im
1319 if(cnvflg(i)) then
1320 if(k.ge.kb(i) .and. k.lt.ktop(i)) then
1321 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
1322 endif
1323 endif
1324 enddo
1325 enddo
1327 do i = 1, im
1328 if(cnvflg(i)) then
1329 k = ktop(i)-1
1330 dt_mf(i,k) = ud_mf(i,k)
1331 endif
1332 enddo
1333!!
1334 return
1335
1336 end subroutine shalcnv_run
1337
1338 end module shalcnv
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
Definition sflx.f:3
The Mass-Flux shallow convection scheme parameterizes the effect of shallow convection on the environ...
Definition shalcnv.F:10