7 !-------------------------------------------------------------------------------
8 subroutine cu_ksas(dt,dx,p3di,p3d,pi3d,qc3d,qi3d,rho3d,itimestep,stepcu, &
9 hbot,htop,cu_act_flag, &
10 rthcuten,rqvcuten,rqccuten,rqicuten, &
12 qv3d,t3d,raincv,pratec,xland,dz8w,w,u3d,v3d, &
15 mp_physics,dx_factor_nsas, &
16 p_qc,p_qi,p_first_scalar, &
18 cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2, &
19 cice,xls,psat,f_qi,f_qc, &
20 ids,ide, jds,jde, kds,kde, &
21 ims,ime, jms,jme, kms,kme, &
22 its,ite, jts,jte, kts,kte)
23 !-------------------------------------------------------------------------------
25 !-------------------------------------------------------------------------------
27 !-- dx grid interval (m)
28 !-- p3di 3d pressure (pa) at interface level
29 !-- p3d 3d pressure (pa)
30 !-- pi3d 3d exner function (dimensionless)
31 !-- z height above sea level (m)
32 !-- qc3d cloud water mixing ratio (kg/kg)
33 !-- qi3d cloud ice mixing ratio (kg/kg)
34 !-- qv3d 3d water vapor mixing ratio (kg/kg)
35 !-- t3d temperature (k)
36 !-- raincv cumulus scheme precipitation (mm)
37 !-- w vertical velocity (m/s)
38 !-- dz8w dz between full levels (m)
39 !-- u3d 3d u-velocity interpolated to theta points (m/s)
40 !-- v3d 3d v-velocity interpolated to theta points (m/s)
41 !-- ids start index for i in domain
42 !-- ide end index for i in domain
43 !-- jds start index for j in domain
44 !-- jde end index for j in domain
45 !-- kds start index for k in domain
46 !-- kde end index for k in domain
47 !-- ims start index for i in memory
48 !-- ime end index for i in memory
49 !-- jms start index for j in memory
50 !-- jme end index for j in memory
51 !-- kms start index for k in memory
52 !-- kme end index for k in memory
53 !-- its start index for i in tile
54 !-- ite end index for i in tile
55 !-- jts start index for j in tile
56 !-- jte end index for j in tile
57 !-- kts start index for k in tile
58 !-- kte end index for k in tile
59 !-------------------------------------------------------------------------------
60 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
61 ims,ime, jms,jme, kms,kme, &
62 its,ite, jts,jte, kts,kte, &
64 p_qc,p_qi,p_first_scalar
65 real, intent(in ) :: cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2, &
67 real, intent(in ) :: dt,dx
68 real, optional, intent(in ) :: pgcon
69 real, dimension( ims:ime, kms:kme, jms:jme ),optional ,&
70 intent(inout) :: rthcuten,&
76 logical, optional :: F_QC,F_QI
77 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
85 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
87 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
90 real, dimension( ims:ime, jms:jme ) ,&
91 intent(inout) :: raincv,&
93 real, dimension( ims:ime, jms:jme ) ,&
97 real, dimension( ims:ime, jms:jme ) ,&
100 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
103 logical, dimension( ims:ime, jms:jme ) ,&
104 intent(inout) :: cu_act_flag
106 real, dimension( ims:ime, jms:jme ) ,&
107 intent(in ) :: hpbl,&
111 real, dimension( kms:kme ) ,&
113 integer, intent(in ) :: mp_physics
114 integer, intent(in ) :: dx_factor_nsas
119 real, dimension( its:ite, jts:jte ) :: raincv1,&
123 real, dimension( its:ite, kts:kte ) :: del,&
132 real, dimension( its:ite, kts:kte+1 ) :: prsii,&
134 real, dimension( its:ite, kts:kte ) :: zll
135 real, dimension( its:ite) :: rain
137 integer, dimension (its:ite) :: kbot,&
141 integer :: i,j,k,kp, kbmax,kbm,kmax
143 ! microphysics scheme --> ncloud
145 if (mp_physics .eq. 0) then
147 elseif ( mp_physics .eq. 1 .or. mp_physics .eq. 3 ) then
153 if(present(pgcon)) then
156 ! pgcon_use = 0.7 ! Gregory et al. (1997, QJRMS)
157 pgcon_use = 0.55 ! Zhang & Wu (2003,JAS)
158 ! 0.55 is a physically-based value used by GFS
159 ! HWRF uses 0.2, for model tuning purposes
164 cu_act_flag(i,j)=.TRUE.
174 if(znu(k).gt.0.45) kbmax = k + 1
175 if(znu(k).gt.0.70) kbm = k + 1
176 if(znu(k).gt.0.05) kmax = k + 1
185 dot(i,k) = -5.0e-4*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
186 prsll(i,k)=p3d(i,k,j)*0.001
187 prsii(i,k)=p3di(i,k,j)*0.001
192 prsii(i,kte+1)=p3di(i,kte+1,j)*0.001
201 zii(i,k+1)=zii(i,k)+dz8w(i,k,j)
207 zll(i,k)=0.5*(zii(i,k)+zii(i,k+1))
213 del(i,k)=prsll(i,k)*g/r_d*dz8w(i,k,j)/t3d(i,k,j)
217 ! q1(i,k)=qv3d(i,k,j)/(1.+qv3d(i,k,j))
219 qi2(i,k) = qi3d(i,k,j)
220 qc2(i,k) = qc3d(i,k,j)
226 call nsas2d(delt=delt,delx=dx,del=del(its,kts), &
227 prsl=prsll(its,kts),prsi=prsii(its,kts),prslk=pi3d(ims,kms,j), &
229 ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
230 q1=q1(its,kts),t1=t1(its,kts),rain=rain(its), &
231 kbot=kbot(its),ktop=ktop(its), &
233 lat=j,slimsk=xland(ims,j),dot=dot(its,kts), &
234 u1=u1(its,kts), v1=v1(its,kts), &
235 cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv, &
236 rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2, &
237 cice=cice,xls=xls,psat=psat, &
238 dx_factor_nsas=dx_factor_nsas, &
239 hpbl=hpbl(ims,j),hpbl_hold=hpbl_hold(ims,j), &
240 kbmax=kbmax,kbm=kbm,kmax=kmax, &
241 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
242 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
243 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
246 pratec1(i,j)=rain(i)*1000./(stepcu*dt)
247 raincv1(i,j)=rain(i)*1000./(stepcu)
251 raincv(i,j) = raincv1(i,j)
252 pratec(i,j) = pratec1(i,j)
257 IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
260 rthcuten(i,k,j)=(t1(i,k)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
261 rqvcuten(i,k,j)=(q1(i,k)-qv3d(i,k,j))*rdelt
266 IF(PRESENT(rucuten).AND.PRESENT(rvcuten)) THEN
269 rucuten(i,k,j)=(u1(i,k)-u3d(i,k,j))*rdelt
270 rvcuten(i,k,j)=(v1(i,k)-v3d(i,k,j))*rdelt
275 IF(PRESENT( rqicuten )) THEN
279 rqicuten(i,k,j)=(qi2(i,k)-qi3d(i,k,j))*rdelt
285 IF(PRESENT( rqccuten )) THEN
289 rqccuten(i,k,j)=(qc2(i,k)-qc3d(i,k,j))*rdelt
295 enddo ! outer most J_loop
298 end subroutine cu_ksas
300 !-------------------------------------------------------------------------------
301 ! NCEP SAS (Deep Convection Scheme)
302 !-------------------------------------------------------------------------------
303 subroutine nsas2d(delt,delx,del,prsl,prsi,prslk,zl, &
306 q1,t1,rain,kbot,ktop, &
308 lat,slimsk,dot,u1,v1,cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2, &
313 ids,ide, jds,jde, kds,kde, &
314 ims,ime, jms,jme, kms,kme, &
315 its,ite, jts,jte, kts,kte)
316 !-------------------------------------------------------------------------------
318 ! subprogram: phys_cps_sas computes convective heating and moistening
319 ! and momentum transport
321 ! abstract: computes convective heating and moistening using a one
322 ! cloud type arakawa-schubert convection scheme originally developed
323 ! by georg grell. the scheme has been revised at ncep since initial
324 ! implementation in 1993. it includes updraft and downdraft effects.
325 ! the closure is the cloud work function. both updraft and downdraft
326 ! are assumed to be saturated and the heating and moistening are
327 ! accomplished by the compensating environment. convective momentum transport
328 ! is taken into account. the name comes from "simplified arakawa-schubert
329 ! convection parameterization".
331 ! developed by hua-lu pan, wan-shu wu, songyou hong, and jongil han
332 ! implemented into wrf by kyosun sunny lim and songyou hong
333 ! module with cpp-based options is available in grims
335 ! program history log:
336 ! 92-03-01 hua-lu pan operational development
337 ! 96-03-01 song-you hong revised closure, and trigger
338 ! 99-03-01 hua-lu pan multiple clouds
339 ! 06-03-01 young-hwa byun closure based on moisture convergence (optional)
340 ! 09-10-01 jung-eun kim f90 format with standard physics modules
341 ! 10-07-01 jong-il han revised cloud model,trigger, as in gfs july 2010
342 ! 10-12-01 kyosun sunny lim wrf compatible version
343 ! 14-01-09 song-you hong dx dependent trigger, closure, and mass flux
344 ! 15-02-26 song-you hong negative qv generation suppressed
345 ! 15-06-01 jongil han gfs sas
346 ! 15-06-07 song-you hong scale-aware cps
347 ! 15-06-28 song-you hong diurnal cycle with cldwrk in pbl
348 ! 15-10-01 ji-young han bug fix (dz to dz1)
349 ! 16-01-01 ji-young han enhanced entrainment at lower RH
350 ! 16-06-01 ji-young han moisture-based trigger threshold
351 ! 16-09-01 ji-young han bug fix & eliminate inconsistency
352 ! 16-09-23 ji-young han code clean-up
353 ! 16-11-01 ji-young han revised pgcon & bug fix in vshear
354 ! 17-02-23 ji-young han revised xlamb
355 ! 18-03-01 ji-young han kim sas
356 ! 19-03-01 ji-yeon jang revised triggering ftn for higher resolution<10km
359 ! eun-jeong lee bug fix in c0 for overshooting layer
360 ! seung-bu park bug fix in xmb trigger
362 ! usage: call phys_cps_sas(delt,delx,del,prsl,prsi,prslk,prsik,zl, &
363 ! q2,q1,t1,u1,v1,rcs,slimsk,dot,cldwrk,rain, &
364 ! jcap,ncloud,lat,kbot,ktop,icps, &
365 ! ids,ide, jds,jde, kds,kde, &
366 ! ims,ime, jms,jme, kms,kme, &
367 ! its,ite, jts,jte, kts,kte)
369 ! delt - real model integration time step
370 ! delx - real model grid interval
371 ! del - real (kms:kme) sigma layer thickness
372 ! prsl - real (ims:ime,kms:kme) pressure values
373 ! prsi - real (ims:ime,kms:kme) pressure values at interface level
374 ! prslk - real (ims:ime,kms:kme) pressure values to the kappa
375 ! prsik - real (ims:ime,kms:kme) pressure values to the kappa at interface lev.
376 ! zl - real (ims:ime,kms:kme) height above sea level
377 ! zi - real (ims:ime,kms:kme) height above sea level at interface level
379 ! slimsk - real (ims:ime) land(1),sea(0), ice(2) flag
380 ! dot - real (ims:ime,kms:kme) vertical velocity
381 ! jcap - integer wave number
382 ! ncloud - integer no_cloud(0),no_ice(1),cloud+ice(2)
383 ! lat - integer current latitude index
385 ! output argument list:
386 ! q2 - real (ims:ime,kms:kme) detrained hydrometeors in kg/kg
387 ! - in case of the --> qc2(cloud), qi2(ice)
388 ! q1 - real (ims:ime,kms:kme) adjusted specific humidity in kg/kg
389 ! t1 - real (ims:ime,kms:kme) adjusted temperature in kelvin
390 ! u1 - real (ims:ime,kms:kme) adjusted zonal wind in m/s
391 ! v1 - real (ims:ime,kms:kme) adjusted meridional wind in m/s
392 ! cldwrk - real (ims:ime) cloud work function
393 ! rain - real (ims:ime) convective rain in meters
394 ! kbot - integer (ims:ime) cloud bottom level
395 ! ktop - integer (ims:ime) cloud top level
396 ! icps - integer (ims:ime) bit flag indicating deep convection
398 ! subprograms called:
399 ! fpvs - function to compute saturation vapor pressure
401 ! remarks: function fpvs is inlined by fpp.
402 ! nonstandard automatic arrays are used.
405 ! pan and wu (1995, ncep office note)
406 ! hong and pan (1998, mon wea rev)
407 ! park and hong (2007,jmsj)
408 ! byun and hong (2007, mon wea rev)
409 ! han and pan (2011, wea. forecasting)
410 ! lim et al. (2014, wea. forecasting)
411 ! han et al. (2016, mon wea rev)
412 ! kwon and hong (2017, mon wea rev)
413 ! han and hong (2018, wea. forecasting)
415 !-------------------------------------------------------------------------------
417 !-------------------------------------------------------------------------------
419 ! model tunable parameters
421 real,parameter :: betal = 0.05, betas = 0.05
422 real,parameter :: c0 = 0.002, c1 = 0.002
423 real,parameter :: xlamdd = 1.0e-4, xlamde = 1.0e-4
424 real,parameter :: clam = 0.1, cxlamu = 1.0e-3
425 real,parameter :: aafac = 0.1
426 real,parameter :: dthk = 25.
427 real,parameter :: cinpcrmx = 240.,cinpcrmn = 120.
428 real,parameter :: cinacrmx = -120.
429 real,parameter :: bet1 = 1.875, cd1 = 0.506
430 real,parameter :: f1 = 2.0, gam1 = 0.5, tfac = 1.0
431 real,parameter :: dx1km = 1000., dx5km = 5000., dx250m = 250.
432 real,parameter :: edtmaxl = 0.3, edtmaxs = 0.3
433 real,parameter :: evfacts = 0.3, evfactl = 0.3
434 real,parameter :: tf=233.16,tcr=273.16,tcrf=1.0/(tcr-tf)
438 real :: cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
439 real :: pi_,qmin_,t0c_,cice,xlv0,xls,psat
440 integer :: dx_factor_nsas
443 ids,ide, jds,jde, kds,kde, &
444 ims,ime, jms,jme, kms,kme, &
445 its,ite, jts,jte, kts,kte
448 real :: del(its:ite,kts:kte), &
449 prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme), &
450 prsi(its:ite,kts:kte+1), &
451 zl(its:ite,kts:kte), &
452 q1(its:ite,kts:kte),t1(its:ite,kts:kte), &
453 u1(its:ite,kts:kte),v1(its:ite,kts:kte), &
454 qci(its:ite,kts:kte),qrs(its:ite,kts:kte), &
456 real :: qi2(its:ite,kts:kte)
457 real :: qc2(its:ite,kts:kte)
459 real :: rain(its:ite)
460 real :: hpbl(ims:ime),hpbl_hold(ims:ime)
461 integer :: kbot(its:ite),ktop(its:ite),icps(its:ite)
462 real :: slimsk(ims:ime)
464 ! local variables and arrays
466 integer :: i,k,kmax,kbmax,kbm,jmn,indx, kts1,kte1,kmax1,kk
467 real :: p(its:ite,kts:kte),pdot(its:ite),acrtfct(its:ite)
468 real :: pden(its:ite)
469 real :: zi(its:ite,kts:kte+1)
470 real :: uo(its:ite,kts:kte),vo(its:ite,kts:kte)
471 real :: to(its:ite,kts:kte),qo(its:ite,kts:kte)
472 real :: hcko(its:ite,kts:kte)
473 real :: qcko(its:ite,kts:kte),eta(its:ite,kts:kte)
474 real :: etad(its:ite,kts:kte)
475 real :: qrcdo(its:ite,kts:kte)
476 real :: pwo(its:ite,kts:kte),pwdo(its:ite,kts:kte)
477 real :: c0t(its:ite,kts:kte)
478 real :: c1t(its:ite,kts:kte)
479 real :: bb1, bb2, wucb
483 ! for updraft velocity calculation
485 real :: po1(its:ite,kts:kte),wu2(its:ite,kts:kte), &
486 buo(its:ite,kts:kte),drag(its:ite,kts:kte)
487 real :: wbar(its:ite),wc(its:ite),clear(its:ite)
489 real :: sigma, sigma_con
491 real :: dtconv(its:ite)
492 real :: deltv(its:ite),acrt(its:ite)
493 real :: apbl(its:ite),dtpbl(its:ite)
494 real :: qeso(its:ite,kts:kte)
495 real :: tvo(its:ite,kts:kte),dbyo(its:ite,kts:kte)
496 real :: heo(its:ite,kts:kte),heso(its:ite,kts:kte)
497 real :: qrcd(its:ite,kts:kte)
498 real :: dellah(its:ite,kts:kte),dellaq(its:ite,kts:kte)
500 integer :: kb(its:ite),kbcon(its:ite)
501 integer :: kbcon1(its:ite)
502 real :: hmax(its:ite),delq(its:ite)
503 real :: hkbo(its:ite),qkbo(its:ite)
504 integer :: lmin(its:ite),jmin(its:ite)
505 integer :: ktcon(its:ite)
506 integer :: ktcon1(its:ite)
507 integer :: kbdtr(its:ite)
508 real :: hmin(its:ite),pwavo(its:ite)
509 real :: aa1(its:ite),vshear(its:ite)
510 real :: qevap(its:ite)
512 real :: edt_s(its:ite)
513 real :: edto(its:ite),pwevo(its:ite)
514 real :: qcond(its:ite)
515 real :: hcdo(its:ite,kts:kte)
516 real :: qcdo(its:ite,kts:kte)
517 real :: xhkb(its:ite),xqkb(its:ite)
518 real :: xpwav(its:ite),xpwev(its:ite),xhcd(its:ite,kts:kte)
519 real :: xaa0(its:ite),f(its:ite),xk(its:ite)
521 real :: edtx(its:ite),xqcd(its:ite,kts:kte)
522 real :: hsbar(its:ite),xmbmax(its:ite)
523 real :: xlamb(its:ite,kts:kte),xlamd(its:ite)
524 real :: cina(its:ite)
525 real :: delhbar(its:ite),delqbar(its:ite),deltbar(its:ite)
526 real :: qcirs(its:ite,kts:kte)
527 real :: dellal(its:ite,kts:kte)
528 real :: rntot(its:ite),delqev(its:ite),delq2(its:ite)
530 real :: fent1(its:ite,kts:kte),fent2(its:ite,kts:kte)
531 real :: frh(its:ite,kts:kte)
532 real :: xlamud(its:ite),sumx(its:ite)
533 real :: frh_sum(its:ite),cinpcri(its:ite)
535 real :: ucko(its:ite,kts:kte),vcko(its:ite,kts:kte)
536 real :: ucdo(its:ite,kts:kte),vcdo(its:ite,kts:kte)
537 real :: dellau(its:ite,kts:kte),dellav(its:ite,kts:kte)
538 real :: delubar(its:ite),delvbar(its:ite)
539 real :: qlko_ktcon(its:ite)
541 real :: alpha,beta, &
543 el2orc,eps,fact1,fact2, &
545 real :: dz,dp,es,pprime,qs, &
546 dqsdp,desdt,dqsdt,gamma, &
547 c0fac,alpha1,beta1,ccn_f, &
549 factor,onemf,dz1,qrch,etah,qlk,qc,rfact,shear, &
550 e1,dh, edtmax,dhh,dg,aup,adw, &
551 dv1,dv2,dv3,dv1q,dv2q,dv3q, &
552 dv1u,dv2u,dv3u,dv1v,dv2v,dv3v, &
553 dellat,xdby,xqrch, xpw,xpwd, &
554 qrsk(its:ite,kts:kte),evef,ptem,ptem1
556 logical :: totflg, cnvflg(its:ite),flg(its:ite)
557 real :: pgcon(its:ite,kts:kte)
559 !-----------------------------------------------------------------------
561 ! define miscellaneous values
568 el2orc = hvap_*hvap_/(rv_*cp_)
570 fact1 = (cvap_-cliq_)/rv_
571 fact2 = hvap_/rv_-fact1*t0c_
575 dtmin = max(dt2,600.)
576 dtmax = max(dt2,10800.)
578 sigma_con = tan(0.4*pi_)/(dx5km-dx1km) ! 7.7 e-4 m-1
579 sigma = (1.-1./pi_*(atan(sigma_con*(delx-dx5km))+pi_/2.)) ! 1(1km),0.1(10km)
581 if (delx.lt.dx5km) then
582 sigma = min(sigma - 0.01684 * delx/1000. + 0.0842, 1.0)
585 cinpcr = (cinpcrmn + 0.5*(cinpcrmx-cinpcrmn)) * (1.-sigma)
615 ! Define top layer for search of the downdraft originating layer
616 ! and the maximum thetae for updraft
622 ! convert surface pressure to mb from cb
627 c1t(i,k) = c1 * sigma
647 p(i,k) = prsl(i,k) * 10.
658 uo(i,k) = u1(i,k) * rcs
659 vo(i,k) = v1(i,k) * rcs
665 ! p is pressure of the layer (mb)
666 ! t is temperature at t-dt (k)..tn
667 ! q is mixing ratio at t-dt (kg/kg)..qn
668 ! to is temperature at t+dt (k)... this is after advection and turbulan
669 ! qo is mixing ratio at t+dt (kg/kg)..q1
673 qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
674 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
675 qeso(i,k) = max(qeso(i,k),qmin_)
676 qo(i,k) = max(qo(i,k), 1.e-10 )
677 ! tvo(i,k) = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
681 ! compute moist static energy
685 heo(i,k) = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qo(i,k)
686 heso(i,k) = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qeso(i,k)
690 ! Determine level with largest moist static energy
691 ! This is the level where updraft starts
700 if(heo(i,k).gt.hmax(i)) then
708 if (qo(i,kb(i)).lt.qmin_) cnvflg(i) = .false.
714 dz = .5 * (zl(i,k+1) - zl(i,k))
715 dp = .5 * (p(i,k+1) - p(i,k))
716 es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
717 pprime = p(i,k+1) + (eps-1.) * es
718 qs = eps * es / pprime
719 dqsdp = - qs / pprime
720 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
721 dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
722 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
723 dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
724 dq = dqsdt * dt + dqsdp * dp
725 to(i,k) = to(i,k+1) + dt
726 qo(i,k) = qo(i,k+1) + dq
727 po = .5 * (p(i,k) + p(i,k+1))
728 qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
729 qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
730 qeso(i,k) = max(qeso(i,k),qmin_)
731 qo(i,k) = max(qo(i,k), 1.e-10)
732 frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), 1.)
733 heo(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
734 cp_ * to(i,k) + hvap_ * qo(i,k)
735 heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
736 cp_ * to(i,k) + hvap_ * qeso(i,k)
737 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
738 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
743 ! look for convective cloud base as the level of free convection
748 hkbo(i) = heo(i,indx)
760 if(flg(i).and.k.gt.kb(i)) then
762 if(hkbo(i).gt.hsbar(i)) then
771 if(kbcon(i).eq.kmax) cnvflg(i) = .false.
776 totflg = totflg .and. (.not. cnvflg(i))
783 if (k.ge.kb(i).and.k.le.kbcon(i)) then
784 frh_sum(i) = frh_sum(i) + (1-frh(i,k))
792 tem1 = p(i,kb(i)) - p(i,kbcon(i))
793 cinpcri(i) = cinpcr * frh_sum(i)/(kbcon(i)-kb(i)+1)
794 if (tem1.gt.cinpcri(i)) cnvflg(i) = .false.
800 totflg = totflg .and. (.not. cnvflg(i))
806 zi(i,k) = 0.5*(zl(i,k-1)+zl(i,k))
812 xlamb(i,k) = clam / zi(i,k+1)
816 ! assume that updraft entrainment rate above cloud base is
817 ! same as that at cloud base
821 if(cnvflg(i).and.(k.gt.kbcon(i))) then
822 xlamb(i,k) = xlamb(i,kbcon(i))
827 ! assume the detrainment rate for the updrafts to be same as
828 ! the entrainment rate at cloud base
832 xlamud(i) = xlamb(i,kbcon(i))
836 ! functions rapidly decreasing with height, mimicking a cloud ensemble
837 ! (Bechtold et al., 2008)
841 if(cnvflg(i).and.(k.gt.kbcon(i))) then
842 tem = qeso(i,k)/qeso(i,kbcon(i))
849 ! final entrainment rate as the sum of turbulent part and organized entrainment
850 ! depending on the environmental relative humidity
851 ! (Bechtold et al., 2008)
855 if(cnvflg(i).and.(k.ge.kbcon(i))) then
856 tem = cxlamu * frh(i,k) * fent2(i,k)
857 xlamb(i,k) = xlamb(i,k)*fent1(i,k) + tem
862 ! determine updraft mass flux
874 if(cnvflg(i).and.k.lt.kbcon(i).and.k.ge.kb(i)) then
875 dz = zi(i,k+2) - zi(i,k+1)
876 ptem = 0.5*(xlamb(i,k)+xlamb(i,k+1))-xlamud(i)
877 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
883 if(cnvflg(i).and.k.gt.kbcon(i)) then
884 dz = zi(i,k+1) - zi(i,k)
885 ptem = 0.5*(xlamb(i,k)+xlamb(i,k-1))-xlamud(i)
886 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
892 dz = zi(i,3) - zi(i,2)
893 ptem = 0.5*(xlamb(i,1)+xlamb(i,2))-xlamud(i)
894 eta(i,1) = eta(i,2) / (1. + ptem * dz)
898 ! work up updraft cloud properties
903 hcko(i,indx) = hkbo(i)
904 qcko(i,indx) = qkbo(i)
905 ucko(i,indx) = uo(i,indx)
906 vcko(i,indx) = vo(i,indx)
911 ! cloud property below cloud base is modified by the entrainment proces
915 if(cnvflg(i).and.k.gt.kb(i)) then
916 dz = zi(i,k+1) - zi(i,k)
917 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
918 tem1 = 0.5 * xlamud(i) * dz
919 factor = 1. + tem - tem1
920 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
921 (heo(i,k)+heo(i,k-1)))/factor
922 dbyo(i,k) = hcko(i,k) - heso(i,k)
927 ! taking account into convection inhibition due to existence of
928 ! dry layers below cloud base
937 if(flg(i).and.k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
946 if(kbcon1(i).eq.kmax) cnvflg(i) = .false.
952 tem = p(i,kbcon(i)) - p(i,kbcon1(i))
961 totflg = totflg .and. (.not. cnvflg(i))
965 ! calculate convective inhibition
970 if (k.gt.kb(i).and.k.lt.kbcon1(i)) then
971 dz1 = (zi(i,k+1) - zi(i,k))
972 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
973 rfact = 1. + fv_ * cp_ * gamma * to(i,k) / hvap_
974 cina(i) = cina(i) + dz1 * (g_ / (cp_ * to(i,k))) &
975 * dbyo(i,k) / (1. + gamma) * rfact
976 cina(i) = cina(i) + dz1 * g_ * fv_ * max(0.,(qeso(i,k) - qo(i,k)))
985 if (cina(i).lt.cinacr) cnvflg(i) = .false.
991 totflg = totflg .and. (.not. cnvflg(i))
995 ! determine cloud top
1006 if(dbyo(i,k).lt.0..and.flg(i).and.k.gt. kbcon1(i)) then
1016 if(cnvflg(i).and.(p(i,kbcon(i)) - p(i,ktcon(i))).lt.150.) &
1022 totflg = totflg .and. (.not. cnvflg(i))
1026 ! search for downdraft originating level above theta-e minimum
1030 hmin(i) = heo(i,kbcon1(i))
1038 if(cnvflg(i).and.k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
1045 ! make sure that jmin is within the cloud
1049 jmin(i) = min(lmin(i),ktcon(i)-1)
1050 jmin(i) = max(jmin(i),kbcon1(i)+1)
1051 if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
1052 if(jmin(i).le.kbcon(i)) cnvflg(i) = .false.
1056 ! specify upper limit of mass flux at cloud base
1061 dp = 1000. * del(i,k)
1062 xmbmax(i) = dp / (g_ * dt2)
1066 ! compute cloud moisture property and precipitation
1070 if (cnvflg(i).and.k.gt.kb(i)) then
1071 alpha1 = min((-0.7*log(100.)+24.)*0.0001,c0)
1074 if (to(i,k).gt.t0c_) then
1077 c0fac = alpha1*exp(beta1*(to(i,k)-t0c_))
1080 c0fac = max(0.0,c0fac)
1084 if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1085 dz1 = (zi(i,k+1) - zi(i,k))
1086 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1088 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1089 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz1
1090 tem1 = 0.5 * xlamud(i) * dz1
1091 factor = 1. + tem - tem1
1092 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
1093 (qo(i,k)+qo(i,k-1)))/factor
1094 qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1096 ! check if there is excess moisture to release latent heat
1098 if(qcirs(i,k).gt.0. .and. k.ge.kbcon(i)) then
1099 etah = .5 * (eta(i,k) + eta(i,k-1))
1100 if(ncloud.gt.0..and.k.gt.jmin(i)) then
1101 dp = 1000. * del(i,k)
1102 ptem = c0t(i,k) + c1t(i,k)
1103 qlk = qcirs(i,k) / (eta(i,k) + etah * ptem * dz1)
1104 dellal(i,k) = etah * c1t(i,k) * dz1 * qlk * g_ / dp
1106 qlk = qcirs(i,k) / (eta(i,k) + etah * c0t(i,k) * dz1)
1108 pwo(i,k) = etah * c0t(i,k) * dz1 * qlk
1111 pwavo(i) = pwavo(i) + pwo(i,k)
1112 buo(i,k) = buo(i,k) - g_ * qlk
1114 ! compute buoyancy and drag for updraft velocity
1116 if (k.ge.kbcon(i)) then
1117 rfact = 1. + fv_ * cp_ * gamma &
1119 buo(i,k) = buo(i,k) + (g_ / (cp_ * to(i,k))) &
1120 * dbyo(i,k) / (1. + gamma) &
1122 buo(i,k) = buo(i,k) + g_ * fv_ * &
1123 max(0.,(qeso(i,k) - qo(i,k)))
1124 drag(i,k) = max(xlamb(i,k),xlamud(i))
1131 ! calculate cloud work function at t+dt
1142 if (k.ge.kbcon(i) .and. k.lt.ktcon(i)) then
1143 dz1 = zl(i,k+1) - zl(i,k)
1144 aa1(i) = aa1(i) + buo(i,k) * dz1
1151 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1156 totflg = totflg .and. (.not. cnvflg(i))
1160 ! estimate the convective overshooting as the level
1161 ! where the [aafac * cloud work function] becomes zero,
1162 ! which is the final cloud top
1166 aa2(i) = aafac * aa1(i)
1178 if(k.ge.ktcon(i).and.k.lt.kmax) then
1179 dz1 = zl(i,k+1) - zl(i,k)
1180 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1181 rfact = 1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1182 aa2(i) = aa2(i) +dz1 * (g_ / (cp_ * to(i,k))) &
1183 * dbyo(i,k) / (1. + gamma)* rfact
1184 if(aa2(i).lt.0.) then
1193 ! compute cloud moisture property, detraining cloud water
1194 ! and precipitation in overshooting layers
1199 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
1200 dz = (zi(i,k+1) - zi(i,k))
1201 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1202 qrch = qeso(i,k)+ gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1203 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1204 tem1 = 0.5 * xlamud(i) * dz
1205 factor = 1. + tem - tem1
1206 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
1207 (qo(i,k)+qo(i,k-1)))/factor
1208 qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1210 ! check if there is excess moisture to release latent heat
1212 if(qcirs(i,k).gt.0.) then
1213 etah = .5 * (eta(i,k) + eta(i,k-1))
1214 if(ncloud.gt.0.) then
1215 dp = 1000. * del(i,k)
1216 qlk = qcirs(i,k) / (eta(i,k) + etah * (c0t(i,k) + c1t(i,k)) * dz)
1217 dellal(i,k) = etah * c1t(i,k) * dz * qlk * g_ / dp
1219 qlk = qcirs(i,k) / (eta(i,k) + etah * c0t(i,k) * dz)
1221 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1224 pwavo(i) = pwavo(i) + pwo(i,k)
1231 ! compute updraft velocity square(wu2)
1233 bb1 = 2. * (1.+bet1*cd1)
1234 bb2 = 2. / (f1*(1.+gam1))
1242 po = .5 * (p(i,k) + p(i,k+1))
1243 tem = po / (rd_ * to(i,k))
1244 wucb = -10.*dot(i,k) / (tem * g_)
1245 if (wucb.gt.0.) then
1246 wu2(i,k) = wucb * wucb
1256 if (k.gt.kbcon1(i) .and. k.lt.ktcon(i)) then
1257 dz = zi(i,k+1) - zi(i,k)
1258 tem = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz
1259 tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz
1260 ptem = (1. - tem) * wu2(i,k-1)
1262 wu2(i,k) = (ptem + tem1) / ptem1
1263 wu2(i,k) = max(wu2(i,k), 0.)
1269 ! compute mean updraft velocity and mean grid-scale vertical velocity
1271 wc = 0. ; wbar = 0. ; sumx = 0.
1273 ptem = -0.5 * rd_ / g_
1276 po1(i,k) = .5 * (p(i,k) + p(i,k+1))
1283 if (k.gt.kbcon1(i) .and. k.lt.ktcon(i)) then
1284 dz = zi(i,k+1) - zi(i,k)
1285 tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1286 wc(i) = wc(i) + tem * dz
1287 tem = 10. * dot(i,k) * to(i,k) / po1(i,k)
1288 tem1 = 10. * dot(i,k-1) * to(i,k-1) / po1(i,k-1)
1289 wbar(i) = wbar(i) + ptem * (tem + tem1) * dz
1290 sumx(i) = sumx(i) + dz
1298 if (sumx(i) == 0.) then
1301 wc(i) = wc(i) / sumx(i)
1302 wbar(i) = wbar(i) / sumx(i)
1304 if (wc(i).lt.1.e-4) cnvflg(i) = .false.
1308 ! compute mean cloud core fraction
1309 ! assume mean cloud core fraction to be the ratio of
1310 ! mean grid-scale vertical velocity (wbar) and mean updraft velocity
1314 tem = wbar(i) / wc(i)
1317 clear(i) = max(min(clear(i), 1.0), 0.)
1318 if (wbar(i).gt.0. .and. wbar(i).gt.wc(i)) cnvflg(i) = .false.
1322 ! exchange ktcon with ktcon1
1327 ktcon(i) = ktcon1(i)
1333 if (ktcon(i).le.kb(i)) cnvflg(i) = .false.
1339 if (k.le.ktcon(i)) then
1340 pgcon(i,k) = 0.5+0.5*exp(3.*(zl(i,k)-zl(i,ktcon(i)))/zl(i,ktcon(i)))
1345 pgcon(i,k) = min(pgcon(i,k), 1.0)
1346 pgcon(i,k) = max(pgcon(i,k), 0.5)
1352 if (cnvflg(i) .and. k.gt.kb(i)) then
1353 dz = zi(i,k+1) - zi(i,k)
1354 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1355 tem1 = 0.5 * xlamud(i) * dz
1356 factor = 1. + tem - tem1
1357 ptem = 0.5 * tem + pgcon(i,k)
1358 ptem1 = 0.5 * tem - pgcon(i,k)
1359 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
1360 +ptem1*uo(i,k-1))/factor
1361 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
1362 +ptem1*vo(i,k-1))/factor
1367 ! this section is ready for cloud water
1369 if (ncloud.gt.0) then
1371 ! compute liquid and vapor separation at cloud top
1376 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1378 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1379 dq = qcko(i,k) - qrch
1381 ! check if there is excess moisture to release latent heat
1384 qlko_ktcon(i) = dq * sigma
1385 qcko(i,k) = qrch + dq * (1.-sigma)
1391 ! ..... downdraft calculations .....
1393 ! determine downdraft strength in terms of wind shear
1403 if(k.gt.kb(i).and.k.le.ktcon(i).and.cnvflg(i)) then
1404 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
1405 + (vo(i,k)-vo(i,k-1)) ** 2)
1406 vshear(i) = vshear(i) + shear
1413 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i)+1)-zi(i,kb(i)+1))
1414 e1 = 1.591-.639*vshear(i) &
1415 +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1420 edt_s(i) = min(edt_s(i),.9)
1421 edt_s(i) = max(edt_s(i),.0)
1422 edt(i) = min(edt(i),.9)
1423 edt(i) = max(edt(i),.0)
1424 edt(i) = edt(i) * ccn_f
1430 ! determine detrainment rate between 1 and kbdtr
1440 if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
1441 dz = zi(i,k+2) - zi(i,k+1)
1442 sumx(i) = sumx(i) + dz
1449 if(slimsk(i).eq.1.) beta = betal
1452 kbdtr(i) = max(kbdtr(i),1)
1453 dz =(sumx(i)+zi(i,2))/float(kbcon(i))
1454 tem = 1./float(kbcon(i))
1455 xlamd(i) = (1.-beta**tem)/dz
1459 ! determine downdraft mass flux
1474 if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
1475 dz = (zi(i,k+2) - zi(i,k+1))
1476 ptem = xlamdd-xlamde
1477 etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1478 elseif(k.lt.kbcon(i)) then
1479 dz = (zi(i,k+2) - zi(i,k+1))
1480 ptem = xlamd(i)+xlamdd-xlamde
1481 etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1487 ! downdraft moisture properties
1498 hcdo(i,jmn) = heo(i,jmn)
1499 qcdo(i,jmn) = qo(i,jmn)
1500 qrcdo(i,jmn) = qeso(i,jmn)
1501 ucdo(i,jmn) = uo(i,jmn)
1502 vcdo(i,jmn) = vo(i,jmn)
1508 if (cnvflg(i) .and. k.lt.jmin(i)) then
1509 dz = zi(i,k+2) - zi(i,k+1)
1510 if(k.ge.kbcon(i)) then
1512 tem1 = 0.5 * xlamdd * dz
1515 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1517 factor = 1. + tem - tem1
1518 ptem = 0.5 * tem - pgcon(i,k)
1519 ptem1 = 0.5 * tem + pgcon(i,k)
1520 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
1521 (heo(i,k)+heo(i,k+1)))/factor
1522 ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1) &
1523 +ptem1*uo(i,k))/factor
1524 vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1) &
1525 +ptem1*vo(i,k))/factor
1526 dbyo(i,k) = hcdo(i,k) - heso(i,k)
1533 if(cnvflg(i).and.k.lt.jmin(i)) then
1536 gamma = el2orc * dq / dt**2
1537 qrcdo(i,k) = dq+(1./hvap_)*(gamma/(1.+gamma))*dbyo(i,k)
1538 dz = zi(i,k+2) - zi(i,k+1)
1539 if(k.ge.kbcon(i)) then
1541 tem1 = 0.5 * xlamdd * dz
1544 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1546 factor = 1. + tem - tem1
1547 qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
1548 (qo(i,k)+qo(i,k+1)))/factor
1549 pwdo(i,k) = etad(i,k+1) * qcdo(i,k) -etad(i,k+1) * qrcdo(i,k)
1550 qcdo(i,k) = qrcdo(i,k)
1551 pwevo(i) = pwevo(i) + pwdo(i,k)
1556 ! final downdraft strength dependent on precip
1557 ! efficiency (edt), normalized condensate (pwav), and
1562 if(slimsk(i).eq.2.) edtmax = edtmaxs
1564 if(pwevo(i).lt.0.) then
1565 edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1566 edto(i) = min(edto(i),edtmax)
1573 ! downdraft cloudwork functions
1577 if(cnvflg(i).and.k.lt.jmin(i)) then
1578 gamma = el2orc * qeso(i,k) / to(i,k)**2
1583 dz = -1.*(zl(i,k+1)-zl(i,k))
1584 aa1(i) = aa1(i)+edto(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg)) &
1585 *(1.+fv_*cp_*dg*dt/hvap_)
1586 aa1(i) = aa1(i)+edto(i)*dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1592 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1597 totflg = totflg .and. (.not. cnvflg(i))
1601 ! what would the change be, that a cloud with unit mass
1602 ! will do to the environment?
1617 dp = 1000. * del(i,1)
1618 dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1) &
1619 - heo(i,1)) * g_ / dp
1620 dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1) &
1621 - qo(i,1)) * g_ / dp
1622 dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1) &
1623 - uo(i,1)) * g_ / dp
1624 dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1) &
1625 - vo(i,1)) * g_ / dp
1629 ! changed due to subsidence and entrainment
1633 if(cnvflg(i).and.k.lt.ktcon(i)) then
1635 if(k.le.kb(i)) aup = 0.
1637 if(k.gt.jmin(i)) adw = 0.
1639 dv2 = .5 * (heo(i,k) + heo(i,k-1))
1642 dv2q = .5 * (qo(i,k) + qo(i,k-1))
1645 dv2u = .5 * (uo(i,k) + uo(i,k-1))
1648 dv2v = .5 * (vo(i,k) + vo(i,k-1))
1650 dp = 1000. * del(i,k)
1651 dz = zi(i,k+1) - zi(i,k)
1652 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1))
1654 if(k.le.kbcon(i)) then
1656 ptem1 = xlamd(i)+xlamdd
1662 dellah(i,k) = dellah(i,k) + &
1663 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1 &
1664 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3 &
1665 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2*dz &
1666 + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
1667 + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz) *g_/dp
1669 dellaq(i,k) = dellaq(i,k) + &
1670 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1q &
1671 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3q &
1672 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz &
1673 + aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
1674 + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1))*dz) *g_/dp
1676 dellau(i,k) = dellau(i,k) + &
1677 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1u &
1678 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3u &
1679 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz &
1680 + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
1681 + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz &
1682 - pgcon(i,k)*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u))*g_/dp
1684 dellav(i,k) = dellav(i,k) + &
1685 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1v &
1686 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3v &
1687 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz &
1688 + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
1689 + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz &
1690 - pgcon(i,k)*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v))*g_/dp
1700 dp = 1000. * del(i,indx)
1702 dellah(i,indx) = eta(i,indx-1) * &
1703 (hcko(i,indx-1) - dv1) * g_ / dp
1705 dellaq(i,indx) = eta(i,indx-1) * &
1706 (qcko(i,indx-1) - dv1q) * g_ / dp
1708 dellau(i,indx) = eta(i,indx-1) * &
1709 (ucko(i,indx-1) - dv1u) * g_ / dp
1711 dellav(i,indx) = eta(i,indx-1) * &
1712 (vcko(i,indx-1) - dv1v) * g_ / dp
1716 dellal(i,indx) = eta(i,indx-1) * qlko_ktcon(i) * g_ / dp
1720 ! final changed variable per unit mass flux
1724 if(cnvflg(i).and.k.gt.ktcon(i)) then
1728 if(cnvflg(i).and.k.le.ktcon(i)) then
1729 qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
1730 dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1731 to(i,k) = dellat * mbdt + t1(i,k)
1732 qo(i,k) = max(qo(i,k),1.0e-10)
1737 !------------------------------------------------------------------------------
1739 ! the above changed environment is now used to calulate the
1740 ! effect the arbitrary cloud (with unit mass flux)
1741 ! which then is used to calculate the real mass flux,
1742 ! necessary to keep this change in balance with the large-scale
1745 ! environmental conditions again
1747 !------------------------------------------------------------------------------
1752 qeso(i,k)=0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1753 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1754 qeso(i,k) = max(qeso(i,k),qmin_)
1755 ! tvo(i,k) = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
1760 ! moist static energy
1765 dz = .5 * (zl(i,k+1) - zl(i,k))
1766 dp = .5 * (p(i,k+1) - p(i,k))
1767 es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1768 pprime = p(i,k+1) + (eps-1.) * es
1769 qs = eps * es / pprime
1770 dqsdp = - qs / pprime
1771 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1772 dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
1773 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1774 dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
1775 dq = dqsdt * dt + dqsdp * dp
1776 to(i,k) = to(i,k+1) + dt
1777 qo(i,k) = qo(i,k+1) + dq
1778 po = .5 * (p(i,k) + p(i,k+1))
1779 qeso(i,k) = 0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1780 qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
1781 qeso(i,k) = max(qeso(i,k),qmin_)
1782 qo(i,k) = max(qo(i,k), 1.0e-10)
1783 heo(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
1784 cp_ * to(i,k) + hvap_ * qo(i,k)
1785 heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
1786 cp_ * to(i,k) + hvap_ * qeso(i,k)
1794 heo(i,k) = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qo(i,k)
1795 heso(i,k) = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qeso(i,k)
1804 xhkb(i) = heo(i,indx)
1805 xqkb(i) = qo(i,indx)
1806 hcko(i,indx) = xhkb(i)
1807 qcko(i,indx) = xqkb(i)
1811 ! ..... static control .....
1813 ! moisture and cloud work functions
1817 if(cnvflg(i).and.k.gt.kb(i).and.k.le.ktcon(i)) then
1818 dz = zi(i,k+1) - zi(i,k)
1819 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1820 tem1 = 0.5 * xlamud(i) * dz
1821 factor = 1. + tem - tem1
1822 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
1823 (heo(i,k)+heo(i,k-1)))/factor
1830 if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1831 alpha1 = min((-0.7*log(100.)+24.)*0.0001,c0)
1834 if (to(i,k).gt.t0c_) then
1837 c0fac = alpha1*exp(beta1*(to(i,k)-t0c_))
1840 c0fac = max(0.0,c0fac)
1843 dz = zi(i,k+1) - zi(i,k)
1844 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1845 xdby = hcko(i,k) - heso(i,k)
1847 + gamma * xdby / (hvap_ * (1. + gamma))
1848 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1849 tem1 = 0.5 * xlamud(i) * dz
1850 factor = 1. + tem - tem1
1851 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*(qo(i,k)+qo(i,k-1)))/factor
1852 dq = eta(i,k) * qcko(i,k) - eta(i,k) * xqrch
1853 if(k.ge.kbcon(i).and.dq.gt.0.) then
1854 etah = .5 * (eta(i,k) + eta(i,k-1))
1855 if(ncloud.gt.0..and.k.gt.jmin(i)) then
1856 qlk = dq / (eta(i,k) + etah * (c0t(i,k) + c1t(i,k)) * dz)
1858 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1860 xpw = etah * c0t(i,k) * dz * qlk
1861 if(k.lt.ktcon1(i)) then
1862 xaa0(i) = xaa0(i) - (zl(i,k+1) - zl(i,k)) * g_ * qlk
1864 qcko(i,k) = qlk + xqrch
1865 xpwav(i) = xpwav(i) + xpw
1869 if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
1870 dz1 = zl(i,k+1) - zl(i,k)
1871 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1872 rfact = 1. + fv_ * cp_ * gamma &
1874 xdby = hcko(i,k) - heso(i,k)
1876 + dz1 * (g_ / (cp_ * to(i,k))) &
1877 * xdby / (1. + gamma) &
1881 max(0.,(qeso(i,k) - qo(i,k)))
1886 ! ..... downdraft calculations .....
1888 ! downdraft moisture properties
1897 xhcd(i,jmn) = heo(i,jmn)
1898 xqcd(i,jmn) = qo(i,jmn)
1899 qrcd(i,jmn) = qeso(i,jmn)
1905 if(cnvflg(i).and.k.lt.jmin(i)) then
1906 dz = zi(i,k+2) - zi(i,k+1)
1907 if(k.ge.kbcon(i)) then
1909 tem1 = 0.5 * xlamdd * dz
1912 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1914 factor = 1. + tem - tem1
1915 xhcd(i,k) = ((1.-tem1)*xhcd(i,k+1)+tem*0.5* &
1916 (heo(i,k)+heo(i,k+1)))/factor
1923 if(cnvflg(i).and.k.lt.jmin(i)) then
1926 gamma = el2orc * dq / dt**2
1927 dh = xhcd(i,k) - heso(i,k)
1928 qrcd(i,k) = dq+(1./hvap_)*(gamma/(1.+gamma))*dh
1929 dz = zi(i,k+2) - zi(i,k+1)
1930 if(k.ge.kbcon(i)) then
1932 tem1 = 0.5 * xlamdd * dz
1935 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1937 factor = 1. + tem - tem1
1938 xqcd(i,k) = ((1.-tem1)*xqcd(i,k+1)+tem*0.5* &
1939 (qo(i,k)+qo(i,k+1)))/factor
1940 xpwd = etad(i,k+1) * (xqcd(i,k) - qrcd(i,k))
1941 xqcd(i,k)= qrcd(i,k)
1942 xpwev(i) = xpwev(i) + xpwd
1949 if(slimsk(i).eq.2.) edtmax = edtmaxs
1951 if(xpwev(i).ge.0.) then
1954 edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1955 edtx(i) = min(edtx(i),edtmax)
1960 ! downdraft cloudwork functions
1964 if(cnvflg(i).and.k.lt.jmin(i)) then
1965 gamma = el2orc * qeso(i,k) / to(i,k)**2
1970 dz =-1.*(zl(i,k+1)-zl(i,k))
1971 xaa0(i) = xaa0(i)+edtx(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg)) &
1972 *(1.+fv_*cp_*dg*dt/hvap_)
1973 xaa0(i) = xaa0(i)+edtx(i)* &
1974 dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1979 ! calculate critical cloud work function
1982 tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
1983 dtconv(i) = tfac * tem / wc(i)
1984 dtconv(i) = max(dtconv(i),dtmin)
1985 dtconv(i) = min(dtconv(i),dtmax)
1986 ! dtconv(i) = max(1800., dt2)
1989 ! large scale forcing
1991 ! compute mean velocity between kb kbcon
1994 if (cnvflg(i).and.slimsk(i).eq.2.) then
2002 if (cnvflg(i).and.slimsk(i).eq.2.) then
2003 if (k.gt.kb(i) .and. k.le.kbcon(i)) then
2004 dz = zi(i,k+1) - zi(i,k)
2005 tem = sqrt(uo(i,k)**2 + vo(i,k)**2)
2006 wc(i) = wc(i) + tem * dz
2007 sumx(i) = sumx(i) + dz
2014 if (cnvflg(i).and.slimsk(i).eq.2.) then
2015 wc(i) = wc(i) / sumx(i)
2016 dtpbl(i) = zl(i,kbcon(i))/wc(i)
2018 dtpbl(i) = dtconv(i)
2024 apbl(i) = g_ * (hpbl(i)-hpbl_hold(i))/dt2 * dtpbl(i)
2025 apbl(i) = min(max(apbl(i),0.0),aa1(i))
2031 f(i) = (aa1(i) - apbl(i)) / dtconv(i)
2032 if(f(i).le.0.) cnvflg(i) = .false.
2035 xk(i) = (xaa0(i) - aa1(i)) / mbdt
2036 if(xk(i).ge.0.) cnvflg(i) = .false.
2039 ! kernel, cloud base mass flux
2042 xmb(i) = -f(i) / xk(i)
2043 xmb(i) = xmb(i) * clear(i) * (1.-sigma)
2044 xmb(i) = min(xmb(i),xmbmax(i))
2046 pden(i) = p(i,kbcon(i))/to(i,kbcon(i))/rd_
2047 pdot(i) = 10.* dot(i,kbcon(i))
2048 ! if (pden(i)*pdot(i).gt.xmb(i)) cnvflg(i) = .false.
2049 if (xmb(i).lt.pdot(i)/g_) cnvflg(i) = .false.
2053 totflg = totflg .and. (.not. cnvflg(i))
2057 ! restore t0 and qo to t1 and q1 in case convection stops
2066 qeso(i,k) = 0.01*fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2067 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
2068 qeso(i,k) = max(qeso(i,k),qmin_)
2073 ! feedback: simply the changes from the cloud with unit mass flux
2074 ! multiplied by the mass flux necessary to keep the
2075 ! equilibrium with the larger-scale.
2088 if (cnvflg(i).and.k.le.ktcon(i).and.dellaq(i,k).le.0.) then
2089 if (q1(i,k).gt.0.) then
2090 tem = dellaq(i,k) * xmb(i) * dt2
2091 dellaq(i,k) = max(tem,-q1(i,k))/(xmb(i)*dt2)
2101 if(cnvflg(i).and.k.le.ktcon(i)) then
2102 dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
2103 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
2104 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
2106 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
2107 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
2108 dp = 1000. * del(i,k)
2109 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
2110 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
2111 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
2112 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
2113 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
2120 if (cnvflg(i) .and. k.le.ktcon(i)) then
2121 qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2122 qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
2123 qeso(i,k) = max(qeso(i,k), qmin_)
2139 if(cnvflg(i).and.k.lt.ktcon(i)) then
2141 if(k.le.kb(i)) aup = 0.
2143 if(k.ge.jmin(i)) adw = 0.
2144 rntot(i) = rntot(i) &
2145 + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)) &
2146 * xmb(i) * .001 * dt2
2151 ! conversion rainfall (m) and compute the evaporation of falling raindrops
2158 if(cnvflg(i).and.k.lt.ktcon(i)) then
2160 if(k.le.kb(i)) aup = 0.
2162 if(k.ge.jmin(i)) adw = 0.
2164 + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)) &
2165 * xmb(i) * .001 * dt2
2166 qrsk(i,k) = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
2169 if(cnvflg(i).and.flg(i).and.k.lt.ktcon(i)) then
2171 evef = edt_s(i) * evfacts * ccn_f
2172 if(slimsk(i).eq.1.) evef = edt_s(i) * evfactl * ccn_f
2173 qcond(i) = evef * (q1(i,k) - qeso(i,k)) / (1. + el2orc * &
2174 qeso(i,k) / t1(i,k)**2)
2175 dp = 1000. * del(i,k)
2176 if(rain(i).gt.0..and.qcond(i).lt.0.) then
2177 qevap(i) = -qcond(i) * (1. - exp(-.32 * sqrt(dt2 * rain(i))))
2178 qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
2179 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
2180 if (delq2(i).gt.rntot(i)) then
2181 qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
2185 if(rain(i).gt.0..and.qevap(i).gt.0.) then
2186 q1(i,k) = q1(i,k) + qevap(i)
2187 t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
2188 rain(i) = rain(i) - .001 * qevap(i) * dp / g_
2189 delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
2190 deltv(i) = - (hvap_/cp_)*qevap(i)/dt2
2191 delq(i) = + qevap(i)/dt2
2193 dellaq(i,k) = dellaq(i,k) + delq(i)/xmb(i)
2194 delqbar(i) = delqbar(i) + delq(i)*dp/g_
2195 deltbar(i) = deltbar(i) + deltv(i)*dp/g_
2198 if (cnvflg(i).and.k.lt.ktcon(i)) then
2199 qrs(i,k) = max(qrsk(i,k) - max(qevap(i),0.),0.)
2200 qci(i,k) = max(qcirs(i,k) - aup*pwo(i,k),0.)
2205 ! consider the negative rain in the event of rain evaporation and downdrafts
2209 if(rain(i).lt.0..and..not.flg(i)) rain(i) = 0.
2210 if(rain(i).le.0.) then
2222 if(cnvflg(i).and.rain(i).le.0.) then
2231 ! detrainment of cloud water and ice
2233 if (ncloud.gt.0) then
2236 if (cnvflg(i) .and. rain(i).gt.0.) then
2237 if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
2238 tem = dellal(i,k) * xmb(i) * dt2
2239 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2240 if (ncloud.ge.2) then
2241 qi2(i,k) = qi2(i,k) + tem * tem1 ! ice
2242 qc2(i,k) = qc2(i,k) + tem *(1.0-tem1) ! water
2244 qc2(i,k) = qc2(i,k) + tem
2252 end subroutine nsas2d
2253 !-------------------------------------------------------------------------------
2254 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
2255 !-------------------------------------------------------------------------------
2257 !-------------------------------------------------------------------------------
2258 REAL :: t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
2268 xbi=xai+hsub/(rv*ttp)
2270 if(t.lt.ttp.and.ice.eq.1) then
2271 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2273 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2278 if(t.lt.ttp.and.ice.eq.1) then
2279 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2281 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2287 if(t.lt.ttp.and.ice.eq.1) then
2288 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2290 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2295 !-------------------------------------------------------------------------------
2297 !-------------------------------------------------------------------------------
2298 subroutine nsasinit(rthcuten,rqvcuten,rqccuten,rqicuten, &
2300 restart,p_qc,p_qi,p_first_scalar, &
2302 ids, ide, jds, jde, kds, kde, &
2303 ims, ime, jms, jme, kms, kme, &
2304 its, ite, jts, jte, kts, kte )
2305 !-------------------------------------------------------------------------------
2307 !-------------------------------------------------------------------------------
2308 logical , intent(in) :: allowed_to_read,restart
2309 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
2310 ims, ime, jms, jme, kms, kme, &
2311 its, ite, jts, jte, kts, kte
2312 integer , intent(in) :: p_first_scalar, p_qi, p_qc
2313 real, dimension( ims:ime , kms:kme , jms:jme ) , intent(out) :: &
2320 integer :: i, j, k, itf, jtf, ktf
2326 if(.not.restart)then
2338 if (p_qc .ge. p_first_scalar) then
2348 if (p_qi .ge. p_first_scalar) then
2359 end subroutine nsasinit
2360 !-------------------------------------------------------------------------------
2362 END MODULE module_cu_ksas