7 !-------------------------------------------------------------------------------
8 subroutine cu_nsas(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, &
14 mp_physics,dx_factor_nsas, &
15 p_qc,p_qi,p_first_scalar, &
17 cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2, &
18 cice,xls,psat,f_qi,f_qc, &
19 ids,ide, jds,jde, kds,kde, &
20 ims,ime, jms,jme, kms,kme, &
21 its,ite, jts,jte, kts,kte)
22 !-------------------------------------------------------------------------------
24 !-------------------------------------------------------------------------------
26 !-- dx grid interval (m)
27 !-- p3di 3d pressure (pa) at interface level
28 !-- p3d 3d pressure (pa)
29 !-- pi3d 3d exner function (dimensionless)
30 !-- z height above sea level (m)
31 !-- qc3d cloud water mixing ratio (kg/kg)
32 !-- qi3d cloud ice mixing ratio (kg/kg)
33 !-- qv3d 3d water vapor mixing ratio (kg/kg)
34 !-- t3d temperature (k)
35 !-- raincv cumulus scheme precipitation (mm)
36 !-- w vertical velocity (m/s)
37 !-- dz8w dz between full levels (m)
38 !-- u3d 3d u-velocity interpolated to theta points (m/s)
39 !-- v3d 3d v-velocity interpolated to theta points (m/s)
40 !-- ids start index for i in domain
41 !-- ide end index for i in domain
42 !-- jds start index for j in domain
43 !-- jde end index for j in domain
44 !-- kds start index for k in domain
45 !-- kde end index for k in domain
46 !-- ims start index for i in memory
47 !-- ime end index for i in memory
48 !-- jms start index for j in memory
49 !-- jme end index for j in memory
50 !-- kms start index for k in memory
51 !-- kme end index for k in memory
52 !-- its start index for i in tile
53 !-- ite end index for i in tile
54 !-- jts start index for j in tile
55 !-- jte end index for j in tile
56 !-- kts start index for k in tile
57 !-- kte end index for k in tile
58 !-------------------------------------------------------------------------------
59 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
60 ims,ime, jms,jme, kms,kme, &
61 its,ite, jts,jte, kts,kte, &
63 p_qc,p_qi,p_first_scalar
64 real, intent(in ) :: cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2, &
66 real, intent(in ) :: dt,dx
67 real, optional, intent(in ) :: pgcon
68 real, dimension( ims:ime, kms:kme, jms:jme ),optional ,&
69 intent(inout) :: rthcuten,&
75 logical, optional :: F_QC,F_QI
76 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
84 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
86 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
89 real, dimension( ims:ime, jms:jme ) ,&
90 intent(inout) :: raincv,&
92 real, dimension( ims:ime, jms:jme ) ,&
96 real, dimension( ims:ime, jms:jme ) ,&
99 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
102 logical, dimension( ims:ime, jms:jme ) ,&
103 intent(inout) :: cu_act_flag
105 real, dimension( ims:ime, jms:jme ) ,&
106 intent(in ) :: hpbl,&
109 integer, intent(in ) :: mp_physics
110 integer, intent(in ) :: dx_factor_nsas
115 real, dimension( its:ite, jts:jte ) :: raincv1,&
119 real, dimension( its:ite, kts:kte ) :: del,&
128 real, dimension( its:ite, kts:kte+1 ) :: prsii,&
130 real, dimension( its:ite, kts:kte ) :: zll
131 real, dimension( its:ite) :: rain
133 integer, dimension (its:ite) :: kbot,&
139 ! microphysics scheme --> ncloud
141 if (mp_physics .eq. 0) then
143 elseif ( mp_physics .eq. 1 .or. mp_physics .eq. 3 ) then
149 if(present(pgcon)) then
152 ! pgcon_use = 0.7 ! Gregory et al. (1997, QJRMS)
153 pgcon_use = 0.55 ! Zhang & Wu (2003,JAS)
154 ! 0.55 is a physically-based value used by GFS
155 ! HWRF uses 0.2, for model tuning purposes
160 cu_act_flag(i,j)=.TRUE.
172 dot(i,k) = -5.0e-4*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
173 prsll(i,k)=p3d(i,k,j)*0.001
174 prsii(i,k)=p3di(i,k,j)*0.001
179 prsii(i,kte+1)=p3di(i,kte+1,j)*0.001
188 zii(i,k+1)=zii(i,k)+dz8w(i,k,j)
194 zll(i,k)=0.5*(zii(i,k)+zii(i,k+1))
200 del(i,k)=prsll(i,k)*g/r_d*dz8w(i,k,j)/t3d(i,k,j)
204 ! q1(i,k)=qv3d(i,k,j)/(1.+qv3d(i,k,j))
206 qi2(i,k) = qi3d(i,k,j)
207 qc2(i,k) = qc3d(i,k,j)
213 call nsas2d(delt=delt,delx=dx,del=del(its,kts), &
214 prsl=prsll(its,kts),prsi=prsii(its,kts),prslk=pi3d(ims,kms,j), &
216 ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
217 q1=q1(its,kts),t1=t1(its,kts),rain=rain(its), &
218 kbot=kbot(its),ktop=ktop(its), &
220 lat=j,slimsk=xland(ims,j),dot=dot(its,kts), &
221 u1=u1(its,kts), v1=v1(its,kts), &
222 cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv, &
223 rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2, &
224 cice=cice,xls=xls,psat=psat, &
226 dx_factor_nsas=dx_factor_nsas, &
227 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
228 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
229 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
232 pratec1(i,j)=rain(i)*1000./(stepcu*dt)
233 raincv1(i,j)=rain(i)*1000./(stepcu)
238 call nscv2d(delt=delt,del=del(its,kts),prsl=prsll(its,kts), &
239 prsi=prsii(its,kts),prslk=pi3d(ims,kms,j),zl=zll(its,kts), &
240 ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
241 q1=q1(its,kts),t1=t1(its,kts),rain=rain(its), &
242 kbot=kbot(its),ktop=ktop(its), &
244 slimsk=xland(ims,j),dot=dot(its,kts), &
245 u1=u1(its,kts), v1=v1(its,kts), &
246 cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv, &
247 rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2, &
248 cice=cice,xls=xls,psat=psat, &
249 hpbl=hpbl(ims,j),hfx=hfx(ims,j),qfx=qfx(ims,j), &
251 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
252 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
253 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
256 pratec2(i,j)=rain(i)*1000./(stepcu*dt)
257 raincv2(i,j)=rain(i)*1000./(stepcu)
261 raincv(i,j) = raincv1(i,j) + raincv2(i,j)
262 pratec(i,j) = pratec1(i,j) + pratec2(i,j)
267 IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
270 rthcuten(i,k,j)=(t1(i,k)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
271 rqvcuten(i,k,j)=(q1(i,k)-qv3d(i,k,j))*rdelt
276 IF(PRESENT(rucuten).AND.PRESENT(rvcuten)) THEN
279 rucuten(i,k,j)=(u1(i,k)-u3d(i,k,j))*rdelt
280 rvcuten(i,k,j)=(v1(i,k)-v3d(i,k,j))*rdelt
285 IF(PRESENT( rqicuten )) THEN
289 rqicuten(i,k,j)=(qi2(i,k)-qi3d(i,k,j))*rdelt
295 IF(PRESENT( rqccuten )) THEN
299 rqccuten(i,k,j)=(qc2(i,k)-qc3d(i,k,j))*rdelt
305 enddo ! outer most J_loop
308 end subroutine cu_nsas
310 !-------------------------------------------------------------------------------
311 ! NCEP SAS (Deep Convection Scheme)
312 !-------------------------------------------------------------------------------
313 subroutine nsas2d(delt,delx,del,prsl,prsi,prslk,zl, &
316 q1,t1,rain,kbot,ktop, &
318 lat,slimsk,dot,u1,v1,cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2, &
322 ids,ide, jds,jde, kds,kde, &
323 ims,ime, jms,jme, kms,kme, &
324 its,ite, jts,jte, kts,kte)
325 !-------------------------------------------------------------------------------
327 ! subprogram: phys_cps_sas computes convective heating and moistening
328 ! and momentum transport
330 ! abstract: computes convective heating and moistening using a one
331 ! cloud type arakawa-schubert convection scheme originally developed
332 ! by georg grell. the scheme has been revised at ncep since initial
333 ! implementation in 1993. it includes updraft and downdraft effects.
334 ! the closure is the cloud work function. both updraft and downdraft
335 ! are assumed to be saturated and the heating and moistening are
336 ! accomplished by the compensating environment. convective momentum transport
337 ! is taken into account. the name comes from "simplified arakawa-schubert
338 ! convection parameterization".
340 ! developed by hua-lu pan, wan-shu wu, songyou hong, and jongil han
341 ! implemented into wrf by kyosun sunny lim and songyou hong
342 ! module with cpp-based options is available in grims
344 ! program history log:
345 ! 92-03-01 hua-lu pan operational development
346 ! 96-03-01 song-you hong revised closure, and trigger
347 ! 99-03-01 hua-lu pan multiple clouds
348 ! 06-03-01 young-hwa byun closure based on moisture convergence (optional)
349 ! 09-10-01 jung-eun kim f90 format with standard physics modules
350 ! 10-07-01 jong-il han revised cloud model,trigger, as in gfs july 2010
351 ! 10-12-01 kyosun sunny lim wrf compatible version
352 ! 14-01-09 song-you hong dx dependent trigger, closure, and mass flux
355 ! usage: call phys_cps_sas(delt,delx,del,prsl,prsi,prslk,prsik,zl, &
356 ! q2,q1,t1,u1,v1,rcs,slimsk,dot,cldwrk,rain, &
357 ! jcap,ncloud,lat,kbot,ktop,icps, &
358 ! ids,ide, jds,jde, kds,kde, &
359 ! ims,ime, jms,jme, kms,kme, &
360 ! its,ite, jts,jte, kts,kte)
362 ! delt - real model integration time step
363 ! delx - real model grid interval
364 ! del - real (kms:kme) sigma layer thickness
365 ! prsl - real (ims:ime,kms:kme) pressure values
366 ! prsi - real (ims:ime,kms:kme) pressure values at interface level
367 ! prslk - real (ims:ime,kms:kme) pressure values to the kappa
368 ! prsik - real (ims:ime,kms:kme) pressure values to the kappa at interface lev.
369 ! zl - real (ims:ime,kms:kme) height above sea level
370 ! zi - real (ims:ime,kms:kme) height above sea level at interface level
372 ! slimsk - real (ims:ime) land(1),sea(0), ice(2) flag
373 ! dot - real (ims:ime,kms:kme) vertical velocity
374 ! jcap - integer wave number
375 ! ncloud - integer no_cloud(0),no_ice(1),cloud+ice(2)
376 ! lat - integer current latitude index
378 ! output argument list:
379 ! q2 - real (ims:ime,kms:kme) detrained hydrometeors in kg/kg
380 ! - in case of the --> qc2(cloud), qi2(ice)
381 ! q1 - real (ims:ime,kms:kme) adjusted specific humidity in kg/kg
382 ! t1 - real (ims:ime,kms:kme) adjusted temperature in kelvin
383 ! u1 - real (ims:ime,kms:kme) adjusted zonal wind in m/s
384 ! v1 - real (ims:ime,kms:kme) adjusted meridional wind in m/s
385 ! cldwrk - real (ims:ime) cloud work function
386 ! rain - real (ims:ime) convective rain in meters
387 ! kbot - integer (ims:ime) cloud bottom level
388 ! ktop - integer (ims:ime) cloud top level
389 ! icps - integer (ims:ime) bit flag indicating deep convection
391 ! subprograms called:
392 ! fpvs - function to compute saturation vapor pressure
394 ! remarks: function fpvs is inlined by fpp.
395 ! nonstandard automatic arrays are used.
398 ! pan and wu (1995, ncep office note)
399 ! hong and pan (1998, mon wea rev)
400 ! park and hong (2007,jmsj)
401 ! byun and hong (2007, mon wea rev)
402 ! han and pan (2011, wea. forecasting)
404 !-------------------------------------------------------------------------------
405 !-------------------------------------------------------------------------------
407 !-------------------------------------------------------------------------------
409 ! model tunable parameters
411 real,parameter :: alphal = 0.5, alphas = 0.5
412 real,parameter :: betal = 0.05, betas = 0.05
413 real,parameter :: pdpdwn = 0.0, pdetrn = 200.0
414 real,parameter :: c0 = 0.002, c1 = 0.002
415 real,parameter :: xlamdd = 1.0e-4, xlamde = 1.0e-4
416 real,parameter :: clam = 0.1, cxlamu = 1.0e-4
417 real,parameter :: aafac = 0.1
418 real,parameter :: dthk=25.
419 real,parameter :: cincrmax = 180.,cincrmin = 120.
420 real,parameter :: mbdt = 10., edtmaxl = 0.3, edtmaxs = 0.3
421 real,parameter :: evfacts = 0.3, evfactl = 0.3
423 real,parameter :: tf=233.16,tcr=263.16,tcrf=1.0/(tcr-tf)
424 real,parameter :: xk1=2.e-5,xlhor=3.e4,xhver=5000.,theimax=1.
425 real,parameter :: xc1=1.e-7,xc2=1.e4,xc3=3.e3,ecesscr=3.0,edtk1=3.e4
429 real :: cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
430 real :: pi_,qmin_,t0c_,cice,xlv0,xls,psat
432 integer :: dx_factor_nsas
435 ids,ide, jds,jde, kds,kde, &
436 ims,ime, jms,jme, kms,kme, &
437 its,ite, jts,jte, kts,kte
440 real :: del(its:ite,kts:kte), &
441 prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme), &
442 prsi(its:ite,kts:kte+1), &
443 zl(its:ite,kts:kte), &
444 q1(its:ite,kts:kte),t1(its:ite,kts:kte), &
445 u1(its:ite,kts:kte),v1(its:ite,kts:kte), &
447 real :: qi2(its:ite,kts:kte)
448 real :: qc2(its:ite,kts:kte)
450 real :: rain(its:ite)
451 integer :: kbot(its:ite),ktop(its:ite),icps(its:ite)
452 real :: slimsk(ims:ime)
455 ! local variables and arrays
457 integer :: i,k,kmax,kbmax,kbm,jmn,indx,indp,kts1,kte1,kmax1,kk
458 real :: p(its:ite,kts:kte),pdot(its:ite),acrtfct(its:ite)
459 real :: zi(its:ite,kts:kte+1)
460 real :: uo(its:ite,kts:kte),vo(its:ite,kts:kte)
461 real :: to(its:ite,kts:kte),qo(its:ite,kts:kte)
462 real :: hcko(its:ite,kts:kte)
463 real :: qcko(its:ite,kts:kte),eta(its:ite,kts:kte)
464 real :: etad(its:ite,kts:kte)
465 real :: qrcdo(its:ite,kts:kte)
466 real :: pwo(its:ite,kts:kte),pwdo(its:ite,kts:kte)
467 real :: dtconv(its:ite)
468 real :: deltv(its:ite),acrt(its:ite)
469 real :: qeso(its:ite,kts:kte)
470 real :: tvo(its:ite,kts:kte),dbyo(its:ite,kts:kte)
471 real :: heo(its:ite,kts:kte),heso(its:ite,kts:kte)
472 real :: qrcd(its:ite,kts:kte)
473 real :: dellah(its:ite,kts:kte),dellaq(its:ite,kts:kte)
475 integer :: kb(its:ite),kbcon(its:ite)
476 integer :: kbcon1(its:ite)
477 real :: hmax(its:ite),delq(its:ite)
478 real :: hkbo(its:ite),qkbo(its:ite),pbcdif(its:ite)
479 integer :: kbds(its:ite),lmin(its:ite),jmin(its:ite)
480 integer :: ktcon(its:ite)
481 integer :: ktcon1(its:ite)
482 integer :: kbdtr(its:ite)
483 integer :: klcl(its:ite),ktdown(its:ite)
484 real :: vmax(its:ite)
485 real :: hmin(its:ite),pwavo(its:ite)
486 real :: aa1(its:ite),vshear(its:ite)
487 real :: qevap(its:ite)
489 real :: edto(its:ite),pwevo(its:ite)
490 real :: qcond(its:ite)
491 real :: hcdo(its:ite,kts:kte)
492 real :: ddp(its:ite),pp2(its:ite)
493 real :: qcdo(its:ite,kts:kte)
494 real :: adet(its:ite),aatmp(its:ite)
495 real :: xhkb(its:ite),xqkb(its:ite)
496 real :: xpwav(its:ite),xpwev(its:ite),xhcd(its:ite,kts:kte)
497 real :: xaa0(its:ite),f(its:ite),xk(its:ite)
499 real :: edtx(its:ite),xqcd(its:ite,kts:kte)
500 real :: hsbar(its:ite),xmbmax(its:ite)
501 real :: xlamb(its:ite,kts:kte),xlamd(its:ite)
502 real :: excess(its:ite)
503 real :: plcl(its:ite)
504 real :: delhbar(its:ite),delqbar(its:ite),deltbar(its:ite)
505 real,save :: pcrit(15), acritt(15)
507 real :: qcirs(its:ite,kts:kte),qrski(its:ite)
508 real :: dellal(its:ite,kts:kte)
509 real :: rntot(its:ite),delqev(its:ite),delq2(its:ite)
511 real :: fent1(its:ite,kts:kte),fent2(its:ite,kts:kte)
512 real :: frh(its:ite,kts:kte)
513 real :: xlamud(its:ite),sumx(its:ite)
515 real :: ucko(its:ite,kts:kte),vcko(its:ite,kts:kte)
516 real :: ucdo(its:ite,kts:kte),vcdo(its:ite,kts:kte)
517 real :: dellau(its:ite,kts:kte),dellav(its:ite,kts:kte)
518 real :: delubar(its:ite),delvbar(its:ite)
519 real :: qlko_ktcon(its:ite)
521 real :: alpha,beta, &
522 dt2,dtmin,dtmax,dtmaxl,dtmaxs, &
523 el2orc,eps,fact1,fact2, &
525 real :: dz,dp,es,pprime,qs, &
526 dqsdp,desdt,dqsdt,gamma, &
527 dt,dq,po,thei,delx,delza,dzfac, &
528 thec,theb,thekb,thekh,theavg,thedif, &
529 omgkb,omgkbp1,omgdif,omgfac,heom,rh,thermal,chi, &
530 factor,onemf,dz1,qrch,etah,qlk,qc,rfact,shear, &
531 e1,dh,deta,detad,theom,edtmax,dhh,dg,aup,adw, &
532 dv1,dv2,dv3,dv1q,dv2q,dv3q,dvq1, &
533 dv1u,dv2u,dv3u,dv1v,dv2v,dv3v, &
534 dellat,xdby,xqrch,xqc,xpw,xpwd, &
535 W1l,W2l,W3l,W4l,W1s,W2s,W3s,W4s, &
536 w1,w2,w3,w4,qrsk(its:ite,kts:kte),evef,ptem,ptem1
538 logical :: totflg, cnvflg(its:ite),flg(its:ite),lclflg
541 ! climatological critical cloud work functions for closure
543 data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
544 350.,300.,250.,200.,150./
545 data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
546 .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
548 !-----------------------------------------------------------------------
550 ! define miscellaneous values
557 el2orc = hvap_*hvap_/(rv_*cp_)
559 fact1 = (cvap_-cliq_)/rv_
560 fact2 = hvap_/rv_-fact1*t0c_
564 dtmin = max(dt2,1200.)
565 dtmax = max(dt2,3600.)
567 if (dx_factor_nsas == 1) then
568 dx_factor = 250. / delx ! assume 2.5 ms-1 (1km) and 1.125 cms-1 (200km)
569 W1l = dx_factor * 0.1 * (-1.)
570 W2l = dx_factor * (-1.)
571 W3l = dx_factor * (-1.)
572 W4l = dx_factor * 0.1 * (-1.)
613 acrit(k) = acritt(k) * (975. - pcrit(k))
616 ! Define top layer for search of the downdraft originating layer
617 ! and the maximum thetae for updraft
624 if(prsl(i,k).gt.prsi(i,1)*0.45) kbmax = k + 1
625 if(prsl(i,k).gt.prsi(i,1)*0.70) kbm = k + 1
626 if(prsl(i,k).gt.prsi(i,1)*0.04) kmax = k + 1
633 ! convert surface pressure to mb from cb
649 p(i,k) = prsl(i,k) * 10.
662 uo(i,k) = u1(i,k) * rcs
663 vo(i,k) = v1(i,k) * rcs
669 ! p is pressure of the layer (mb)
670 ! t is temperature at t-dt (k)..tn
671 ! q is mixing ratio at t-dt (kg/kg)..qn
672 ! to is temperature at t+dt (k)... this is after advection and turbulan
673 ! qo is mixing ratio at t+dt (kg/kg)..q1
677 qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
678 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
679 qeso(i,k) = max(qeso(i,k),qmin_)
680 qo(i,k) = max(qo(i,k), 1.e-10 )
681 tvo(i,k) = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
685 ! compute moist static energy
689 heo(i,k) = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qo(i,k)
690 heso(i,k) = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qeso(i,k)
694 ! Determine level with largest moist static energy
695 ! This is the level where updraft starts
704 if(heo(i,k).gt.hmax(i)) then
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
770 if(kbcon(i).eq.kmax) cnvflg(i) = .false.
775 totflg = totflg .and. (.not. cnvflg(i))
782 ! determine critical convective inhibition
783 ! as a function of vertical velocity at cloud base.
785 pdot(i) = 10.* dot(i,kbcon(i))
786 if(slimsk(i).eq.1.) then
797 if(pdot(i).le.w4) then
798 tem = (pdot(i) - w4) / (w3 - w4)
799 elseif(pdot(i).ge.-w4) then
800 tem = - (pdot(i) + w4) / (w4 - w3)
807 tem1= .5*(cincrmax-cincrmin)
808 cincr = cincrmax - tem * tem1
809 pbcdif(i) = -p(i,kbcon(i)) + p(i,kb(i))
810 if(pbcdif(i).gt.cincr) cnvflg(i) = .false.
817 totflg = totflg .and. (.not. cnvflg(i))
823 zi(i,k) = 0.5*(zl(i,k-1)+zl(i,k))
829 xlamb(i,k) = clam / zi(i,k+1)
833 ! assume that updraft entrainment rate above cloud base is
834 ! same as that at cloud base
838 if(cnvflg(i).and.(k.gt.kbcon(i))) then
839 xlamb(i,k) = xlamb(i,kbcon(i))
844 ! assume the detrainment rate for the updrafts to be same as
845 ! the entrainment rate at cloud base
849 xlamud(i) = xlamb(i,kbcon(i))
853 ! functions rapidly decreasing with height, mimicking a cloud ensemble
854 ! (Bechtold et al., 2008)
858 if(cnvflg(i).and.(k.gt.kbcon(i))) then
859 tem = qeso(i,k)/qeso(i,kbcon(i))
866 ! final entrainment rate as the sum of turbulent part and organized entrainment
867 ! depending on the environmental relative humidity
868 ! (Bechtold et al., 2008)
872 if(cnvflg(i).and.(k.ge.kbcon(i))) then
873 tem = cxlamu * frh(i,k) * fent2(i,k)
874 xlamb(i,k) = xlamb(i,k)*fent1(i,k) + tem
879 ! determine updraft mass flux
891 if(cnvflg(i).and.k.lt.kbcon(i).and.k.ge.kb(i)) then
892 dz = zi(i,k+2) - zi(i,k+1)
893 ptem = 0.5*(xlamb(i,k)+xlamb(i,k+1))-xlamud(i)
894 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
900 if(cnvflg(i).and.k.gt.kbcon(i)) then
901 dz = zi(i,k+1) - zi(i,k)
902 ptem = 0.5*(xlamb(i,k)+xlamb(i,k-1))-xlamud(i)
903 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
909 dz = zi(i,3) - zi(i,2)
910 ptem = 0.5*(xlamb(i,1)+xlamb(i,2))-xlamud(i)
911 eta(i,1) = eta(i,2) / (1. + ptem * dz)
915 ! work up updraft cloud properties
920 hcko(i,indx) = hkbo(i)
921 qcko(i,indx) = qkbo(i)
922 ucko(i,indx) = uo(i,indx)
923 vcko(i,indx) = vo(i,indx)
928 ! cloud property below cloud base is modified by the entrainment proces
932 if(cnvflg(i).and.k.gt.kb(i)) then
933 dz = zi(i,k+1) - zi(i,k)
934 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
935 tem1 = 0.5 * xlamud(i) * dz
936 factor = 1. + tem - tem1
937 ptem = 0.5 * tem + pgcon
938 ptem1= 0.5 * tem - pgcon
939 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
940 (heo(i,k)+heo(i,k-1)))/factor
941 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
942 +ptem1*uo(i,k-1))/factor
943 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
944 +ptem1*vo(i,k-1))/factor
945 dbyo(i,k) = hcko(i,k) - heso(i,k)
950 ! taking account into convection inhibition due to existence of
951 ! dry layers below cloud base
960 if(flg(i).and.k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
969 if(kbcon1(i).eq.kmax) cnvflg(i) = .false.
975 tem = p(i,kbcon(i)) - p(i,kbcon1(i))
984 totflg = totflg .and. (.not. cnvflg(i))
989 ! determine cloud top
1001 if(dbyo(i,k).lt.0..and.flg(i).and.k.gt. kbcon1(i)) then
1012 if(cnvflg(i).and.(p(i,kbcon(i)) - p(i,ktcon(i))).lt.150.) &
1018 totflg = totflg .and. (.not. cnvflg(i))
1023 ! search for downdraft originating level above theta-e minimum
1027 hmin(i) = heo(i,kbcon1(i))
1035 if(cnvflg(i).and.k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
1042 ! make sure that jmin is within the cloud
1046 jmin(i) = min(lmin(i),ktcon(i)-1)
1047 jmin(i) = max(jmin(i),kbcon1(i)+1)
1048 if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
1052 ! specify upper limit of mass flux at cloud base
1057 dp = 1000. * del(i,k)
1058 xmbmax(i) = dp / (g_ * dt2)
1063 ! compute cloud moisture property and precipitation
1071 if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1072 dz = .5 * (zl(i,k+1) - zl(i,k-1))
1073 dz1 = (zi(i,k+1) - zi(i,k))
1074 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1076 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1077 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz1
1078 tem1 = 0.5 * xlamud(i) * dz1
1079 factor = 1. + tem - tem1
1080 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
1081 (qo(i,k)+qo(i,k-1)))/factor
1082 qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1084 ! check if there is excess moisture to release latent heat
1086 if(qcirs(i,k).gt.0. .and. k.ge.kbcon(i)) then
1087 etah = .5 * (eta(i,k) + eta(i,k-1))
1088 if(ncloud.gt.0..and.k.gt.jmin(i)) then
1089 dp = 1000. * del(i,k)
1090 qlk = qcirs(i,k) / (eta(i,k) + etah * (c0 + c1) * dz1)
1091 dellal(i,k) = etah * c1 * dz1 * qlk * g_ / dp
1093 qlk = qcirs(i,k) / (eta(i,k) + etah * c0 * dz1)
1095 aa1(i) = aa1(i) - dz1 * g_ * qlk
1097 pwo(i,k) = etah * c0 * dz1 * qlk
1099 pwavo(i) = pwavo(i) + pwo(i,k)
1105 ! calculate cloud work function at t+dt
1109 if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1110 dz1 = zl(i,k+1) - zl(i,k)
1111 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1112 rfact = 1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1113 aa1(i) = aa1(i) +dz1 * (g_ / (cp_ * to(i,k))) &
1114 * dbyo(i,k) / (1. + gamma)* rfact
1115 aa1(i) = aa1(i)+dz1 * g_ * fv_ * &
1116 max(0.,(qeso(i,k) - qo(i,k)))
1122 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1127 totflg = totflg .and. (.not. cnvflg(i))
1131 ! estimate the convective overshooting as the level
1132 ! where the [aafac * cloud work function] becomes zero,
1133 ! which is the final cloud top
1137 aa2(i) = aafac * aa1(i)
1149 if(k.ge.ktcon(i).and.k.lt.kmax) then
1150 dz1 = zl(i,k+1) - zl(i,k)
1151 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1152 rfact = 1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1153 aa2(i) = aa2(i) +dz1 * (g_ / (cp_ * to(i,k))) &
1154 * dbyo(i,k) / (1. + gamma)* rfact
1155 if(aa2(i).lt.0.) then
1164 ! compute cloud moisture property, detraining cloud water
1165 ! and precipitation in overshooting layers
1170 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
1171 dz = (zi(i,k+1) - zi(i,k))
1172 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1173 qrch = qeso(i,k)+ gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1174 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1175 tem1 = 0.5 * xlamud(i) * dz
1176 factor = 1. + tem - tem1
1177 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
1178 (qo(i,k)+qo(i,k-1)))/factor
1179 qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1181 ! check if there is excess moisture to release latent heat
1183 if(qcirs(i,k).gt.0.) then
1184 etah = .5 * (eta(i,k) + eta(i,k-1))
1185 if(ncloud.gt.0.) then
1186 dp = 1000. * del(i,k)
1187 qlk = qcirs(i,k) / (eta(i,k) + etah * (c0 + c1) * dz)
1188 dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
1190 qlk = qcirs(i,k) / (eta(i,k) + etah * c0 * dz)
1193 pwo(i,k) = etah * c0 * dz * qlk
1195 pwavo(i) = pwavo(i) + pwo(i,k)
1202 ! exchange ktcon with ktcon1
1207 ktcon(i) = ktcon1(i)
1212 ! this section is ready for cloud water
1214 if (ncloud.gt.0) then
1216 ! compute liquid and vapor separation at cloud top
1221 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1223 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1224 dq = qcko(i,k) - qrch
1226 ! check if there is excess moisture to release latent heat
1236 ! ..... downdraft calculations .....
1238 ! determine downdraft strength in terms of wind shear
1248 if(k.gt.kb(i).and.k.le.ktcon(i).and.cnvflg(i)) then
1249 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
1250 + (vo(i,k)-vo(i,k-1)) ** 2)
1251 vshear(i) = vshear(i) + shear
1258 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i)+1)-zi(i,kb(i)+1))
1259 e1 = 1.591-.639*vshear(i) &
1260 +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1263 edt(i) = min(edt(i),.9)
1264 edt(i) = max(edt(i),.0)
1270 ! determine detrainment rate between 1 and kbdtr
1280 if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
1281 dz = zi(i,k+2) - zi(i,k+1)
1282 sumx(i) = sumx(i) + dz
1290 if(slimsk(i).eq.1.) beta = betal
1293 kbdtr(i) = max(kbdtr(i),1)
1294 dz =(sumx(i)+zi(i,2))/float(kbcon(i))
1295 tem = 1./float(kbcon(i))
1296 xlamd(i) = (1.-beta**tem)/dz
1300 ! determine downdraft mass flux
1315 if(k.lt.jmin(i).and.k.ge.kbcon(i))then
1316 dz = (zi(i,k+2) - zi(i,k+1))
1317 ptem = xlamdd-xlamde
1318 etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1319 elseif(k.lt.kbcon(i))then
1320 dz = (zi(i,k+2) - zi(i,k+1))
1321 ptem = xlamd(i)+xlamdd-xlamde
1322 etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1329 ! downdraft moisture properties
1340 hcdo(i,jmn) = heo(i,jmn)
1341 qcdo(i,jmn) = qo(i,jmn)
1342 qrcdo(i,jmn) = qeso(i,jmn)
1343 ucdo(i,jmn) = uo(i,jmn)
1344 vcdo(i,jmn) = vo(i,jmn)
1350 if (cnvflg(i) .and. k.lt.jmin(i)) then
1351 dz = zi(i,k+2) - zi(i,k+1)
1352 if(k.ge.kbcon(i)) then
1354 tem1 = 0.5 * xlamdd * dz
1357 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1359 factor = 1. + tem - tem1
1360 ptem = 0.5 * tem - pgcon
1361 ptem1= 0.5 * tem + pgcon
1362 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
1363 (heo(i,k)+heo(i,k+1)))/factor
1364 ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1) &
1365 +ptem1*uo(i,k))/factor
1366 vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1) &
1367 +ptem1*vo(i,k))/factor
1368 dbyo(i,k) = hcdo(i,k) - heso(i,k)
1375 if(cnvflg(i).and.k.lt.jmin(i)) then
1378 gamma = el2orc * dq / dt**2
1379 qrcdo(i,k)=dq+(1./hvap_)*(gamma/(1.+gamma))*dbyo(i,k)
1380 detad = etad(i,k+1) - etad(i,k)
1381 dz = zi(i,k+2) - zi(i,k+1)
1382 if(k.ge.kbcon(i)) then
1384 tem1 = 0.5 * xlamdd * dz
1387 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1389 factor = 1. + tem - tem1
1390 qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
1391 (qo(i,k)+qo(i,k+1)))/factor
1392 pwdo(i,k) = etad(i,k+1) * qcdo(i,k) -etad(i,k+1) * qrcdo(i,k)
1393 qcdo(i,k) = qrcdo(i,k)
1394 pwevo(i) = pwevo(i) + pwdo(i,k)
1399 ! final downdraft strength dependent on precip
1400 ! efficiency (edt), normalized condensate (pwav), and
1405 if(slimsk(i).eq.2.) edtmax = edtmaxs
1407 if(pwevo(i).lt.0.) then
1408 edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1409 edto(i) = min(edto(i),edtmax)
1416 ! downdraft cloudwork functions
1420 if(cnvflg(i).and.k.lt.jmin(i)) then
1421 gamma = el2orc * qeso(i,k) / to(i,k)**2
1426 dz=-1.*(zl(i,k+1)-zl(i,k))
1427 aa1(i)=aa1(i)+edto(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg)) &
1428 *(1.+fv_*cp_*dg*dt/hvap_)
1429 aa1(i)=aa1(i)+edto(i)*dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1435 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1440 totflg = totflg .and. (.not. cnvflg(i))
1444 ! what would the change be, that a cloud with unit mass
1445 ! will do to the environment?
1460 dp = 1000. * del(i,1)
1461 dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1) &
1462 - heo(i,1)) * g_ / dp
1463 dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1) &
1464 - qo(i,1)) * g_ / dp
1465 dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1) &
1466 - uo(i,1)) * g_ / dp
1467 dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1) &
1468 - vo(i,1)) * g_ / dp
1472 ! changed due to subsidence and entrainment
1476 if(cnvflg(i).and.k.lt.ktcon(i)) then
1478 if(k.le.kb(i)) aup = 0.
1480 if(k.gt.jmin(i)) adw = 0.
1482 dv2 = .5 * (heo(i,k) + heo(i,k-1))
1485 dv2q = .5 * (qo(i,k) + qo(i,k-1))
1488 dv2u = .5 * (uo(i,k) + uo(i,k-1))
1491 dv2v = .5 * (vo(i,k) + vo(i,k-1))
1493 dp = 1000. * del(i,k)
1494 dz = zi(i,k+1) - zi(i,k)
1495 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1))
1497 if(k.le.kbcon(i)) then
1499 ptem1 = xlamd(i)+xlamdd
1504 deta = eta(i,k) - eta(i,k-1)
1505 detad = etad(i,k) - etad(i,k-1)
1506 dellah(i,k) = dellah(i,k) + &
1507 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1 &
1508 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3 &
1509 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2*dz &
1510 + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
1511 + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz) *g_/dp
1512 dellaq(i,k) = dellaq(i,k) + &
1513 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1q &
1514 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3q &
1515 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz &
1516 + aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
1517 + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1))*dz) *g_/dp
1518 dellau(i,k) = dellau(i,k) + &
1519 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1u &
1520 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3u &
1521 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz &
1522 + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
1523 + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz &
1524 - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u))*g_/dp
1526 dellav(i,k) = dellav(i,k) + &
1527 ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1v &
1528 - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3v &
1529 - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz &
1530 + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
1531 + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz &
1532 - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v))*g_/dp
1542 dp = 1000. * del(i,indx)
1544 dellah(i,indx) = eta(i,indx-1) * &
1545 (hcko(i,indx-1) - dv1) * g_ / dp
1547 dellaq(i,indx) = eta(i,indx-1) * &
1548 (qcko(i,indx-1) - dvq1) * g_ / dp
1550 dellau(i,indx) = eta(i,indx-1) * &
1551 (ucko(i,indx-1) - dv1u) * g_ / dp
1553 dellav(i,indx) = eta(i,indx-1) * &
1554 (vcko(i,indx-1) - dv1v) * g_ / dp
1558 dellal(i,indx) = eta(i,indx-1) * qlko_ktcon(i) * g_ / dp
1562 ! final changed variable per unit mass flux
1566 if(cnvflg(i).and.k.gt.ktcon(i)) then
1570 if(cnvflg(i).and.k.le.ktcon(i)) then
1571 qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
1572 dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1573 to(i,k) = dellat * mbdt + t1(i,k)
1574 qo(i,k) = max(qo(i,k),1.0e-10)
1579 !------------------------------------------------------------------------------
1581 ! the above changed environment is now used to calulate the
1582 ! effect the arbitrary cloud (with unit mass flux)
1583 ! which then is used to calculate the real mass flux,
1584 ! necessary to keep this change in balance with the large-scale
1587 ! environmental conditions again
1589 !------------------------------------------------------------------------------
1594 qeso(i,k)=0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1595 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1596 qeso(i,k) = max(qeso(i,k),qmin_)
1597 tvo(i,k) = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
1609 ! moist static energy
1614 dz = .5 * (zl(i,k+1) - zl(i,k))
1615 dp = .5 * (p(i,k+1) - p(i,k))
1616 es =0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1617 pprime = p(i,k+1) + (eps-1.) * es
1618 qs = eps * es / pprime
1619 dqsdp = - qs / pprime
1620 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1621 dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
1622 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1623 dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
1624 dq = dqsdt * dt + dqsdp * dp
1625 to(i,k) = to(i,k+1) + dt
1626 qo(i,k) = qo(i,k+1) + dq
1627 po = .5 * (p(i,k) + p(i,k+1))
1628 qeso(i,k) =0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1629 qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
1630 qeso(i,k) = max(qeso(i,k),qmin_)
1631 qo(i,k) = max(qo(i,k), 1.0e-10)
1632 heo(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
1633 cp_ * to(i,k) + hvap_ * qo(i,k)
1634 heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
1635 cp_ * to(i,k) + hvap_ * qeso(i,k)
1643 heo(i,k) = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qo(i,k)
1644 heso(i,k) = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qeso(i,k)
1653 xhkb(i) = heo(i,indx)
1654 xqkb(i) = qo(i,indx)
1655 hcko(i,indx) = xhkb(i)
1656 qcko(i,indx) = xqkb(i)
1660 ! ..... static control .....
1662 ! moisture and cloud work functions
1666 if(cnvflg(i).and.k.gt.kb(i).and.k.le.ktcon(i)) then
1667 dz = zi(i,k+1) - zi(i,k)
1668 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1669 tem1 = 0.5 * xlamud(i) * dz
1670 factor = 1. + tem - tem1
1671 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
1672 (heo(i,k)+heo(i,k-1)))/factor
1679 if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1680 dz = zi(i,k+1) - zi(i,k)
1681 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1682 xdby = hcko(i,k) - heso(i,k)
1684 + gamma * xdby / (hvap_ * (1. + gamma))
1685 tem = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1686 tem1 = 0.5 * xlamud(i) * dz
1687 factor = 1. + tem - tem1
1688 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*(qo(i,k)+qo(i,k-1)))/factor
1689 dq = eta(i,k) * qcko(i,k) - eta(i,k) * xqrch
1690 if(k.ge.kbcon(i).and.dq.gt.0.) then
1691 etah = .5 * (eta(i,k) + eta(i,k-1))
1692 if(ncloud.gt.0..and.k.gt.jmin(i)) then
1693 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
1695 qlk = dq / (eta(i,k) + etah * c0 * dz)
1697 if(k.lt.ktcon1(i)) then
1698 xaa0(i) = xaa0(i) - dz * g_ * qlk
1700 qcko(i,k) = qlk + xqrch
1701 xpw = etah * c0 * dz * qlk
1702 xpwav(i) = xpwav(i) + xpw
1705 if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
1706 dz1 = zl(i,k+1) - zl(i,k)
1707 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1708 rfact = 1. + fv_ * cp_ * gamma &
1710 xdby = hcko(i,k) - heso(i,k)
1712 + dz1 * (g_ / (cp_ * to(i,k))) &
1713 * xdby / (1. + gamma) &
1717 max(0.,(qeso(i,k) - qo(i,k)))
1722 ! ..... downdraft calculations .....
1725 ! downdraft moisture properties
1734 xhcd(i,jmn) = heo(i,jmn)
1735 xqcd(i,jmn) = qo(i,jmn)
1736 qrcd(i,jmn) = qeso(i,jmn)
1742 if(cnvflg(i).and.k.lt.jmin(i)) then
1743 dz = zi(i,k+2) - zi(i,k+1)
1744 if(k.ge.kbcon(i)) then
1746 tem1 = 0.5 * xlamdd * dz
1749 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1751 factor = 1. + tem - tem1
1752 xhcd(i,k) = ((1.-tem1)*xhcd(i,k+1)+tem*0.5* &
1753 (heo(i,k)+heo(i,k+1)))/factor
1760 if(cnvflg(i).and.k.lt.jmin(i)) then
1763 gamma = el2orc * dq / dt**2
1764 dh = xhcd(i,k) - heso(i,k)
1765 qrcd(i,k)=dq+(1./hvap_)*(gamma/(1.+gamma))*dh
1766 dz = zi(i,k+2) - zi(i,k+1)
1767 if(k.ge.kbcon(i)) then
1769 tem1 = 0.5 * xlamdd * dz
1772 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1774 factor = 1. + tem - tem1
1775 xqcd(i,k) = ((1.-tem1)*xqcd(i,k+1)+tem*0.5* &
1776 (qo(i,k)+qo(i,k+1)))/factor
1777 xpwd = etad(i,k+1) * (xqcd(i,k) - qrcd(i,k))
1778 xqcd(i,k)= qrcd(i,k)
1779 xpwev(i) = xpwev(i) + xpwd
1786 if(slimsk(i).eq.2.) edtmax = edtmaxs
1788 if(xpwev(i).ge.0.) then
1791 edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1792 edtx(i) = min(edtx(i),edtmax)
1797 ! downdraft cloudwork functions
1801 if(cnvflg(i).and.k.lt.jmin(i)) then
1802 gamma = el2orc * qeso(i,k) / to(i,k)**2
1807 dz=-1.*(zl(i,k+1)-zl(i,k))
1808 xaa0(i)=xaa0(i)+edtx(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg)) &
1809 *(1.+fv_*cp_*dg*dt/hvap_)
1810 xaa0(i)=xaa0(i)+edtx(i)* &
1811 dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1816 ! calculate critical cloud work function
1821 if(p(i,ktcon(i)).lt.pcrit(15))then
1822 acrt(i)=acrit(15)*(975.-p(i,ktcon(i)))/(975.-pcrit(15))
1823 else if(p(i,ktcon(i)).gt.pcrit(1))then
1826 k = int((850. - p(i,ktcon(i)))/50.) + 2
1829 acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))* &
1830 (p(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
1841 if(slimsk(i).eq.1.) then
1848 if(pdot(i).le.w4) then
1849 acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
1850 elseif(pdot(i).ge.-w4) then
1851 acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
1855 acrtfct(i) = max(acrtfct(i),-1.)
1856 acrtfct(i) = min(acrtfct(i),1.)
1857 acrtfct(i) = 1. - acrtfct(i)
1858 dtconv(i) = dt2 + max((1800. - dt2),0.) * (pdot(i) - w2) / (w1 - w2)
1859 dtconv(i) = max(dtconv(i),dtmin)
1860 dtconv(i) = min(dtconv(i),dtmax)
1865 ! large scale forcing
1869 f(i) = (aa1(i) - acrt(i) * acrtfct(i)) / dtconv(i)
1870 if(f(i).le.0.) cnvflg(i) = .false.
1873 xk(i) = (xaa0(i) - aa1(i)) / mbdt
1874 if(xk(i).ge.0.) cnvflg(i) = .false.
1877 ! kernel, cloud base mass flux
1880 xmb(i) = -f(i) / xk(i)
1881 xmb(i) = min(xmb(i),xmbmax(i))
1890 totflg = totflg .and. (.not. cnvflg(i))
1894 ! restore t0 and qo to t1 and q1 in case convection stops
1903 qeso(i,k) = 0.01*fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1904 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1905 qeso(i,k) = max(qeso(i,k),qmin_)
1910 ! feedback: simply the changes from the cloud with unit mass flux
1911 ! multiplied by the mass flux necessary to keep the
1912 ! equilibrium with the larger-scale.
1926 if(cnvflg(i).and.k.le.ktcon(i)) then
1928 if(k.le.kb(i)) aup = 0.
1930 if(k.gt.jmin(i)) adw = 0.
1931 dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1932 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1933 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1935 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1936 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1937 dp = 1000. * del(i,k)
1938 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
1939 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
1940 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
1941 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
1942 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
1954 if (cnvflg(i) .and. k.le.ktcon(i)) then
1955 qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1956 qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
1957 qeso(i,k) = max(qeso(i,k), qmin_ )
1973 if(cnvflg(i).and.k.lt.ktcon(i)) then
1975 if(k.le.kb(i)) aup = 0.
1977 if(k.ge.jmin(i)) adw = 0.
1978 rntot(i) = rntot(i) &
1979 + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)) &
1980 * xmb(i) * .001 * dt2
1985 ! conversion rainfall (m) and compute the evaporation of falling raindrops
1992 if(cnvflg(i).and.k.lt.ktcon(i)) then
1994 if(k.le.kb(i)) aup = 0.
1996 if(k.ge.jmin(i)) adw = 0.
1998 + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)) &
1999 * xmb(i) * .001 * dt2
2001 if(cnvflg(i).and.flg(i).and.k.lt.ktcon(i)) then
2003 evef = edt(i) * evfacts
2004 if(slimsk(i).eq.1.) evef = edt(i) * evfactl
2005 qcond(i) = evef * (q1(i,k) - qeso(i,k)) / (1. + el2orc * &
2006 qeso(i,k) / t1(i,k)**2)
2007 dp = 1000. * del(i,k)
2008 if(rain(i).gt.0..and.qcond(i).lt.0.) then
2009 qevap(i) = -qcond(i) * (1. - exp(-.32 * sqrt(dt2 * rain(i))))
2010 qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
2011 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
2012 if (delq2(i).gt.rntot(i)) then
2013 qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
2017 if(rain(i).gt.0..and.qevap(i).gt.0.) then
2018 q1(i,k) = q1(i,k) + qevap(i)
2019 t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
2020 rain(i) = rain(i) - .001 * qevap(i) * dp / g_
2021 delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
2022 deltv(i) = - (hvap_/cp_)*qevap(i)/dt2
2023 delq(i) = + qevap(i)/dt2
2025 dellaq(i,k) = dellaq(i,k) + delq(i)/xmb(i)
2026 delqbar(i) = delqbar(i) + delq(i)*dp/g_
2027 deltbar(i) = deltbar(i) + deltv(i)*dp/g_
2033 ! consider the negative rain in the event of rain evaporation and downdrafts
2037 if(rain(i).lt.0..and..not.flg(i)) rain(i) = 0.
2038 if(rain(i).le.0.) then
2050 if(cnvflg(i).and.rain(i).le.0.) then
2059 ! detrainment of cloud water and ice
2061 if (ncloud.gt.0) then
2064 if (cnvflg(i) .and. rain(i).gt.0.) then
2065 if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
2066 tem = dellal(i,k) * xmb(i) * dt2
2067 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2068 if (ncloud.ge.2) then
2069 qi2(i,k) = qi2(i,k) + tem * tem1 ! ice
2070 qc2(i,k) = qc2(i,k) + tem *(1.0-tem1) ! water
2072 qc2(i,k) = qc2(i,k) + tem
2080 end subroutine nsas2d
2081 !-------------------------------------------------------------------------------
2082 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
2083 !-------------------------------------------------------------------------------
2085 !-------------------------------------------------------------------------------
2086 REAL :: t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
2096 xbi=xai+hsub/(rv*ttp)
2098 if(t.lt.ttp.and.ice.eq.1) then
2099 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2101 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2106 if(t.lt.ttp.and.ice.eq.1) then
2107 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2109 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2115 if(t.lt.ttp.and.ice.eq.1) then
2116 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2118 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2123 !-------------------------------------------------------------------------------
2125 !-------------------------------------------------------------------------------
2126 subroutine nsasinit(rthcuten,rqvcuten,rqccuten,rqicuten, &
2128 restart,p_qc,p_qi,p_first_scalar, &
2130 ids, ide, jds, jde, kds, kde, &
2131 ims, ime, jms, jme, kms, kme, &
2132 its, ite, jts, jte, kts, kte )
2133 !-------------------------------------------------------------------------------
2135 !-------------------------------------------------------------------------------
2136 logical , intent(in) :: allowed_to_read,restart
2137 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
2138 ims, ime, jms, jme, kms, kme, &
2139 its, ite, jts, jte, kts, kte
2140 integer , intent(in) :: p_first_scalar, p_qi, p_qc
2141 real, dimension( ims:ime , kms:kme , jms:jme ) , intent(out) :: &
2148 integer :: i, j, k, itf, jtf, ktf
2154 if(.not.restart)then
2166 if (p_qc .ge. p_first_scalar) then
2176 if (p_qi .ge. p_first_scalar) then
2187 end subroutine nsasinit
2189 !-------------------------------------------------------------------------------
2190 ! NCEP SCV (Shallow Convection Scheme)
2191 !-------------------------------------------------------------------------------
2192 subroutine nscv2d(delt,del,prsl,prsi,prslk,zl, &
2193 ncloud,qc2,qi2,q1,t1,rain,kbot,ktop, &
2196 cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2, &
2200 ids,ide, jds,jde, kds,kde, &
2201 ims,ime, jms,jme, kms,kme, &
2202 its,ite, jts,jte, kts,kte)
2203 !-------------------------------------------------------------------------------
2204 ! subprogram: nscv2d computes shallow-convective heating and moisng
2206 ! abstract: computes non-precipitating convective heating and moistening
2207 ! using a one cloud type arakawa-schubert convection scheme as in the ncep
2208 ! sas scheme. the scheme has been operational at ncep gfs model since july 2010
2209 ! the scheme includes updraft and downdraft effects, but the cloud depth is
2210 ! limited less than 150 hpa.
2212 ! developed by jong-il han and hua-lu pan
2213 ! implemented into wrf by jiheon jang and songyou hong
2214 ! module with cpp-based options is available in grims
2216 ! program history log:
2217 ! 10-07-01 jong-il han initial operational implementation at ncep gfs
2218 ! 10-12-01 jihyeon jang implemented into wrf
2220 ! subprograms called:
2221 ! fpvs - function to compute saturation vapor pressure
2224 ! han and pan (2010, wea. forecasting)
2226 !-------------------------------------------------------------------------------
2228 !-------------------------------------------------------------------------------
2232 integer :: ids,ide, jds,jde, kds,kde, &
2233 ims,ime, jms,jme, kms,kme, &
2234 its,ite, jts,jte, kts,kte
2235 real :: cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
2236 real :: pi_,qmin_,t0c_
2237 real :: cice,xlv0,xls,psat
2240 real :: del(its:ite,kts:kte), &
2241 prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme), &
2242 prsi(its:ite,kts:kte+1),zl(its:ite,kts:kte)
2244 real :: slimsk(ims:ime)
2245 real :: dot(its:ite,kts:kte)
2246 real :: hpbl(ims:ime)
2248 real :: hfx(ims:ime),qfx(ims:ime)
2250 real :: qi2(its:ite,kts:kte),qc2(its:ite,kts:kte)
2251 real :: q1(its:ite,kts:kte), &
2252 t1(its:ite,kts:kte), &
2253 u1(its:ite,kts:kte), &
2255 integer :: icps(its:ite)
2257 real :: rain(its:ite)
2258 integer :: kbot(its:ite),ktop(its:ite)
2260 ! local variables and arrays
2262 integer :: i,j,indx, jmn, k, kk, km1
2263 integer :: kpbl(its:ite)
2266 desdt, deta, detad, dg, &
2267 dh, dhh, dlnsig, dp, &
2268 dq, dqsdp, dqsdt, dt, &
2269 dt2, dtmax, dtmin, &
2274 dz, dz1, e1, clam, &
2277 evef, evfact, evfactl, &
2279 gamma, pprime, betaw, &
2281 rfact, shear, tem1, &
2283 val2, w1, w1l, w1s, &
2285 w3l, w3s, w4, w4l, &
2286 w4s, tem, ptem, ptem1, &
2289 integer :: kb(its:ite), kbcon(its:ite), kbcon1(its:ite), &
2290 ktcon(its:ite), ktcon1(its:ite), &
2291 kbm(its:ite), kmax(its:ite)
2293 real :: aa1(its:ite), &
2294 delhbar(its:ite), delq(its:ite), &
2295 delq2(its:ite), delqev(its:ite), rntot(its:ite), &
2296 delqbar(its:ite), deltbar(its:ite), &
2297 deltv(its:ite), edt(its:ite), &
2298 wstar(its:ite), sflx(its:ite), &
2299 pdot(its:ite), po(its:ite,kts:kte), &
2300 qcond(its:ite), qevap(its:ite), hmax(its:ite), &
2302 xlamud(its:ite), xmb(its:ite), xmbmax(its:ite)
2303 real :: delubar(its:ite), delvbar(its:ite)
2307 real :: thx(its:ite, kts:kte)
2308 real :: rhox(its:ite)
2311 real :: p(its:ite,kts:kte), to(its:ite,kts:kte), &
2312 qo(its:ite,kts:kte), qeso(its:ite,kts:kte), &
2313 uo(its:ite,kts:kte), vo(its:ite,kts:kte)
2317 real :: qlko_ktcon(its:ite), dellal(its:ite,kts:kte), &
2318 dbyo(its:ite,kts:kte), &
2319 xlamue(its:ite,kts:kte), &
2320 heo(its:ite,kts:kte), heso(its:ite,kts:kte), &
2321 dellah(its:ite,kts:kte), dellaq(its:ite,kts:kte), &
2322 dellau(its:ite,kts:kte), dellav(its:ite,kts:kte), &
2323 ucko(its:ite,kts:kte), vcko(its:ite,kts:kte), &
2324 hcko(its:ite,kts:kte), qcko(its:ite,kts:kte), &
2325 eta(its:ite,kts:kte), zi(its:ite,kts:kte+1), &
2326 pwo(its:ite,kts:kte)
2328 logical :: totflg, cnvflg(its:ite), flg(its:ite)
2330 ! physical parameters
2332 real,parameter :: c0=.002,c1=5.e-4
2333 real,parameter :: cincrmax=180.,cincrmin=120.,dthk=25.
2334 real :: el2orc,fact1,fact2,eps
2335 real,parameter :: h1=0.33333333
2336 real,parameter :: tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)
2337 !-------------------------------------------------------------------------------
2344 ! compute surface buoyancy flux
2348 thx(i,k) = t1(i,k)/prslk(i,k)
2353 tvcon = (1.+fv_*q1(i,1))
2354 rhox(i) = prsl(i,1)*1.e3/(rd_*t1(i,1)*tvcon)
2358 ! sflx(i) = heat(i)+fv_*t1(i,1)*evap(i)
2359 sflx(i) = hfx(i)/rhox(i)/cp_ + qfx(i)/rhox(i)*fv_*thx(i,1)
2366 if(icps(i).eq.1) cnvflg(i) = .false.
2367 if(sflx(i).le.0.) cnvflg(i) = .false.
2385 totflg = totflg .and. (.not. cnvflg(i))
2391 dtmin = max(dt2, val )
2393 dtmax = max(dt2, val )
2395 ! model tunable parameters are all here
2404 ! define miscellaneous values
2406 el2orc = hvap_*hvap_/(rv_*cp_)
2408 fact1 = (cvap_-cliq_)/rv_
2409 fact2 = hvap_/rv_-fact1*t0c_
2420 ! define top layer for search of the downdraft originating layer
2421 ! and the maximum thetae for updraft
2430 if (prsl(i,k).gt.prsi(i,1)*0.70) kbm(i) = k + 1
2431 if (prsl(i,k).gt.prsi(i,1)*0.60) kmax(i) = k + 1
2436 kbm(i) = min(kbm(i),kmax(i))
2439 ! hydrostatic height assume zero terr and compute
2440 ! updraft entrainment rate as an inverse function of height
2444 zi(i,k) = 0.5*(zl(i,k-1)+zl(i,k))
2450 xlamue(i,k) = clam / zi(i,k+1)
2455 xlamue(i,kte) = xlamue(i,km1)
2467 if (flg(i).and.zl(i,k).le.hpbl(i)) then
2476 kpbl(i)= min(kpbl(i),kbm(i))
2479 ! convert surface pressure to mb from cb
2484 if (cnvflg(i) .and. k .le. kmax(i)) then
2485 p(i,k) = prsl(i,k) * 10.0
2496 uo(i,k) = u1(i,k) * rcs
2497 vo(i,k) = v1(i,k) * rcs
2504 ! p is pressure of the layer (mb)
2505 ! t is temperature at t-dt (k)..tn
2506 ! q is mixing ratio at t-dt (kg/kg)..qn
2507 ! to is temperature at t+dt (k)... this is after advection and turbulan
2508 ! qo is mixing ratio at t+dt (kg/kg)..q1
2512 if (cnvflg(i) .and. k .le. kmax(i)) then
2513 qeso(i,k) = 0.01 * fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2514 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
2516 qeso(i,k) = max(qeso(i,k), val1)
2518 qo(i,k) = max(qo(i,k), val2 )
2523 ! compute moist static energy
2527 if (cnvflg(i) .and. k .le. kmax(i)) then
2528 tem = g_ * zl(i,k) + cp_ * to(i,k)
2529 heo(i,k) = tem + hvap_ * qo(i,k)
2530 heso(i,k) = tem + hvap_ * qeso(i,k)
2535 ! determine level with largest moist static energy within pbl
2536 ! this is the level where updraft starts
2547 if (cnvflg(i).and.k.le.kpbl(i)) then
2548 if(heo(i,k).gt.hmax(i)) then
2558 if (cnvflg(i) .and. k .le. kmax(i)-1) then
2559 dz = .5 * (zl(i,k+1) - zl(i,k))
2560 dp = .5 * (p(i,k+1) - p(i,k))
2561 es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2562 pprime = p(i,k+1) + (eps-1.) * es
2563 qs = eps * es / pprime
2564 dqsdp = - qs / pprime
2565 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
2566 dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
2567 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
2568 dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
2569 dq = dqsdt * dt + dqsdp * dp
2570 to(i,k) = to(i,k+1) + dt
2571 qo(i,k) = qo(i,k+1) + dq
2572 po(i,k) = .5 * (p(i,k) + p(i,k+1))
2579 if (cnvflg(i) .and. k .le. kmax(i)-1) then
2580 qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2581 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + (eps-1.) * qeso(i,k))
2583 qeso(i,k) = max(qeso(i,k), val1)
2585 qo(i,k) = max(qo(i,k), val2 )
2586 heo(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
2587 cp_ * to(i,k) + hvap_ * qo(i,k)
2588 heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
2589 cp_ * to(i,k) + hvap_ * qeso(i,k)
2590 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
2591 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
2596 ! look for the level of free convection as cloud base
2600 if(flg(i)) kbcon(i) = kmax(i)
2605 if (flg(i).and.k.lt.kbm(i)) then
2606 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
2616 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
2622 totflg = totflg .and. (.not. cnvflg(i))
2626 ! determine critical convective inhibition
2627 ! as a function of vertical velocity at cloud base.
2631 pdot(i) = 10.* dot(i,kbcon(i))
2637 if(slimsk(i).eq.1.) then
2648 if(pdot(i).le.w4) then
2649 ptem = (pdot(i) - w4) / (w3 - w4)
2650 elseif(pdot(i).ge.-w4) then
2651 ptem = - (pdot(i) + w4) / (w4 - w3)
2656 ptem = max(ptem,val1)
2658 ptem = min(ptem,val2)
2660 ptem1= .5*(cincrmax-cincrmin)
2661 cincr = cincrmax - ptem * ptem1
2662 tem1 = p(i,kb(i)) - p(i,kbcon(i))
2663 if(tem1.gt.cincr) then
2671 totflg = totflg .and. (.not. cnvflg(i))
2675 ! assume the detrainment rate for the updrafts to be same as
2676 ! the entrainment rate at cloud base
2680 xlamud(i) = xlamue(i,kbcon(i))
2684 ! determine updraft mass flux for the subcloud layers
2689 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
2690 dz = zi(i,k+2) - zi(i,k+1)
2691 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
2692 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
2698 ! compute mass flux above cloud base
2703 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
2704 dz = zi(i,k+1) - zi(i,k)
2705 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
2706 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
2712 ! compute updraft cloud property
2717 hcko(i,indx) = heo(i,indx)
2718 ucko(i,indx) = uo(i,indx)
2719 vcko(i,indx) = vo(i,indx)
2726 if(k.gt.kb(i).and.k.lt.kmax(i)) then
2727 dz = zi(i,k+1) - zi(i,k)
2728 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2729 tem1 = 0.5 * xlamud(i) * dz
2730 factor = 1. + tem - tem1
2731 ptem = 0.5 * tem + pgcon
2732 ptem1= 0.5 * tem - pgcon
2733 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
2734 (heo(i,k)+heo(i,k-1)))/factor
2735 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
2736 +ptem1*uo(i,k-1))/factor
2737 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
2738 +ptem1*vo(i,k-1))/factor
2739 dbyo(i,k) = hcko(i,k) - heso(i,k)
2745 ! taking account into convection inhibition due to existence of
2746 ! dry layers below cloud base
2755 if (flg(i).and.k.lt.kbm(i)) then
2756 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
2766 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
2772 tem = p(i,kbcon(i)) - p(i,kbcon1(i))
2773 if(tem.gt.dthk) then
2781 totflg = totflg .and. (.not. cnvflg(i))
2785 ! determine first guess cloud top as the level of zero buoyancy
2786 ! limited to the level of sigma=0.7
2790 if(flg(i)) ktcon(i) = kbm(i)
2795 if (flg(i).and.k .lt. kbm(i)) then
2796 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
2804 ! specify upper limit of mass flux at cloud base
2809 dp = 1000. * del(i,k)
2810 xmbmax(i) = dp / (g_ * dt2)
2814 ! compute cloud moisture property and precipitation
2819 qcko(i,kb(i)) = qo(i,kb(i))
2826 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
2827 dz = zi(i,k+1) - zi(i,k)
2828 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2829 qrch = qeso(i,k) + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2830 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2831 tem1 = 0.5 * xlamud(i) * dz
2832 factor = 1. + tem - tem1
2833 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
2834 (qo(i,k)+qo(i,k-1)))/factor
2835 dq = eta(i,k) * (qcko(i,k) - qrch)
2837 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
2839 ! below lfc check if there is excess moisture to release latent heat
2841 if(k.ge.kbcon(i).and.dq.gt.0.) then
2842 etah = .5 * (eta(i,k) + eta(i,k-1))
2843 if(ncloud.gt.0) then
2844 dp = 1000. * del(i,k)
2845 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
2846 dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
2848 qlk = dq / (eta(i,k) + etah * c0 * dz)
2850 aa1(i) = aa1(i) - dz * g_ * qlk
2851 qcko(i,k)= qlk + qrch
2852 pwo(i,k) = etah * c0 * dz * qlk
2859 ! calculate cloud work function
2864 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
2865 dz1 = zl(i,k+1) - zl(i,k)
2866 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2867 rfact = 1. + fv_ * cp_ * gamma * to(i,k) / hvap_
2868 aa1(i) = aa1(i) + dz1 * (g_ / (cp_ * to(i,k))) &
2869 * dbyo(i,k) / (1. + gamma) * rfact
2871 aa1(i)=aa1(i)+ dz1 * g_ * fv_ * max(val,(qeso(i,k) - qo(i,k)))
2878 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
2883 totflg = totflg .and. (.not. cnvflg(i))
2887 ! estimate the convective overshooting as the level
2888 ! where the [aafac * cloud work function] becomes zero,
2889 ! which is the final cloud top limited to the level of sigma=0.7
2893 aa1(i) = aafac * aa1(i)
2905 if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
2906 dz1 = zl(i,k+1) - zl(i,k)
2907 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2908 rfact = 1. + fv_ * cp_ * gamma &
2911 dz1 * (g_ / (cp_ * to(i,k))) &
2912 * dbyo(i,k) / (1. + gamma) * rfact
2913 if(aa1(i).lt.0.) then
2922 ! compute cloud moisture property, detraining cloud water
2923 ! and precipitation in overshooting layers
2928 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
2929 dz = zi(i,k+1) - zi(i,k)
2930 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2932 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2933 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2934 tem1 = 0.5 * xlamud(i) * dz
2935 factor = 1. + tem - tem1
2936 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
2937 (qo(i,k)+qo(i,k-1)))/factor
2938 dq = eta(i,k) * (qcko(i,k) - qrch)
2940 ! check if there is excess moisture to release latent heat
2943 etah = .5 * (eta(i,k) + eta(i,k-1))
2944 if(ncloud.gt.0) then
2945 dp = 1000. * del(i,k)
2946 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
2947 dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
2949 qlk = dq / (eta(i,k) + etah * c0 * dz)
2951 qcko(i,k) = qlk + qrch
2952 pwo(i,k) = etah * c0 * dz * qlk
2959 ! exchange ktcon with ktcon1
2964 ktcon(i) = ktcon1(i)
2969 ! this section is ready for cloud water
2971 if(ncloud.gt.0) then
2973 ! compute liquid and vapor separation at cloud top
2978 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2980 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2981 dq = qcko(i,k) - qrch
2983 ! check if there is excess moisture to release latent heat
2994 !--- compute precipitation efficiency in terms of windshear
3005 if(k.gt.kb(i).and.k.le.ktcon(i)) then
3006 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 + (vo(i,k)-vo(i,k-1)) ** 2)
3007 vshear(i) = vshear(i) + shear
3015 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i)+1)-zi(i,kb(i)+1))
3016 e1=1.591-.639*vshear(i) &
3017 +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
3020 edt(i) = min(edt(i),val)
3022 edt(i) = max(edt(i),val)
3026 !--- what would the change be, that a cloud with unit mass
3027 !--- will do to the environment?
3031 if(cnvflg(i) .and. k .le. kmax(i)) then
3040 !--- changed due to subsidence and entrainment
3045 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
3046 dp = 1000. * del(i,k)
3047 dz = zi(i,k+1) - zi(i,k)
3050 dv2h = .5 * (heo(i,k) + heo(i,k-1))
3053 dv2q = .5 * (qo(i,k) + qo(i,k-1))
3056 dv2u = .5 * (uo(i,k) + uo(i,k-1))
3059 dv2v = .5 * (vo(i,k) + vo(i,k-1))
3062 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
3065 dellah(i,k) = dellah(i,k) + &
3066 ( eta(i,k)*dv1h - eta(i,k-1)*dv3h &
3067 - tem*eta(i,k-1)*dv2h*dz &
3068 + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz ) *g_/dp
3070 dellaq(i,k) = dellaq(i,k) + &
3071 ( eta(i,k)*dv1q - eta(i,k-1)*dv3q &
3072 - tem*eta(i,k-1)*dv2q*dz &
3073 + tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz ) *g_/dp
3075 dellau(i,k) = dellau(i,k) + &
3076 ( eta(i,k)*dv1u - eta(i,k-1)*dv3u &
3077 - tem*eta(i,k-1)*dv2u*dz &
3078 + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
3079 - pgcon*eta(i,k-1)*(dv1u-dv3u) ) *g_/dp
3081 dellav(i,k) = dellav(i,k) + &
3082 ( eta(i,k)*dv1v - eta(i,k-1)*dv3v &
3083 - tem*eta(i,k-1)*dv2v*dz &
3084 + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
3085 - pgcon*eta(i,k-1)*(dv1v-dv3v) ) *g_/dp
3097 dp = 1000. * del(i,indx)
3098 dv1h = heo(i,indx-1)
3099 dellah(i,indx) = eta(i,indx-1) * &
3100 (hcko(i,indx-1) - dv1h) * g_ / dp
3102 dellaq(i,indx) = eta(i,indx-1) * &
3103 (qcko(i,indx-1) - dv1q) * g_ / dp
3105 dellau(i,indx) = eta(i,indx-1) * &
3106 (ucko(i,indx-1) - dv1u) * g_ / dp
3108 dellav(i,indx) = eta(i,indx-1) * &
3109 (vcko(i,indx-1) - dv1v) * g_ / dp
3113 dellal(i,indx) = eta(i,indx-1) * &
3114 qlko_ktcon(i) * g_ / dp
3118 ! mass flux at cloud base for shallow convection
3124 ptem = g_*sflx(i)*hpbl(i)/t1(i,1)
3126 tem = po(i,k)*100. / (rd_*t1(i,k))
3127 xmb(i) = betaw*tem*wstar(i)
3128 xmb(i) = min(xmb(i),xmbmax(i))
3134 if (cnvflg(i) .and. k .le. kmax(i)) then
3135 qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
3136 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
3138 qeso(i,k) = max(qeso(i,k), val )
3155 if(k.gt.kb(i).and.k.le.ktcon(i)) then
3156 dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
3157 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
3158 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
3160 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
3161 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
3162 dp = 1000. * del(i,k)
3163 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
3164 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
3165 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
3166 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
3167 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
3176 if(k.gt.kb(i).and.k.le.ktcon(i)) then
3177 qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls &
3179 qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
3181 qeso(i,k) = max(qeso(i,k), val )
3197 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
3198 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
3208 if (k .le. kmax(i)) then
3213 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
3214 rain(i) = rain(i) + pwo(i,k) * xmb(i) * .001 * dt2
3217 if(flg(i).and.k.lt.ktcon(i)) then
3218 evef = edt(i) * evfact
3219 if(slimsk(i).eq.1.) evef=edt(i) * evfactl
3220 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
3221 / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
3222 dp = 1000. * del(i,k)
3223 if(rain(i).gt.0..and.qcond(i).lt.0.) then
3224 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rain(i))))
3225 qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
3226 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
3228 if(rain(i).gt.0..and.qcond(i).lt.0..and.delq2(i).gt.rntot(i)) then
3229 qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
3232 if(rain(i).gt.0..and.qevap(i).gt.0.) then
3233 tem = .001 * dp / g_
3234 tem1 = qevap(i) * tem
3235 if(tem1.gt.rain(i)) then
3236 qevap(i) = rain(i) / tem
3239 rain(i) = rain(i) - tem1
3241 q1(i,k) = q1(i,k) + qevap(i)
3242 t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
3243 deltv(i) = - (hvap_/cp_)*qevap(i)/dt2
3244 delq(i) = + qevap(i)/dt2
3245 delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
3247 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
3248 delqbar(i) = delqbar(i) + delq(i)*dp/g_
3249 deltbar(i) = deltbar(i) + deltv(i)*dp/g_
3257 if(rain(i).lt.0..or..not.flg(i)) rain(i) = 0.
3266 if (ncloud.gt.0) then
3271 if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
3272 tem = dellal(i,k) * xmb(i) * dt2
3273 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
3274 if (ncloud.ge.2) then
3275 qi2(i,k) = qi2(i,k) + tem * tem1 ! ice
3276 qc2(i,k) = qc2(i,k) + tem *(1.0-tem1) ! water
3278 qc2(i,k) = qc2(i,k) + tem
3287 end subroutine nscv2d
3288 !-------------------------------------------------------------------------------
3290 END MODULE module_cu_nsas