From 867dc105197efe592cf9da086c4f9a8ed7929285 Mon Sep 17 00:00:00 2001 From: "jun.wang" Date: Fri, 1 Jun 2018 19:56:46 +0000 Subject: [PATCH] FV3: this commits #refs 48136 --- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 116 +- gfsphysics/GFS_layer/GFS_typedefs.F90 | 14 +- gfsphysics/makefile | 7 +- gfsphysics/physics/mfpblt.f | 440 ++++++ gfsphysics/physics/mfscu.f | 545 ++++++++ gfsphysics/physics/samfdeepcnv.f | 358 ++++- gfsphysics/physics/samfshalcnv.f | 281 +++- gfsphysics/physics/satmedmfvdif.f | 1390 +++++++++++++++++++ 8 files changed, 3026 insertions(+), 125 deletions(-) create mode 100755 gfsphysics/physics/mfpblt.f create mode 100755 gfsphysics/physics/mfscu.f create mode 100644 gfsphysics/physics/satmedmfvdif.f diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index b4d9fae4a..68ff9a413 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -226,7 +226,7 @@ module module_physics_driver !! - determine the index of TKE (ntk) in the convectively transported tracer array (clw) !! - allocate precipitation mixing ratio cloud droplet number concentration arrays !! - Deep Convection: -!! - determine which tracers in the tracer input array undergo convective transport (valid only for the RAS and Chikira-Sugiyama schemes) and allocate a local convective transported tracer array (clw) +!! - determine which tracers in the tracer input array undergo convective transport (valid for the RAS and Chikira-Sugiyama, and SAMF schemes) and allocate a local convective transported tracer array (clw) !! - apply an adjustment to the tracers from the dynamics !! - calculate horizontal grid-related parameters needed for some parameterizations !! - calculate the maxiumum cloud base updraft speed for the Chikira-Sugiyama scheme @@ -257,7 +257,7 @@ module module_physics_driver !! - finally, accumulate surface-related diagnostics and calculate the max/min values of T and q at 2 m height. !! . !! ## Calculate the state variable tendencies due to the PBL (vertical diffusion) scheme. -!! - Call the vertical diffusion scheme (PBL) based on the following logical flags: do_shoc, hybedmf, old_monin, mstrat +!! - Call the vertical diffusion scheme (PBL) based on the following logical flags: do_shoc, hybedmf, satmedmf, old_monin, mstrat !! - the PBL scheme is expected to return tendencies of the state variables !! - If A/O/I coupling and the surface is sea ice, overwrite some surface-related variables to their states before PBL was called !! - For diagnostics, do the following: @@ -448,9 +448,9 @@ subroutine GFS_physics_driver & ims, ime, kms, kme, its, ite, kts, kte, imp_physics, & ntwa, ntia - integer :: i, kk, ic, k, n, k1, iter, levshcm, tracers, & - tottracer, num2, num3, nshocm, nshoc, ntk, nn, nncl, & - seconds + integer :: i, kk, ic, k, n, k1, iter, levshcm, tracers, & + tottracer, nsamftrac, num2, num3, nshocm, nshoc, ntk, & + nn, nncl, seconds integer, dimension(size(Grid%xlon,1)) :: & kbot, ktop, kcnv, soiltyp, vegtype, kpbl, slopetyp, kinver, & @@ -534,7 +534,7 @@ subroutine GFS_physics_driver & !--- ALLOCATABLE ELEMENTS !--- in clw, the first two varaibles are cloud water and ice. !--- from third to ntrac are convective transportable tracers, - !--- third being the ozone, when ntrac=3 (valid only with ras) + !--- third being the ozone, when ntrac=3 (valid with ras, csaw, or samf) !--- Anning Cheng 9/21/2016 leave a hook here for diagnosed snow, !--- rain, and their numbers real(kind=kind_phys), allocatable :: & @@ -579,6 +579,10 @@ subroutine GFS_physics_driver & ntiw = Model%ntiw ncld = Model%ncld ntke = Model%ntke +! +! scal-aware TKE-based moist EDMF (satmedmfvdif) scheme is coded assuming +! ntke=ntrac. If ntrac > ntke, the code needs to be modified. (Jongil Han) +! ntlnc = Model%ntlnc ntinc = Model%ntinc ntrw = Model%ntrw @@ -1445,7 +1449,18 @@ subroutine GFS_physics_driver & ! if (lprnt) write(0,*)'aftmonshoc=',Statein%tgrs(ipr,:) ! if (lprnt) write(0,*)'aftmonshocdtdt=',dtdt(ipr,1:10) else - if (Model%hybedmf) then + if (Model%satmedmf) then + call satmedmfvdif(ix, im, levs, nvdiff, ntcw, ntke, & + dvdt, dudt, dtdt, dqdt, & + Statein%ugrs, Statein%vgrs, Statein%tgrs, Statein%qgrs, & + Radtend%htrsw, Radtend%htrlw, xmu, garea, & + Statein%prsik(1,1), rb, Sfcprop%zorl, Diag%u10m, Diag%v10m, & + Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%tsfc, hflx, evap, & + stress, wind, kpbl, Statein%prsi, del, Statein%prsl, & + Statein%prslk, Statein%phii, Statein%phil, dtp, & + Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl, & + kinver, Model%xkzm_m, Model%xkzm_h, Model%xkzm_s) + elseif (Model%hybedmf) then call moninedmf(ix, im, levs, nvdiff, ntcw, dvdt, dudt, dtdt, dqdt, & Statein%ugrs, Statein%vgrs, Statein%tgrs, Statein%qgrs, & Radtend%htrsw, Radtend%htrlw, xmu, Statein%prsik(1,1), & @@ -1993,12 +2008,12 @@ subroutine GFS_physics_driver & endif -! --- ... for convective tracer transport (while using ras or csaw) +! --- ... for convective tracer transport (while using ras, csaw, or samf) ! (the code here implicitly assumes that ntiw=ntcw+1) ntk = 0 tottracer = 0 - if (Model%cscnv .or. Model%trans_trac ) then + if (Model%cscnv .or. Model%satmedmf .or. Model%trans_trac ) then otspt(:,:) = .true. ! otspt is used only for cscnv otspt(1:3,:) = .false. ! this is for sp.hum, ice and liquid water tracers = 2 @@ -2023,7 +2038,7 @@ subroutine GFS_physics_driver & endif enddo tottracer = tracers - 2 - endif ! end if_ras or cfscnv + endif ! end if_ras or cfscnv or samf ! if (kdt == 1 .and. me == 0) & ! write(0,*)' trans_trac=',Model%trans_trac,' tottracer=', & @@ -2268,16 +2283,20 @@ subroutine GFS_physics_driver & Model%evfact_deep, Model%evfactl_deep, & Model%pgcon_deep) elseif (Model%imfdeepcnv == 2) then - call samfdeepcnv(im, ix, levs, dtp, del, Statein%prsl, & - Statein%pgr, Statein%phil, clw(:,:,1:2), & - Stateout%gq0(:,:,1), & - Stateout%gt0, Stateout%gu0, Stateout%gv0, & - cld1d, rain1, kbot, ktop, kcnv, islmsk, & - garea, Statein%vvl, ncld, ud_mf, dd_mf, & - dt_mf, cnvw, cnvc, & - Model%clam_deep, Model%c0s_deep, & - Model%c1_deep, Model%betal_deep, Model%betas_deep,& - Model%evfact_deep, Model%evfactl_deep, & + if(.not. Model%satmedmf .and. .not. Model%trans_trac) then + nsamftrac = 0 + else + nsamftrac = tottracer + endif + call samfdeepcnv(im, ix, levs, dtp, ntk, nsamftrac, del, & + Statein%prsl, Statein%pgr, Statein%phil, clw, & + Stateout%gq0(:,:,1), Stateout%gt0, & + Stateout%gu0, Stateout%gv0, & + cld1d, rain1, kbot, ktop, kcnv, islmsk, garea, & + Statein%vvl, ncld, ud_mf, dd_mf, dt_mf, cnvw, cnvc, & + Model%clam_deep, Model%c0s_deep, & + Model%c1_deep, Model%betal_deep, Model%betas_deep, & + Model%evfact_deep, Model%evfactl_deep, & Model%pgcon_deep, Model%asolfac_deep) ! if (lprnt) print *,' rain1=',rain1(ipr) elseif (Model%imfdeepcnv == 0) then ! random cloud top @@ -2431,26 +2450,6 @@ subroutine GFS_physics_driver & enddo endif ! if (lgocart) -! --- ... update the tracers due to convective transport -! (except for suspended water and ice) - - if (tottracer > 0) then - tracers = 2 - do n=2,ntrac -! if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt) then - if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & - n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & - n /= ntsnc .and. n /= ntgl .and. n /= ntgnc ) then - tracers = tracers + 1 - do k=1,levs - do i=1,im - Stateout%gq0(i,k,n) = clw(i,k,tracers) - enddo - enddo - endif - enddo - endif - endif ! end if_not_ras else ! no parameterized deep convection cld1d = 0. @@ -2779,10 +2778,15 @@ subroutine GFS_physics_driver & endif elseif (Model%imfshalcnv == 2) then - call samfshalcnv (im, ix, levs, dtp, del, Statein%prsl, & - Statein%pgr, Statein%phil, clw(:,:,1:2), & - Stateout%gq0(:,:,1:1), & - Stateout%gt0, Stateout%gu0, Stateout%gv0, & + if(.not. Model%satmedmf .and. .not. Model%trans_trac) then + nsamftrac = 0 + else + nsamftrac = tottracer + endif + call samfshalcnv (im, ix, levs, dtp, ntk, nsamftrac, del, & + Statein%prsl, Statein%pgr, Statein%phil, clw, & + Stateout%gq0(:,:,1), Stateout%gt0, & + Stateout%gu0, Stateout%gv0, & rain1, kbot, ktop, kcnv, islmsk, garea, & Statein%vvl, ncld, DIag%hpbl, ud_mf, & dt_mf, cnvw, cnvc, & @@ -2977,7 +2981,29 @@ subroutine GFS_physics_driver & ! write(0,*) ' aftshgt0=',gt0(ipr,:) ! write(0,*) ' aftshgq0=',gq0(ipr,:,1) ! endif - +! +!------------------------------------------------------------------------------ +! --- update the tracers due to deep & shallow cumulus convective transport +! (except for suspended water and ice) +! + if (tottracer > 0) then + tracers = 2 + do n=2,ntrac +! if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt) then + if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & + n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & + n /= ntsnc .and. n /= ntgl .and. n /= ntgnc ) then + tracers = tracers + 1 + do k=1,levs + do i=1,im + Stateout%gq0(i,k,n) = clw(i,k,tracers) + enddo + enddo + endif + enddo + endif +!------------------------------------------------------------------------------- +! if (ntcw > 0) then ! for microphysics diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 5dc3f32a0..076e1d37d 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -483,7 +483,7 @@ module GFS_typedefs logical :: ras !< flag for ras convection scheme logical :: flipv !< flag for vertical direction flip (ras) !< .true. implies surface at k=1 - logical :: trans_trac !< flag for convective transport of tracers (RAS only) + logical :: trans_trac !< flag for convective transport of tracers (RAS, CS, or SAMF) logical :: old_monin !< flag for diff monin schemes logical :: cnvgwd !< flag for conv gravity wave drag logical :: mstrat !< flag for moorthi approach for stratus @@ -502,6 +502,8 @@ module GFS_typedefs logical :: shcnvcw !< flag for shallow convective cloud logical :: redrag !< flag for reduced drag coeff. over sea logical :: hybedmf !< flag for hybrid edmf pbl scheme + logical :: satmedmf !< flag for scale-aware TKE-based moist edmf + !< vertical turbulent mixing scheme logical :: dspheat !< flag for tke dissipative heating logical :: cnvcld logical :: random_clds !< flag controls whether clouds are random @@ -1545,7 +1547,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: ras = .false. !< flag for ras convection scheme logical :: flipv = .true. !< flag for vertical direction flip (ras) !< .true. implies surface at k=1 - logical :: trans_trac = .false. !< flag for convective transport of tracers (RAS only) + logical :: trans_trac = .false. !< flag for convective transport of tracers (RAS, CS, or SAMF) logical :: old_monin = .false. !< flag for diff monin schemes logical :: cnvgwd = .false. !< flag for conv gravity wave drag logical :: mstrat = .false. !< flag for moorthi approach for stratus @@ -1563,6 +1565,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: shcnvcw = .false. !< flag for shallow convective cloud logical :: redrag = .false. !< flag for reduced drag coeff. over sea logical :: hybedmf = .false. !< flag for hybrid edmf pbl scheme + logical :: satmedmf = .false. !< flag for scale-aware TKE-based moist edmf + !< vertical turbulent mixing scheme logical :: dspheat = .false. !< flag for tke dissipative heating logical :: cnvcld = .false. logical :: random_clds = .false. !< flag controls whether clouds are random @@ -1702,7 +1706,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- physical parameterizations ras, trans_trac, old_monin, cnvgwd, mstrat, moist_adj, & cscnv, cal_pre, do_aw, do_shoc, shocaftcnv, shoc_cld, & - h2o_phys, pdfcld, shcnvcw, redrag, hybedmf, dspheat, cnvcld, & + h2o_phys, pdfcld, shcnvcw, redrag, hybedmf, satmedmf, & + dspheat, cnvcld, & random_clds, shal_cnv, imfshalcnv, imfdeepcnv, do_deep, jcap,& cs_parm, flgmin, cgwf, ccwf, cdmbgwd, sup, ctei_rm, crtrh, & dlqf, rbcr, shoc_parm, & @@ -1899,6 +1904,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%shcnvcw = shcnvcw Model%redrag = redrag Model%hybedmf = hybedmf + Model%satmedmf = satmedmf Model%dspheat = dspheat Model%cnvcld = cnvcld Model%random_clds = random_clds @@ -2060,6 +2066,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%shal_cnv = .false. Model%imfshalcnv = -1 Model%hybedmf = .false. + Model%satmedmf = .false. if (Model%me == Model%master) print *,' Simplified Higher Order Closure Model used for', & ' Boundary layer and Shallow Convection', & ' nshoc_3d=',Model%nshoc_3d, & @@ -2455,6 +2462,7 @@ subroutine control_print(Model) print *, ' shcnvcw : ', Model%shcnvcw print *, ' redrag : ', Model%redrag print *, ' hybedmf : ', Model%hybedmf + print *, ' satmedmf : ', Model%satmedmf print *, ' dspheat : ', Model%dspheat print *, ' cnvcld : ', Model%cnvcld print *, ' random_clds : ', Model%random_clds diff --git a/gfsphysics/makefile b/gfsphysics/makefile index c31a45e73..4708b33ef 100644 --- a/gfsphysics/makefile +++ b/gfsphysics/makefile @@ -52,9 +52,9 @@ SRCS_f = \ ./physics/iounitdef.f \ ./physics/lrgsclr.f \ ./physics/mersenne_twister.f \ - ./physics/samfdeepcnv.f \ ./physics/mfpbl.f \ - ./physics/samfshalcnv.f \ + ./physics/mfpblt.f \ + ./physics/mfscu.f \ ./physics/module_bfmicrophysics.f \ ./physics/moninedmf.f \ ./physics/moninp.f \ @@ -92,8 +92,11 @@ SRCS_f = \ ./physics/rascnvv2.f \ ./physics/rayleigh_damp.f \ ./physics/rayleigh_damp_mesopause.f \ + ./physics/samfdeepcnv.f \ + ./physics/samfshalcnv.f \ ./physics/sascnv.f \ ./physics/sascnvn.f \ + ./physics/satmedmfvdif.f \ ./physics/set_soilveg.f \ ./physics/sfc_cice.f \ ./physics/sfc_diag.f \ diff --git a/gfsphysics/physics/mfpblt.f b/gfsphysics/physics/mfpblt.f new file mode 100755 index 000000000..3a09ad13a --- /dev/null +++ b/gfsphysics/physics/mfpblt.f @@ -0,0 +1,440 @@ + subroutine mfpblt(im,ix,km,kmpbl,ntcw,ntrac1,delt, + & cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx, + & gdx,hpbl,kpbl,vpert,buo,xmf, + & tcko,qcko,ucko,vcko,xlamue) +! + use machine , only : kind_phys + use funcphys , only : fpvs + use physcons, grav => con_g, cp => con_cp + &, rv => con_rv, hvap => con_hvap + &, fv => con_fvirt + &, eps => con_eps, epsm1 => con_epsm1 +! + implicit none +! + integer im, ix, km, kmpbl, ntcw, ntrac1 +! &, me + integer kpbl(im) + logical cnvflg(im) + real(kind=kind_phys) delt + real(kind=kind_phys) q1(ix,km,ntrac1), + & t1(ix,km), u1(ix,km), v1(ix,km), + & plyr(im,km),pix(im,km),thlx(im,km), + & thvx(im,km),zl(im,km), zm(im,km), + & gdx(im), + & hpbl(im), vpert(im), + & buo(im,km), xmf(im,km), + & tcko(im,km),qcko(im,km,ntrac1), + & ucko(im,km),vcko(im,km), + & xlamue(im,km-1) +! +c local variables and arrays +! + integer i, j, k, n, ndc + integer kpblx(im), kpbly(im) +! + real(kind=kind_phys) dt2, dz, ce0, cm, + & factor, gocp, + & g, b1, f1, + & bb1, bb2, + & alp, a1, pgcon, + & qmin, qlmin, xmmx, rbint, + & tem, tem1, tem2, + & ptem, ptem1, ptem2 +! + real(kind=kind_phys) elocp, el2orc, qs, es, + & tlu, gamma, qlu, + & thup, thvu, dq +! + real(kind=kind_phys) rbdn(im), rbup(im), hpblx(im), + & xlamuem(im,km-1) +! + real(kind=kind_phys) wu2(im,km), thlu(im,km), + & qtx(im,km), qtu(im,km) +! + real(kind=kind_phys) xlamavg(im), sigma(im), + & scaldfunc(im), sumx(im) +! + logical totflg, flg(im) +! +! physical parameters + parameter(g=grav) + parameter(gocp=g/cp) + parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) + parameter(ce0=0.4,cm=1.0) + parameter(qmin=1.e-8,qlmin=1.e-12) + parameter(alp=1.0,pgcon=0.55) + parameter(a1=0.13,b1=0.5,f1=0.15) +! +!************************************************************************ +!! + totflg = .true. + do i=1,im + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +! + dt2 = delt +!! + do k = 1, km + do i=1,im + if (cnvflg(i)) then + buo(i,k) = 0. + wu2(i,k) = 0. + qtx(i,k) = q1(i,k,1) + q1(i,k,ntcw) + endif + enddo + enddo +! +! compute thermal excess +! + do i=1,im + if(cnvflg(i)) then + ptem = alp * vpert(i) + ptem = min(ptem, 3.0) + thlu(i,1)= thlx(i,1) + ptem + qtu(i,1) = qtx(i,1) + buo(i,1) = g * ptem / thvx(i,1) + endif + enddo +! +! compute entrainment rate +! + do k = 1, kmpbl + do i=1,im + if(cnvflg(i)) then + dz = zl(i,k+1) - zl(i,k) + if(k < kpbl(i)) then + ptem = 1./(zm(i,k)+dz) + tem = max((hpbl(i)-zm(i,k)+dz) ,dz) + ptem1 = 1./tem + xlamue(i,k) = ce0 * (ptem+ptem1) + else + xlamue(i,k) = ce0 / dz + endif + xlamuem(i,k) = cm * xlamue(i,k) + endif + enddo + enddo +! +! compute buoyancy for updraft air parcel +! + do k = 2, kmpbl + do i=1,im + if(cnvflg(i)) then + dz = zl(i,k) - zl(i,k-1) + tem = 0.5 * xlamue(i,k-1) * dz + factor = 1. + tem +! + thlu(i,k) = ((1.-tem)*thlu(i,k-1)+tem* + & (thlx(i,k-1)+thlx(i,k)))/factor + qtu(i,k) = ((1.-tem)*qtu(i,k-1)+tem* + & (qtx(i,k-1)+qtx(i,k)))/factor +! + tlu = thlu(i,k) / pix(i,k) + es = 0.01 * fpvs(tlu) ! fpvs in pa + qs = max(qmin, eps * es / (plyr(i,k)+epsm1*es)) + dq = qtu(i,k) - qs +! + if (dq > 0.) then + gamma = el2orc * qs / (tlu**2) + qlu = dq / (1. + gamma) + qtu(i,k) = qs + qlu + tem1 = 1. + fv * qs - qlu + thup = thlu(i,k) + pix(i,k) * elocp * qlu + thvu = thup * tem1 + else + tem1 = 1. + fv * qtu(i,k) + thvu = thlu(i,k) * tem1 + endif + buo(i,k) = g * (thvu / thvx(i,k) - 1.) +! + endif + enddo + enddo +! +! compute updraft velocity square(wu2) +! +! tem = 1.-2.*f1 +! bb1 = 2. * b1 / tem +! bb2 = 2. / tem +! from Soares et al. (2004,QJRMS) +! bb1 = 2. +! bb2 = 4. +! +! from Bretherton et al. (2004, MWR) +! bb1 = 4. +! bb2 = 2. +! +! from our tuning + bb1 = 2.0 + bb2 = 4.0 +! + do i = 1, im + if(cnvflg(i)) then + dz = zm(i,1) + tem = 0.5*bb1*xlamue(i,1)*dz + tem1 = bb2 * buo(i,1) * dz + ptem1 = 1. + tem + wu2(i,1) = tem1 / ptem1 + endif + enddo + do k = 2, kmpbl + do i = 1, im + if(cnvflg(i)) then + dz = zm(i,k) - zm(i,k-1) + tem = 0.25*bb1*(xlamue(i,k)+xlamue(i,k-1))*dz + tem1 = bb2 * buo(i,k) * dz + ptem = (1. - tem) * wu2(i,k-1) + ptem1 = 1. + tem + wu2(i,k) = (ptem + tem1) / ptem1 + endif + enddo + enddo +! +! update pbl height as the height where updraft velocity vanishes +! + do i=1,im + flg(i) = .true. + kpbly(i) = kpbl(i) + if(cnvflg(i)) then + flg(i) = .false. + rbup(i) = wu2(i,1) + endif + enddo + do k = 2, kmpbl + do i = 1, im + if(.not.flg(i)) then + rbdn(i) = rbup(i) + rbup(i) = wu2(i,k) + kpblx(i)= k + flg(i) = rbup(i).le.0. + endif + enddo + enddo + do i = 1,im + if(cnvflg(i)) then + k = kpblx(i) + if(rbdn(i) <= 0.) then + rbint = 0. + elseif(rbup(i) >= 0.) then + rbint = 1. + else + rbint = rbdn(i)/(rbdn(i)-rbup(i)) + endif + hpblx(i) = zm(i,k-1) + rbint*(zm(i,k)-zm(i,k-1)) + endif + enddo +! + do i = 1,im + if(cnvflg(i)) then + if(kpbl(i) > kpblx(i)) then + kpbl(i) = kpblx(i) + hpbl(i) = hpblx(i) + endif + endif + enddo +! +! update entrainment rate +! + do k = 1, kmpbl + do i=1,im + if(cnvflg(i) .and. kpbly(i) > kpblx(i)) then + dz = zl(i,k+1) - zl(i,k) + if(k < kpbl(i)) then + ptem = 1./(zm(i,k)+dz) + tem = max((hpbl(i)-zm(i,k)+dz) ,dz) + ptem1 = 1./tem + xlamue(i,k) = ce0 * (ptem+ptem1) + else + xlamue(i,k) = ce0 / dz + endif + xlamuem(i,k) = cm * xlamue(i,k) + endif + enddo + enddo +! +! compute entrainment rate averaged over the whole pbl +! + do i = 1, im + xlamavg(i) = 0. + sumx(i) = 0. + enddo + do k = 1, kmpbl + do i = 1, im + if (cnvflg(i) .and. k < kpbl(i)) then + dz = zl(i,k+1) - zl(i,k) + xlamavg(i) = xlamavg(i) + xlamue(i,k) * dz + sumx(i) = sumx(i) + dz + endif + enddo + enddo + do i = 1, im + if(cnvflg(i)) then + xlamavg(i) = xlamavg(i) / sumx(i) + endif + enddo +! +! updraft mass flux as a function of updraft velocity profile +! + do k = 1, kmpbl + do i = 1, im + if (cnvflg(i) .and. k < kpbl(i)) then + if(wu2(i,k) > 0.) then + tem = sqrt(wu2(i,k)) + else + tem = 0. + endif + xmf(i,k) = a1 * tem + endif + enddo + enddo +! +!--- compute updraft fraction as a function of mean entrainment rate +! (Grell & Freitas, 2014) +! + do i = 1, im + if(cnvflg(i)) then + tem = 0.2 / xlamavg(i) + tem1 = 3.14 * tem * tem + sigma(i) = tem1 / (gdx(i) * gdx(i)) + sigma(i) = max(sigma(i), 0.001) + sigma(i) = min(sigma(i), 0.999) + endif + enddo +! +!--- compute scale-aware function based on Arakawa & Wu (2013) +! + do i = 1, im + if(cnvflg(i)) then + if (sigma(i) > a1) then + scaldfunc(i) = (1.-sigma(i)) * (1.-sigma(i)) + scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) + else + scaldfunc(i) = 1.0 + endif + endif + enddo +! +! final scale-aware updraft mass flux +! + do k = 1, kmpbl + do i = 1, im + if (cnvflg(i) .and. k < kpbl(i)) then + xmf(i,k) = scaldfunc(i) * xmf(i,k) + dz = zl(i,k+1) - zl(i,k) + xmmx = dz / dt2 + xmf(i,k) = min(xmf(i,k),xmmx) + endif + enddo + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute updraft property using updated entranment rate +! + do i=1,im + if(cnvflg(i)) then + thlu(i,1)= thlx(i,1) + endif + enddo +! +! do i=1,im +! if(cnvflg(i)) then +! ptem1 = max(qcko(i,1,ntcw), 0.) +! tlu = thlu(i,1) / pix(i,1) +! tcko(i,1) = tlu + elocp * ptem1 +! endif +! enddo +! + do k = 2, kmpbl + do i=1,im + if(cnvflg(i) .and. k <= kpbl(i)) then + dz = zl(i,k) - zl(i,k-1) + tem = 0.5 * xlamue(i,k-1) * dz + factor = 1. + tem +! + thlu(i,k) = ((1.-tem)*thlu(i,k-1)+tem* + & (thlx(i,k-1)+thlx(i,k)))/factor + qtu(i,k) = ((1.-tem)*qtu(i,k-1)+tem* + & (qtx(i,k-1)+qtx(i,k)))/factor +! + tlu = thlu(i,k) / pix(i,k) + es = 0.01 * fpvs(tlu) ! fpvs in pa + qs = max(qmin, eps * es / (plyr(i,k)+epsm1*es)) + dq = qtu(i,k) - qs +! + if (dq > 0.) then + gamma = el2orc * qs / (tlu**2) + qlu = dq / (1. + gamma) + qtu(i,k) = qs + qlu + qcko(i,k,1) = qs + qcko(i,k,ntcw) = qlu + tcko(i,k) = tlu + elocp * qlu + else + qcko(i,k,1) = qtu(i,k) + qcko(i,k,ntcw) = 0. + tcko(i,k) = tlu + endif +! + endif + enddo + enddo +! + do k = 2, kmpbl + do i = 1, im + if (cnvflg(i) .and. k <= kpbl(i)) then + dz = zl(i,k) - zl(i,k-1) + tem = 0.5 * xlamuem(i,k-1) * dz + factor = 1. + tem + ptem = tem + pgcon + ptem1= tem - pgcon + ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*u1(i,k) + & +ptem1*u1(i,k-1))/factor + vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*v1(i,k) + & +ptem1*v1(i,k-1))/factor + endif + enddo + enddo +! + if(ntcw > 2) then +! + do n = 2, ntcw-1 + do k = 2, kmpbl + do i = 1, im + if (cnvflg(i) .and. k <= kpbl(i)) then + dz = zl(i,k) - zl(i,k-1) + tem = 0.5 * xlamue(i,k-1) * dz + factor = 1. + tem +! + qcko(i,k,n) = ((1.-tem)*qcko(i,k-1,n)+tem* + & (q1(i,k,n)+q1(i,k-1,n)))/factor + endif + enddo + enddo + enddo +! + endif +! + ndc = ntrac1 - ntcw +! + if(ndc > 0) then +! + do n = ntcw+1, ntrac1 + do k = 2, kmpbl + do i = 1, im + if (cnvflg(i) .and. k <= kpbl(i)) then + dz = zl(i,k) - zl(i,k-1) + tem = 0.5 * xlamue(i,k-1) * dz + factor = 1. + tem +! + qcko(i,k,n) = ((1.-tem)*qcko(i,k-1,n)+tem* + & (q1(i,k,n)+q1(i,k-1,n)))/factor + endif + enddo + enddo + enddo +! + endif +! + return + end diff --git a/gfsphysics/physics/mfscu.f b/gfsphysics/physics/mfscu.f new file mode 100755 index 000000000..692950bd2 --- /dev/null +++ b/gfsphysics/physics/mfscu.f @@ -0,0 +1,545 @@ + subroutine mfscu(im,ix,km,kmscu,ntcw,ntrac1,delt, + & cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix, + & thlx,thvx,thlvx,gdx,thetae,radj, + & krad,mrad,radmin,buo,xmfd, + & tcdo,qcdo,ucdo,vcdo,xlamde) +! + use machine , only : kind_phys + use funcphys , only : fpvs + use physcons, grav => con_g, cp => con_cp + &, rv => con_rv, hvap => con_hvap + &, fv => con_fvirt + &, eps => con_eps, epsm1 => con_epsm1 +! + implicit none +! + integer im, ix, km, kmscu, ntcw, ntrac1 +! &, me + integer krad(im), mrad(im) +! + logical cnvflg(im) + real(kind=kind_phys) delt + real(kind=kind_phys) q1(ix,km,ntrac1),t1(ix,km), + & u1(ix,km), v1(ix,km), + & plyr(im,km), pix(im,km), + & thlx(im,km), + & thvx(im,km), thlvx(im,km), + & gdx(im), radj(im), + & zl(im,km), zm(im,km), + & thetae(im,km), radmin(im), + & buo(im,km), xmfd(im,km), + & tcdo(im,km), qcdo(im,km,ntrac1), + & ucdo(im,km), vcdo(im,km), + & xlamde(im,km-1) +! +! local variables and arrays +! +! + integer i,j,indx, k, n, kk, ndc + integer krad1(im), mradx(im), mrady(im) +! + real(kind=kind_phys) dt2, dz, ce0, cm, + & gocp, factor, g, tau, + & b1, f1, bb1, bb2, + & a1, a2, a11, a22, + & cteit, pgcon, + & qmin, qlmin, + & xmmx, tem, tem1, tem2, + & ptem, ptem1, ptem2 +! + real(kind=kind_phys) elocp, el2orc, qs, es, + & tld, gamma, qld, thdn, + & thvd, dq +! + real(kind=kind_phys) wd2(im,km), thld(im,km), + & qtx(im,km), qtd(im,km), + & thlvd(im), hrad(im), + & xlamdem(im,km-1), ra1(im), ra2(im) +! + real(kind=kind_phys) xlamavg(im), sigma(im), + & scaldfunc(im), sumx(im) +! + logical totflg, flg(im) +! + real(kind=kind_phys) actei, cldtime +! +c physical parameters + parameter(g=grav) + parameter(gocp=g/cp) + parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) + parameter(ce0=0.4,cm=1.0,pgcon=0.55) + parameter(qmin=1.e-8,qlmin=1.e-12) + parameter(b1=0.45,f1=0.15) + parameter(a1=0.12,a2=0.5) + parameter(a11=0.2,a22=1.0) + parameter(cldtime=500.) + parameter(actei = 0.7) +! parameter(actei = 0.23) +! +!************************************************************************ +!! + totflg = .true. + do i=1,im + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +! + dt2 = delt +!! + do k = 1, km + do i=1,im + if(cnvflg(i)) then + buo(i,k) = 0. + wd2(i,k) = 0. + qtx(i,k) = q1(i,k,1) + q1(i,k,ntcw) + endif + enddo + enddo +! + do i = 1, im + if(cnvflg(i)) then + hrad(i) = zm(i,krad(i)) + krad1(i) = krad(i)-1 + endif + enddo +! + do i = 1, im + if(cnvflg(i)) then + k = krad(i) + tem = zm(i,k+1)-zm(i,k) + tem1 = cldtime*radmin(i)/tem + tem1 = max(tem1, -3.0) + thld(i,k)= thlx(i,k) + tem1 + qtd(i,k) = qtx(i,k) + thlvd(i) = thlvx(i,k) + tem1 + buo(i,k) = - g * tem1 / thvx(i,k) + endif + enddo +! +! specify downdraft fraction +! + do i=1,im + if(cnvflg(i)) then + ra1(i) = a1 + ra2(i) = a11 + endif + enddo +! +! if the condition for cloud-top instability is met, +! increase downdraft fraction +! + do i = 1, im + if(cnvflg(i)) then + k = krad(i) + tem = thetae(i,k) - thetae(i,k+1) + tem1 = qtx(i,k) - qtx(i,k+1) + if (tem > 0. .and. tem1 > 0.) then + cteit= cp*tem/(hvap*tem1) + if(cteit > actei) then + ra1(i) = a2 + ra2(i) = a22 + endif + endif + endif + enddo +! +! compute radiative flux jump at stratocumulus top +! + do i = 1, im + if(cnvflg(i)) then + radj(i) = -ra2(i) * radmin(i) + endif + enddo +! +! first-quess level of downdraft extension (mrad) +! + do i = 1, im + flg(i) = cnvflg(i) + mrad(i) = krad(i) + enddo + do k = kmscu,1,-1 + do i = 1, im + if(flg(i) .and. k < krad(i)) then + if(thlvd(i) <= thlvx(i,k)) then + mrad(i) = k + else + flg(i)=.false. + endif + endif + enddo + enddo + do i=1,im + if (cnvflg(i)) then + kk = krad(i)-mrad(i) + if(kk < 1) cnvflg(i)=.false. + endif + enddo +!! + totflg = .true. + do i=1,im + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +!! +! +! compute entrainment rate +! + do k = 1, kmscu + do i=1,im + if(cnvflg(i)) then + dz = zl(i,k+1) - zl(i,k) + if(k >= mrad(i) .and. k < krad(i)) then + if(mrad(i) == 1) then + ptem = 1./(zm(i,k)+dz) + else + ptem = 1./(zm(i,k)-zm(i,mrad(i)-1)+dz) + endif + tem = max((hrad(i)-zm(i,k)+dz) ,dz) + ptem1 = 1./tem + xlamde(i,k) = ce0 * (ptem+ptem1) + else + xlamde(i,k) = ce0 / dz + endif + xlamdem(i,k) = cm * xlamde(i,k) + endif + enddo + enddo +! +! compute buoyancy for downdraft air parcel +! + do k = kmscu,1,-1 + do i=1,im + if(cnvflg(i) .and. k < krad(i)) then + dz = zl(i,k+1) - zl(i,k) + tem = 0.5 * xlamde(i,k) * dz + factor = 1. + tem +! + thld(i,k) = ((1.-tem)*thld(i,k+1)+tem* + & (thlx(i,k)+thlx(i,k+1)))/factor + qtd(i,k) = ((1.-tem)*qtd(i,k+1)+tem* + & (qtx(i,k)+qtx(i,k+1)))/factor +! + tld = thld(i,k) / pix(i,k) + es = 0.01 * fpvs(tld) ! fpvs in pa + qs = max(qmin, eps * es / (plyr(i,k)+epsm1*es)) + dq = qtd(i,k) - qs +! + if (dq > 0.) then + gamma = el2orc * qs / (tld**2) + qld = dq / (1. + gamma) + qtd(i,k) = qs + qld + tem1 = 1. + fv * qs - qld + thdn = thld(i,k) + pix(i,k) * elocp * qld + thvd = thdn * tem1 + else + tem1 = 1. + fv * qtd(i,k) + thvd = thld(i,k) * tem1 + endif + buo(i,k) = g * (1. - thvd / thvx(i,k)) +! + endif + enddo + enddo +! +! compute downdraft velocity square(wd2) +! +! tem = 1.-2.*f1 +! bb1 = 2. * b1 / tem +! bb2 = 2. / tem +! from Soares et al. (2004,QJRMS) +! bb1 = 2. +! bb2 = 4. +! +! from Bretherton et al. (2004, MWR) +! bb1 = 4. +! bb2 = 2. +! +! from our tuning + bb1 = 2.0 + bb2 = 4.0 +! + do i = 1, im + if(cnvflg(i)) then + k = krad1(i) + dz = zm(i,k+1) - zm(i,k) +! tem = 0.25*bb1*(xlamde(i,k)+xlamde(i,k+1))*dz + tem = 0.5*bb1*xlamde(i,k)*dz + tem1 = bb2 * buo(i,k+1) * dz + ptem1 = 1. + tem + wd2(i,k) = tem1 / ptem1 + endif + enddo + do k = kmscu,1,-1 + do i = 1, im + if(cnvflg(i) .and. k < krad1(i)) then + dz = zm(i,k+1) - zm(i,k) + tem = 0.25*bb1*(xlamde(i,k)+xlamde(i,k+1))*dz + tem1 = bb2 * buo(i,k+1) * dz + ptem = (1. - tem) * wd2(i,k+1) + ptem1 = 1. + tem + wd2(i,k) = (ptem + tem1) / ptem1 + endif + enddo + enddo +c + do i = 1, im + flg(i) = cnvflg(i) + mrady(i) = mrad(i) + if(flg(i)) mradx(i) = krad(i) + enddo + do k = kmscu,1,-1 + do i = 1, im + if(flg(i) .and. k < krad(i)) then + if(wd2(i,k) > 0.) then + mradx(i) = k + else + flg(i)=.false. + endif + endif + enddo + enddo +! + do i = 1,im + if(cnvflg(i)) then + if(mrad(i) < mradx(i)) then + mrad(i) = mradx(i) + endif + endif + enddo +! + do i=1,im + if (cnvflg(i)) then + kk = krad(i)-mrad(i) + if(kk < 1) cnvflg(i)=.false. + endif + enddo +!! + totflg = .true. + do i=1,im + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +!! +! +! update entrainment rate +! + do k = 1, kmscu + do i=1,im + if(cnvflg(i) .and. mrady(i) < mradx(i)) then + dz = zl(i,k+1) - zl(i,k) + if(k >= mrad(i) .and. k < krad(i)) then + if(mrad(i) == 1) then + ptem = 1./(zm(i,k)+dz) + else + ptem = 1./(zm(i,k)-zm(i,mrad(i)-1)+dz) + endif + tem = max((hrad(i)-zm(i,k)+dz) ,dz) + ptem1 = 1./tem + xlamde(i,k) = ce0 * (ptem+ptem1) + else + xlamde(i,k) = ce0 / dz + endif + xlamdem(i,k) = cm * xlamde(i,k) + endif + enddo + enddo +! +! compute entrainment rate averaged over the whole downdraft layers +! + do i = 1, im + xlamavg(i) = 0. + sumx(i) = 0. + enddo + do k = kmscu, 1, -1 + do i = 1, im + if(cnvflg(i) .and. + & (k >= mrad(i) .and. k < krad(i))) then + dz = zl(i,k+1) - zl(i,k) + xlamavg(i) = xlamavg(i) + xlamde(i,k) * dz + sumx(i) = sumx(i) + dz + endif + enddo + enddo + do i = 1, im + if(cnvflg(i)) then + xlamavg(i) = xlamavg(i) / sumx(i) + endif + enddo +! +! compute downdraft mass flux +! + do k = kmscu, 1, -1 + do i = 1, im + if(cnvflg(i) .and. + & (k >= mrad(i) .and. k < krad(i))) then + if(wd2(i,k) > 0.) then + tem = sqrt(wd2(i,k)) + else + tem = 0. + endif + xmfd(i,k) = ra1(i) * tem + endif + enddo + enddo +! +!--- compute downdraft fraction as a function of mean entrainment rate +! (Grell & Freitas, 2014) +! + do i = 1, im + if(cnvflg(i)) then + tem = 0.2 / xlamavg(i) + tem1 = 3.14 * tem * tem + sigma(i) = tem1 / (gdx(i) * gdx(i)) + sigma(i) = max(sigma(i), 0.001) + sigma(i) = min(sigma(i), 0.999) + endif + enddo +! +!--- compute scale-aware function based on Arakawa & Wu (2013) +! + do i = 1, im + if(cnvflg(i)) then + if (sigma(i) > ra1(i)) then + scaldfunc(i) = (1.-sigma(i)) * (1.-sigma(i)) + scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) + else + scaldfunc(i) = 1.0 + endif + endif + enddo +! +! final scale-aware downdraft mass flux +! + do k = kmscu, 1, -1 + do i = 1, im + if(cnvflg(i) .and. + & (k >= mrad(i) .and. k < krad(i))) then + xmfd(i,k) = scaldfunc(i) * xmfd(i,k) + dz = zl(i,k+1) - zl(i,k) + xmmx = dz / dt2 + xmfd(i,k) = min(xmfd(i,k),xmmx) + endif + enddo + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute downdraft property using updated entranment rate +! + do i = 1, im + if(cnvflg(i)) then + k = krad(i) + thld(i,k)= thlx(i,k) + endif + enddo +! +! do i = 1, im +! if(cnvflg(i)) then +! k = krad(i) +! ptem1 = max(qcdo(i,k,ntcw), 0.) +! tld = thld(i,k) / pix(i,k) +! tcdo(i,k) = tld + elocp * ptem1 +! qcdo(i,k,1) = qcdo(i,k,1)+0.2*qcdo(i,k,1) +! qcdo(i,k,ntcw) = qcdo(i,k,ntcw)+0.2*qcdo(i,k,ntcw) +! endif +! enddo +! + do k = kmscu,1,-1 + do i=1,im + if(cnvflg(i) .and. + & (k >= mrad(i) .and. k < krad(i))) then + dz = zl(i,k+1) - zl(i,k) + tem = 0.5 * xlamde(i,k) * dz + factor = 1. + tem +! + thld(i,k) = ((1.-tem)*thld(i,k+1)+tem* + & (thlx(i,k)+thlx(i,k+1)))/factor + qtd(i,k) = ((1.-tem)*qtd(i,k+1)+tem* + & (qtx(i,k)+qtx(i,k+1)))/factor +! + tld = thld(i,k) / pix(i,k) + es = 0.01 * fpvs(tld) ! fpvs in pa + qs = max(qmin, eps * es / (plyr(i,k)+epsm1*es)) + dq = qtd(i,k) - qs +! + if (dq > 0.) then + gamma = el2orc * qs / (tld**2) + qld = dq / (1. + gamma) + qtd(i,k) = qs + qld + qcdo(i,k,1) = qs + qcdo(i,k,ntcw) = qld + tcdo(i,k) = tld + elocp * qld + else + qcdo(i,k,1) = qtd(i,k) + qcdo(i,k,ntcw) = 0. + tcdo(i,k) = tld + endif +! + endif + enddo + enddo +! + do k = kmscu, 1, -1 + do i = 1, im + if (cnvflg(i) .and. k < krad(i)) then + if(k >= mrad(i)) then + dz = zl(i,k+1) - zl(i,k) + tem = 0.5 * xlamdem(i,k) * dz + factor = 1. + tem + ptem = tem - pgcon + ptem1= tem + pgcon +! + ucdo(i,k) = ((1.-tem)*ucdo(i,k+1)+ptem*u1(i,k+1) + & +ptem1*u1(i,k))/factor + vcdo(i,k) = ((1.-tem)*vcdo(i,k+1)+ptem*v1(i,k+1) + & +ptem1*v1(i,k))/factor + endif + endif + enddo + enddo +! + if(ntcw > 2) then +! + do n = 2, ntcw-1 + do k = kmscu, 1, -1 + do i = 1, im + if (cnvflg(i) .and. k < krad(i)) then + if(k >= mrad(i)) then + dz = zl(i,k+1) - zl(i,k) + tem = 0.5 * xlamde(i,k) * dz + factor = 1. + tem +! + qcdo(i,k,n) = ((1.-tem)*qcdo(i,k+1,n)+tem* + & (q1(i,k,n)+q1(i,k+1,n)))/factor + endif + endif + enddo + enddo + enddo +! + endif +! + ndc = ntrac1 - ntcw +! + if(ndc > 0) then +! + do n = ntcw+1, ntrac1 + do k = kmscu, 1, -1 + do i = 1, im + if (cnvflg(i) .and. k < krad(i)) then + if(k >= mrad(i)) then + dz = zl(i,k+1) - zl(i,k) + tem = 0.5 * xlamde(i,k) * dz + factor = 1. + tem +! + qcdo(i,k,n) = ((1.-tem)*qcdo(i,k+1,n)+tem* + & (q1(i,k,n)+q1(i,k+1,n)))/factor + endif + endif + enddo + enddo + enddo +! + endif +! + return + end diff --git a/gfsphysics/physics/samfdeepcnv.f b/gfsphysics/physics/samfdeepcnv.f index eb86428bf..20569e0a4 100755 --- a/gfsphysics/physics/samfdeepcnv.f +++ b/gfsphysics/physics/samfdeepcnv.f @@ -6,6 +6,8 @@ !! !! The SAMF scheme updates the SAS scheme with scale- and aerosol-aware parameterizations from Han et al. (2017) \cite han_et_al_2017 based on the studies by Arakawa and Wu (2013) \cite arakawa_and_wu_2013 and Grell and Freitas (2014) \cite grell_and_freitus_2014 for scale awareness and by Lim (2011) \cite lim_2011 for aerosol awareness. The ratio of advective time to convective turnover time is also taken into account for the scale-aware parameterization. Along with the scale- and aerosol-aware parameterizations, more changes are made to the SAMF scheme. The cloud base mass-flux computation is modified to use convective turnover time as the convective adjustment time scale. The rain conversion rate is modified to decrease with decreasing air temperature above the freezing level. Convective inhibition in the sub-cloud layer is used as an additional trigger condition. Convective cloudiness is enhanced by considering suspended cloud condensate in the updraft. The lateral entrainment is also enhanced to more strongly suppress convection in a drier environment. !! +!! In further update for FY19 GFS implementation, interaction with turbulent kinetic energy (TKE), which is a prognostic variable used in a scale-aware TKE-based moist EDMF vertical turbulent mixing scheme, is included. Entrainment rates in updrafts and downdrafts are proportional to sub-cloud mean TKE. TKE is transported by cumulus convection. TKE contribution from cumulus convection is deduced from cumulus mass flux. On the other hand, tracers such as ozone and aerosol are also transported by cumulus convection. +!! !! \section diagram Calling Hierarchy Diagram !! \image html SAMF_Flowchart.png "Diagram depicting how the SAMF deep convection scheme is called from the FV3GFS physics time loop" height=2cm !! \section intraphysics Intraphysics Communication @@ -24,10 +26,13 @@ !! \param[in] ix horizontal dimension !! \param[in] km vertical layer dimension !! \param[in] delt physics time step in seconds +!! \param[in] ntk index for TKE +!! \param[in] ntr total number of tracers including TKE !! \param[in] delp pressure difference between level k and k+1 (Pa) !! \param[in] prslp mean layer presure (Pa) !! \param[in] psp surface pressure (Pa) !! \param[in] phil layer geopotential (\f$m^2/s^2\f$) +!! \param[in] qtr tracer array including cloud condensate (\f$kg/kg\f$) !! \param[inout] ql cloud water or ice (kg/kg) !! \param[inout] q1 updated tracers (kg/kg) !! \param[inout] t1 updated temperature (K) @@ -72,8 +77,9 @@ !! !! \section detailed Detailed Algorithm !! @{ - subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, - & q1,t1,u1,v1,cldwrk,rn,kbot,ktop,kcnv,islimsk,garea, + subroutine samfdeepcnv(im,ix,km,delt,ntk,ntr,delp, + & prslp,psp,phil,qtr,q1,t1,u1,v1, + & cldwrk,rn,kbot,ktop,kcnv,islimsk,garea, & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & clam,c0s,c1,betal,betas,evfact,evfactl,pgcon,asolfac) ! @@ -85,34 +91,37 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, &, eps => con_eps, epsm1 => con_epsm1 implicit none ! - integer, intent(in) :: im, ix, km, ncloud + integer, intent(in) :: im, ix, km, ntk, ntr, ncloud integer, intent(in) :: islimsk(im) real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(im), delp(ix,km), & prslp(ix,km), garea(im), dot(ix,km), phil(ix,km) integer, intent(inout) :: kcnv(im) - real(kind=kind_phys), intent(inout) :: ql(ix,km,2), + real(kind=kind_phys), intent(inout) :: qtr(ix,km,ntr+2), & q1(ix,km), t1(ix,km), u1(ix,km), v1(ix,km) integer, intent(out) :: kbot(im), ktop(im) real(kind=kind_phys), intent(out) :: cldwrk(im), & rn(im), cnvw(ix,km), cnvc(ix,km), & ud_mf(im,km),dd_mf(im,km), dt_mf(im,km) + + real(kind=kind_phys) clam, c0s, c1, + & betal, betas, asolfac, + & evfact, evfactl, pgcon ! !------local variables integer i, indx, jmn, k, kk, km1, n ! integer latd,lond ! - real(kind=kind_phys) clam, cxlamu, cxlamd, + real(kind=kind_phys) clamd, tkemx, tkemn, dtke, + & beta, dbeta, betamx, betamn, + & cxlame, cxlamd, & xlamde, xlamdd, & crtlamu, crtlamd ! ! real(kind=kind_phys) detad - real(kind=kind_phys) adw, aup, aafac, - & beta, betal, betas, - & c0l, c0s, d0, - & c1, asolfac, + real(kind=kind_phys) adw, aup, aafac, d0, & dellat, delta, desdt, dg, & dh, dhh, dp, & dq, dqsdp, dqsdt, dt, @@ -124,8 +133,7 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & edtmaxl, edtmaxs, el2orc, elocp, & es, etah, & cthk, dthk, - & evef, evfact, evfactl, fact1, - & fact2, factor, + & evef, fact1, fact2, factor, & g, gamma, pprime, cm, & qlk, qrch, qs, & rain, rfact, shear, tfac, @@ -137,8 +145,7 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & xdby, xpw, xpwd, ! & xqrch, mbdt, tem, & xqrch, tem, tem1, tem2, - & ptem, ptem1, ptem2, - & pgcon + & ptem, ptem1, ptem2 ! integer kb(im), kbcon(im), kbcon1(im), & ktcon(im), ktcon1(im), ktconn(im), @@ -146,8 +153,8 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & kbm(im), kmax(im) ! ! real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im), - real(kind=kind_phys) aa1(im), - & ps(im), del(ix,km), prsl(ix,km), + real(kind=kind_phys) aa1(im), tkemean(im),clamt(im), + & ps(im), del(ix,km), prsl(ix,km), & umean(im), tauadv(im), gdx(im), & delhbar(im), delq(im), delq2(im), & delqbar(im), delqev(im), deltbar(im), @@ -155,14 +162,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & edto(im), edtx(im), fld(im), & hcdo(im,km), hmax(im), hmin(im), & ucdo(im,km), vcdo(im,km),aa2(im), + & ecdo(im,km,ntr), & pdot(im), po(im,km), & pwavo(im), pwevo(im), mbdt(im), & qcdo(im,km), qcond(im), qevap(im), & rntot(im), vshear(im), xaa0(im), - & xk(im), xlamd(im), cina(im), + & xlamd(im), xk(im), cina(im), & xmb(im), xmbmax(im), xpwav(im), - & xpwev(im), xlamx(im), - & delubar(im),delvbar(im) + & xpwev(im), xlamx(im), delebar(im,ntr), + & delubar(im), delvbar(im) ! real(kind=kind_phys) c0(im) cj @@ -192,6 +200,9 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, ! parameter(cm=1.0,delta=fv) parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) + parameter(clamd=0.03,tkemx=0.65,tkemn=0.05) + parameter(dtke=tkemx-tkemn) + parameter(dbeta=0.1) parameter(cthk=200.,dthk=25.) parameter(cinpcrmx=180.,cinpcrmn=120.) ! parameter(cinacrmx=-120.,cinacrmn=-120.) @@ -201,7 +212,8 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, ! ! local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), - & uo(im,km), vo(im,km), qeso(im,km) + & uo(im,km), vo(im,km), qeso(im,km), + & ctr(im,km,ntr), ctro(im,km,ntr) ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) @@ -214,8 +226,10 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & fent1(im,km), fent2(im,km), frh(im,km), & heo(im,km), heso(im,km), & qrcd(im,km), dellah(im,km), dellaq(im,km), + & dellae(im,km,ntr), & dellau(im,km), dellav(im,km), hcko(im,km), & ucko(im,km), vcko(im,km), qcko(im,km), + & ecko(im,km,ntr), & eta(im,km), etad(im,km), zi(im,km), & qrcko(im,km), qrcdo(im,km), & pwo(im,km), pwdo(im,km), c0t(im,km), @@ -343,15 +357,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, ! betas = .15 ! betal = .05 ! betas = .05 -c evef = 0.07 +! evef = 0.07 ! evfact = 0.3 ! evfactl = 0.3 ! crtlamu = 1.0e-4 crtlamd = 1.0e-4 ! - cxlamu = 1.0e-3 - cxlamd = 1.0e-4 + cxlame = 1.0e-3 + cxlamd = 1.0e-3 xlamde = 1.0e-4 xlamdd = 1.0e-4 ! @@ -400,12 +414,10 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, zo(i,k) = phil(i,k) / g enddo enddo -!> - Calculate interface height and the initial entrainment rate as an inverse function of height. +!> - Calculate interface height do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) - xlamue(i,k) = clam / zi(i,k) -! xlamue(i,k) = max(xlamue(i,k), crtlamu) enddo enddo c @@ -450,6 +462,23 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo +! +! initialize tracer variables +! + do n = 3, ntr+2 + kk = n-2 + do k = 1, km + do i = 1, im + if (k <= kmax(i)) then + ctr(i,k,kk) = qtr(i,k,n) + ctro(i,k,kk) = qtr(i,k,n) + ecko(i,k,kk) = 0. + ecdo(i,k,kk) = 0. + endif + enddo + enddo + enddo +! !> - Calculate saturation specific humidity and enforce minimum moisture values. do k = 1, km do i=1,im @@ -534,7 +563,8 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, val2 = 1.e-10 qo(i,k) = max(qo(i,k), val2 ) ! qo(i,k) = min(qo(i,k),qeso(i,k)) - frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), 1.) + tem = min(qo(i,k)/qeso(i,k), 1.) + frh(i,k) = 1. - tem heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + & cp * to(i,k) + hvap * qo(i,k) heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + @@ -544,6 +574,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 1, km1 + do i=1,im + if (k <= kmax(i)-1) then + ctro(i,k,n) = .5 * (ctro(i,k,n) + ctro(i,k+1,n)) + endif + enddo + enddo + enddo c c look for the level of free convection as cloud base c @@ -625,6 +664,67 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, enddo if(totflg) return !! +! +! turbulent entrainment rate assumed to be proportional +! to subcloud mean TKE +! + if(ntk > 0) then +! + do i= 1, im + if(cnvflg(i)) then + sumx(i) = 0. + tkemean(i) = 0. + endif + enddo + do k = 1, km1 + do i = 1, im + if(cnvflg(i)) then + if(k >= kb(i) .and. k < kbcon(i)) then + dz = zo(i,k+1) - zo(i,k) + tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk)) + tkemean(i) = tkemean(i) + tem * dz + sumx(i) = sumx(i) + dz + endif + endif + enddo + enddo +! + do i= 1, im + if(cnvflg(i)) then + tkemean(i) = tkemean(i) / sumx(i) + if(tkemean(i) > tkemx) then + clamt(i) = clam + clamd + else if(tkemean(i) < tkemn) then + clamt(i) = clam - clamd + else + tem = tkemx - tkemean(i) + tem1 = 1. - 2. * tem / dtke + clamt(i) = clam + clamd * tem1 + endif + endif + enddo +! + else +! + do i= 1, im + if(cnvflg(i)) then + clamt(i) = clam + endif + enddo +! + endif +! +! also initially assume updraft entrainment rate +! is an inverse function of height +! + do k = 1, km1 + do i=1,im + if(cnvflg(i)) then + xlamue(i,k) = clamt(i) / zi(i,k) +! xlamue(i,k) = max(xlamue(i,k), crtlamu) + endif + enddo + enddo c c assume that updraft entrainment rate above cloud base is c same as that at cloud base @@ -660,7 +760,8 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, enddo enddo c -c functions rapidly decreasing with height, mimicking a cloud ensemble +c entrainment functions decreasing with height (fent), +c mimicking a cloud ensemble c (Bechtold et al., 2008) c do k = 2, km1 @@ -674,15 +775,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, enddo enddo c -c final entrainment and detrainment rates as the sum of turbulent part and -c organized entrainment depending on the environmental relative humidity -c (Bechtold et al., 2008) +c final entrainment rate as the sum of turbulent part and +c organized one depending on the environmental relative humidity +c (Bechtold et al., 2008; Derbyshire et al., 2011) c do k = 2, km1 do i=1,im if(cnvflg(i) .and. & (k > kbcon(i) .and. k < kmax(i))) then - tem = cxlamu * frh(i,k) * fent2(i,k) + tem = cxlame * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem ! tem1 = cxlamd * frh(i,k) ! xlamud(i,k) = xlamud(i,k) + tem1 @@ -747,6 +848,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, pwavo(i) = 0. endif enddo +! for tracers + do n = 1, ntr + do i = 1, im + if(cnvflg(i)) then + indx = kb(i) + ecko(i,indx,n) = ctro(i,indx,n) + endif + enddo + enddo c c cloud property is modified by the entrainment process c @@ -777,6 +887,21 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < kmax(i)) then + dz = zi(i,k) - zi(i,k-1) + tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + factor = 1. + tem + ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem* + & (ctro(i,k,n)+ctro(i,k-1,n)))/factor + endif + endif + enddo + enddo + enddo c c taking account into convection inhibition due to existence of c dry layers below cloud base @@ -1375,8 +1500,21 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, enddo enddo do i = 1, im - beta = betas - if(islimsk(i) == 1) beta = betal + betamn = betas + if(islimsk(i) == 1) betamn = betal + if(ntk > 0) then + betamx = betamn + dbeta + if(tkemean(i) > tkemx) then + beta = betamn + else if(tkemean(i) < tkemn) then + beta = betamx + else + tem = (betamx - betamn) * (tkemean(i) - tkemn) + beta = betamx - tem / dtke + endif + else + beta = betamn + endif if(cnvflg(i)) then dz = (sumx(i)+zi(i,1))/float(kbcon(i)) tem = 1./float(kbcon(i)) @@ -1417,6 +1555,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, pwevo(i) = 0. endif enddo +! for tracers + do n = 1, ntr + do i = 1, im + if(cnvflg(i)) then + jmn = jmin(i) + ecdo(i,jmn,n) = ctro(i,jmn,n) + endif + enddo + enddo cj !> - Calculate the cloud properties as a parcel descends, modified by entrainment and detrainment. Discretization follows Appendix B of Grell (1993) \cite grell_1993 . do k = km1, 1, -1 @@ -1446,6 +1593,19 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = km1, 1, -1 + do i = 1, im + if (cnvflg(i) .and. k < jmin(i)) then + dz = zi(i,k+1) - zi(i,k) + tem = 0.5 * xlamde * dz + factor = 1. + tem + ecdo(i,k,n) = ((1.-tem)*ecdo(i,k+1,n)+tem* + & (ctro(i,k,n)+ctro(i,k+1,n)))/factor + endif + enddo + enddo + enddo c !> - Compute the amount of moisture that is necessary to keep the downdraft saturated. do k = km1, 1, -1 @@ -1548,6 +1708,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 1, km + do i = 1, im + if(cnvflg(i) .and. k <= kmax(i)) then + dellae(i,k,n) = 0. + endif + enddo + enddo + enddo do i = 1, im if(cnvflg(i)) then dp = 1000. * del(i,1) @@ -1561,6 +1730,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & - vo(i,1)) * g / dp endif enddo + do n = 1, ntr + do i = 1, im + if(cnvflg(i)) then + dp = 1000. * del(i,1) + dellae(i,1,n) = edto(i) * etad(i,1) * (ecdo(i,1,n) + & - ctro(i,1,n)) * g / dp + endif + enddo + enddo c c--- changed due to subsidence and entrainment c @@ -1625,6 +1803,27 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 2, km1 + do i = 1, im + if (cnvflg(i) .and. k < ktcon(i)) then + aup = 1. + if(k <= kb(i)) aup = 0. + adw = 1. + if(k > jmin(i)) adw = 0. + dp = 1000. * del(i,k) +cj + tem1=eta(i,k)*(ctro(i,k,n)-ecko(i,k,n)) + tem2=eta(i,k-1)*(ctro(i,k-1,n)-ecko(i,k-1,n)) + ptem1=etad(i,k)*(ctro(i,k,n)-ecdo(i,k,n)) + ptem2=etad(i,k-1)*(ctro(i,k-1,n)-ecdo(i,k-1,n)) + dellae(i,k,n) = dellae(i,k,n) + + & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*g/dp +cj + endif + enddo + enddo + enddo c c------- cloud top c @@ -1649,6 +1848,16 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & qlko_ktcon(i) * g / dp endif enddo + do n = 1, ntr + do i = 1, im + if(cnvflg(i)) then + indx = ktcon(i) + dp = 1000. * del(i,indx) + dellae(i,indx,n) = eta(i,indx-1) * + & (ecko(i,indx-1,n) - ctro(i,indx-1,n)) * g / dp + endif + enddo + enddo c c------- final changed variable per unit mass flux c @@ -2157,6 +2366,15 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 1, km + do i = 1, im + if (cnvflg(i) .and. k <= kmax(i)) then + ctro(i,k,n) = ctr(i,k,n) + endif + enddo + enddo + enddo c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c--- feedback: simply the changes from the cloud with unit mass flux @@ -2175,6 +2393,11 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, delvbar(i) = 0. qcond(i) = 0. enddo + do n = 1, ntr + do i = 1, im + delebar(i,n) = 0. + enddo + enddo do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then @@ -2197,6 +2420,20 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + kk = n+2 + do k = 1, km + do i = 1, im + if (cnvflg(i) .and. k <= kmax(i)) then + if(k <= ktcon(i)) then + ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2 + delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/g + qtr(i,k,kk) = ctr(i,k,n) + endif + endif + enddo + enddo + enddo !> - Recalculate saturation specific humidity using the updated temperature. do k = 1, km do i = 1, im @@ -2355,11 +2592,11 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, if (k >= kbcon(i) .and. k <= ktcon(i)) then tem = dellal(i,k) * xmb(i) * dt2 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf)) - if (ql(i,k,2) > -999.0) then - ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice - ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water + if (qtr(i,k,2) > -999.0) then + qtr(i,k,1) = qtr(i,k,1) + tem * tem1 ! ice + qtr(i,k,2) = qtr(i,k,2) + tem *(1.0-tem1) ! water else - ql(i,k,1) = ql(i,k,1) + tem + qtr(i,k,1) = qtr(i,k,1) + tem endif endif endif @@ -2381,6 +2618,19 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + kk = n+2 + do k = 1, km + do i = 1, im + if(cnvflg(i) .and. rn(i) <= 0.) then + if (k <= kmax(i)) then + ctr(i,k,n)= ctro(i,k,n) + qtr(i,k,kk)= ctr(i,k,n) + endif + endif + enddo + enddo + enddo ! ! hchuang code change ! @@ -2413,6 +2663,40 @@ subroutine samfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo +! +! include TKE contribution from deep convection +! + if (ntk > 0) then +! + do k = 2, km1 + do i = 1, im + if(cnvflg(i) .and. rn(i) > 0.) then + if(k > kb(i) .and. k < ktop(i)) then + tem = 0.5 * (eta(i,k-1) + eta(i,k)) * xmb(i) + tem1 = pfld(i,k) * 100. / (rd * t1(i,k)) + sigmagfm(i) = max(sigmagfm(i), betaw) + ptem = tem / (sigmagfm(i) * tem1) + qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*sigmagfm(i)*ptem*ptem + endif + endif + enddo + enddo +! + do k = 2, km1 + do i = 1, im + if(cnvflg(i) .and. rn(i) > 0.) then + if(k > 1 .and. k <= jmin(i)) then + tem = 0.5*edto(i)*(etad(i,k-1)+etad(i,k))*xmb(i) + tem1 = pfld(i,k) * 100. / (rd * t1(i,k)) + sigmagfm(i) = max(sigmagfm(i), betaw) + ptem = tem / (sigmagfm(i) * tem1) + qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*sigmagfm(i)*ptem*ptem + endif + endif + enddo + enddo +! + endif !! return end diff --git a/gfsphysics/physics/samfshalcnv.f b/gfsphysics/physics/samfshalcnv.f index b3b9672ab..8cdb9364e 100755 --- a/gfsphysics/physics/samfshalcnv.f +++ b/gfsphysics/physics/samfshalcnv.f @@ -4,6 +4,8 @@ !! !! The previous version of the shallow convection scheme (shalcnv.f) is described in Han and Pan (2011) \cite han_and_pan_2011 and differences between the shallow and deep convection schemes are presented in Han and Pan (2011) \cite han_and_pan_2011 and Han et al. (2017) \cite han_et_al_2017 . Details of scale- and aerosol-aware parameterizations are described in Han et al. (2017) \cite han_et_al_2017 . !! +!! In further update for FY19 GFS implementation, interaction with turbulent kinetic energy (TKE), which is a prognostic variable used in a scale-aware TKE-based moist EDMF vertical turbulent mixing scheme, is included. Entrainment rates in updrafts are proportional to sub-cloud mean TKE. TKE is transported by cumulus convection. TKE contribution from cumulus convection is deduced from cumulus mass flux. On the other hand, tracers such as ozone and aerosol are also transported by cumulus convection. +!! !! \section diagram Calling Hierarchy Diagram !! \image html SAMF_shal_Flowchart.png "Diagram depicting how the SAMF shallow convection scheme is called from the FV3GFS physics time loop" height=2cm !! \section intraphysics Intraphysics Communication @@ -19,12 +21,14 @@ !! \param[in] im number of used points !! \param[in] ix horizontal dimension !! \param[in] km vertical layer dimension -!! \param[in] jcap number of spectral wave trancation !! \param[in] delt physics time step in seconds +!! \param[in] ntk index for TKE +!! \param[in] ntr total number of tracers including TKE !! \param[in] delp pressure difference between level k and k+1 (Pa) !! \param[in] prslp mean layer presure (Pa) !! \param[in] psp surface pressure (Pa) !! \param[in] phil layer geopotential (\f$m^s/s^2\f$) +!! \param[in] qtr tracer array including cloud condensate (\f$kg/kg\f$) !! \param[inout] ql cloud water or ice (kg/kg) !! \param[inout] q1 updated tracers (kg/kg) !! \param[inout] t1 updated temperature (K) @@ -58,10 +62,10 @@ !! -# For the "feedback control", calculate updated values of the state variables by multiplying the cloud base mass flux and the tendencies calculated per unit cloud base mass flux from the static control. !! \section detailed Detailed Algorithm !! @{ - subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, - & q1,t1,u1,v1,rn,kbot,ktop,kcnv,islimsk,garea, + subroutine samfshalcnv(im,ix,km,delt,ntk,ntr,delp, + & prslp,psp,phil,qtr,q1,t1,u1,v1, + & rn,kbot,ktop,kcnv,islimsk,garea, & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, -! & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc,me) & clam,c0s,c1,pgcon,asolfac) ! use machine , only : kind_phys @@ -72,34 +76,37 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, &, eps => con_eps, epsm1 => con_epsm1 implicit none ! - integer im, ix, km, ncloud, - & kbot(im), ktop(im), kcnv(im) -! &, me - real(kind=kind_phys) delt - real(kind=kind_phys) psp(im), delp(ix,km), prslp(ix,km) - real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), - & ql(ix,km,2),q1(ix,km), t1(ix,km), - & u1(ix,km), v1(ix,km), -! & u1(ix,km), v1(ix,km), rcs(im), - & rn(im), garea(im), - & dot(ix,km), phil(ix,km), hpbl(im), - & cnvw(ix,km),cnvc(ix,km) -! hchuang code change mass flux output - &, ud_mf(im,km),dt_mf(im,km) + integer, intent(in) :: im, ix, km, ntk, ntr, ncloud + integer, intent(in) :: islimsk(im) + real(kind=kind_phys), intent(in) :: delt + real(kind=kind_phys), intent(in) :: psp(im), delp(ix,km), + & prslp(ix,km), garea(im), hpbl(im), dot(ix,km), phil(ix,km) +! + integer, intent(inout) :: kcnv(im) + real(kind=kind_phys), intent(inout) :: qtr(ix,km,ntr+2), + & q1(ix,km), t1(ix,km), u1(ix,km), v1(ix,km) +! + integer, intent(out) :: kbot(im), ktop(im) + real(kind=kind_phys), intent(out) :: rn(im), + & cnvw(ix,km), cnvc(ix,km), ud_mf(im,km), dt_mf(im,km) +! + real(kind=kind_phys) clam, c0s, c1, + & asolfac, pgcon ! +! local variables integer i,j,indx, k, kk, km1, n integer kpbl(im) - integer, dimension(im), intent(in) :: islimsk +! + real(kind=kind_phys) clamd, tkemx, tkemn, dtke ! real(kind=kind_phys) dellat, delta, - & c0l, c0s, d0, - & c1, asolfac, + & c0l, d0, & desdt, dp, & dq, dqsdp, dqsdt, dt, & dt2, dtmax, dtmin, dxcrt, & dv1h, dv2h, dv3h, & dv1q, dv2q, dv3q, - & dz, dz1, e1, clam, + & dz, dz1, e1, & el2orc, elocp, aafac, cm, & es, etah, h1, & evef, evfact, evfactl, fact1, @@ -112,15 +119,16 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & w2l, w2s, w3, w3l, & w3s, w4, w4l, w4s, & rho, tem, tem1, tem2, - & ptem, ptem1, - & pgcon + & ptem, ptem1 ! integer kb(im), kbcon(im), kbcon1(im), & ktcon(im), ktcon1(im), ktconn(im), & kbm(im), kmax(im) ! real(kind=kind_phys) aa1(im), cina(im), - & umean(im), tauadv(im), gdx(im), + & tkemean(im), clamt(im), + & ps(im), del(ix,km), prsl(ix,km), + & umean(im), tauadv(im), gdx(im), & delhbar(im), delq(im), delq2(im), & delqbar(im), delqev(im), deltbar(im), & deltv(im), dtconv(im), edt(im), @@ -128,6 +136,7 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & qcond(im), qevap(im), hmax(im), & rntot(im), vshear(im), & xlamud(im), xmb(im), xmbmax(im), + & delebar(im,ntr), & delubar(im), delvbar(im) ! real(kind=kind_phys) c0(im) @@ -156,10 +165,12 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, ! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s) ! Nccn: CCN number concentration in cm^(-3) ! Until a realistic Nccn is provided, Nccns are assumed -! as Nccn=100 for sea and Nccn=7000 for land +! as Nccn=100 for sea and Nccn=1000 for land ! parameter(cm=1.0,delta=fv) parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) + parameter(clamd=0.1,tkemx=0.65,tkemn=0.05) + parameter(dtke=tkemx-tkemn) parameter(dthk=25.) parameter(cinpcrmx=180.,cinpcrmn=120.) ! parameter(cinacrmx=-120.,cinacrmn=-120.) @@ -171,7 +182,8 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, parameter(h1=0.33333333) c local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), - & uo(im,km), vo(im,km), qeso(im,km) + & uo(im,km), vo(im,km), qeso(im,km), + & ctr(im,km,ntr), ctro(im,km,ntr) ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) @@ -182,9 +194,11 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & dbyo(im,km), zo(im,km), xlamue(im,km), & heo(im,km), heso(im,km), & dellah(im,km), dellaq(im,km), + & dellae(im,km,ntr), & dellau(im,km), dellav(im,km), hcko(im,km), & ucko(im,km), vcko(im,km), qcko(im,km), - & qrcko(im,km), eta(im,km), + & qrcko(im,km), ecko(im,km,ntr), + & eta(im,km), & zi(im,km), pwo(im,km), c0t(im,km), & sumx(im), tx1(im), cnvwt(im,km) ! @@ -324,16 +338,12 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, zo(i,k) = phil(i,k) / g enddo enddo -!> - Calculate interface height and the entrainment rate as an inverse function of height. +!> - Calculate interface height do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) - xlamue(i,k) = clam / zi(i,k) enddo enddo - do i=1,im - xlamue(i,km) = xlamue(i,km1) - enddo c c pbl height c @@ -385,6 +395,21 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo +! +! initialize tracer variables +! + do n = 3, ntr+2 + kk = n-2 + do k = 1, km + do i = 1, im + if (cnvflg(i) .and. k <= kmax(i)) then + ctr(i,k,kk) = qtr(i,k,n) + ctro(i,k,kk) = qtr(i,k,n) + ecko(i,k,kk) = 0. + endif + enddo + enddo + enddo !> - Calculate saturation specific humidity and enforce minimum moisture values. do k = 1, km do i=1,im @@ -480,6 +505,15 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 1, km1 + do i=1,im + if (cnvflg(i) .and. k <= kmax(i)-1) then + ctro(i,k,n) = .5 * (ctro(i,k,n) + ctro(i,k+1,n)) + endif + enddo + enddo + enddo c c look for the level of free convection as cloud base c @@ -562,7 +596,73 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, totflg = totflg .and. (.not. cnvflg(i)) enddo if(totflg) return +! +! turbulent entrainment rate assumed to be proportional +! to subcloud mean TKE +! + if(ntk > 0) then +! + do i= 1, im + if(cnvflg(i)) then + sumx(i) = 0. + tkemean(i) = 0. + endif + enddo + do k = 1, km1 + do i = 1, im + if(cnvflg(i)) then + if(k >= kb(i) .and. k < kbcon(i)) then + dz = zo(i,k+1) - zo(i,k) + tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk)) + tkemean(i) = tkemean(i) + tem * dz + sumx(i) = sumx(i) + dz + endif + endif + enddo + enddo +! + do i= 1, im + if(cnvflg(i)) then + tkemean(i) = tkemean(i) / sumx(i) + if(tkemean(i) > tkemx) then + clamt(i) = clam + clamd + else if(tkemean(i) < tkemn) then + clamt(i) = clam - clamd + else + tem = tkemx - tkemean(i) + tem1 = 1. - 2. * tem / dtke + clamt(i) = clam + clamd * tem1 + endif + endif + enddo +! + else +! + do i= 1, im + if(cnvflg(i)) then + clamt(i) = clam + endif + enddo +! + endif !! +! +! assume updraft entrainment rate +! is an inverse function of height +! + do k = 1, km1 + do i=1,im + if(cnvflg(i)) then + xlamue(i,k) = clamt(i) / zi(i,k) +! xlamue(i,k) = max(xlamue(i,k), crtlamu) + endif + enddo + enddo + do i=1,im + if(cnvflg(i)) then + xlamue(i,km) = xlamue(i,km1) + endif + enddo c c specify the detrainment rate for the updrafts c @@ -627,6 +727,15 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, vcko(i,indx) = vo(i,indx) endif enddo +! for tracers + do n = 1, ntr + do i = 1, im + if(cnvflg(i)) then + indx = kb(i) + ecko(i,indx,n) = ctro(i,indx,n) + endif + enddo + enddo c ! cm is an enhancement factor in entrainment rates for momentum ! @@ -655,6 +764,21 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < kmax(i)) then + dz = zi(i,k) - zi(i,k-1) + tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + factor = 1. + tem + ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem* + & (ctro(i,k,n)+ctro(i,k-1,n)))/factor + endif + endif + enddo + enddo + enddo c c taking account into convection inhibition due to existence of c dry layers below cloud base @@ -1174,6 +1298,15 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 1, km + do i = 1, im + if(cnvflg(i) .and. k <= kmax(i)) then + dellae(i,k,n) = 0. + endif + enddo + enddo + enddo c c--- changed due to subsidence and entrainment c @@ -1218,6 +1351,22 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < ktcon(i)) then + dp = 1000. * del(i,k) +cj + tem1=eta(i,k)*(ctro(i,k,n)-ecko(i,k,n)) + tem2=eta(i,k-1)*(ctro(i,k-1,n)-ecko(i,k-1,n)) + dellae(i,k,n) = dellae(i,k,n) + (tem1-tem2) * g/dp +cj + endif + endif + enddo + enddo + enddo c c------- cloud top c @@ -1242,7 +1391,16 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & qlko_ktcon(i) * g / dp endif enddo -! + do n = 1, ntr + do i = 1, im + if(cnvflg(i)) then + indx = ktcon(i) + dp = 1000. * del(i,indx) + dellae(i,indx,n) = eta(i,indx-1) * + & (ecko(i,indx-1,n) - ctro(i,indx-1,n)) * g / dp + endif + enddo + enddo ! ! compute convective turn-over time ! @@ -1354,6 +1512,11 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, delvbar(i) = 0. qcond(i) = 0. enddo + do n = 1, ntr + do i = 1, im + delebar(i,n) = 0. + enddo + enddo do k = 1, km do i = 1, im if (cnvflg(i)) then @@ -1376,6 +1539,20 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, endif enddo enddo + do n = 1, ntr + kk = n+2 + do k = 1, km + do i = 1, im + if (cnvflg(i) .and. k <= kmax(i)) then + if(k <= ktcon(i)) then + ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2 + delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/g + qtr(i,k,kk) = ctr(i,k,n) + endif + endif + enddo + enddo + enddo ! !> - Recalculate saturation specific humidity using the updated temperature. do k = 1, km @@ -1474,6 +1651,14 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i)) ! endif ! enddo +! do n = 1, ntr +! do i = 1, im +! if(me == 31 .and. cnvflg(i)) then +! if(cnvflg(i)) then +! print *, ' tracer delebar = ',delebar(i,n) +! endif +! enddo +! enddo cj do i = 1, im if(cnvflg(i)) then @@ -1525,11 +1710,11 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, if (k >= kbcon(i) .and. k <= ktcon(i)) then tem = dellal(i,k) * xmb(i) * dt2 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf)) - if (ql(i,k,2) > -999.0) then - ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice - ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water + if (qtr(i,k,2) > -999.0) then + qtr(i,k,1) = qtr(i,k,1) + tem * tem1 ! ice + qtr(i,k,2) = qtr(i,k,2) + tem *(1.0-tem1) ! water else - ql(i,k,1) = ql(i,k,1) + tem + qtr(i,k,1) = qtr(i,k,1) + tem endif endif endif @@ -1559,6 +1744,26 @@ subroutine samfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, dt_mf(i,k) = ud_mf(i,k) endif enddo +! +! include TKE contribution from shallow convection +! + if (ntk > 0) then +! + do k = 2, km1 + do i = 1, im + if(cnvflg(i)) then + if(k > kb(i) .and. k < ktop(i)) then + tem = 0.5 * (eta(i,k-1) + eta(i,k)) * xmb(i) + tem1 = pfld(i,k) * 100. / (rd * t1(i,k)) + sigmagfm(i) = max(sigmagfm(i), betaw) + ptem = tem / (sigmagfm(i) * tem1) + qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*sigmagfm(i)*ptem*ptem + endif + endif + enddo + enddo +! + endif !! return end diff --git a/gfsphysics/physics/satmedmfvdif.f b/gfsphysics/physics/satmedmfvdif.f new file mode 100644 index 000000000..909f6e0ac --- /dev/null +++ b/gfsphysics/physics/satmedmfvdif.f @@ -0,0 +1,1390 @@ +!!!!! ================================================================== !!!!! +! subroutine 'satmedmfvdif.f' computes subgrid vertical turbulence mixing +! using scale-aware TKE-based moist eddy-diffusion mass-flux (EDMF) parameterization +! (by Jongil Han) +! +! For the convective boundary layer, the scheme adopts +! EDMF parameterization (Siebesma et al., 2007) to take +! into account nonlocal transport by large eddies (mfpblt.f). +! +! A new mass-flux parameterization for stratocumulus-top-induced turbulence +! mixing has been introduced (previously, it was eddy diffusion form) +! [mfscu.f]. +! +! For local turbulence mixing, a TKE closure model is used. +! +!---------------------------------------------------------------------- + subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntke, + & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea, + & psk,rbsoil,zorl,u10m,v10m,fm,fh, + & tsea,heat,evap,stress,spd1,kpbl, + & prsi,del,prsl,prslk,phii,phil,delt, + & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl, + & kinver,xkzm_m,xkzm_h,xkzm_s) +! + use machine , only : kind_phys + use funcphys , only : fpvs + use physcons, grav => con_g, rd => con_rd, cp => con_cp + &, rv => con_rv, hvap => con_hvap + &, fv => con_fvirt + &, eps => con_eps, epsm1 => con_epsm1 +! + implicit none +! +!---------------------------------------------------------------------- + integer ix, im, km, ntrac, ntcw, ntke + integer kpbl(im), kinver(im) +! + real(kind=kind_phys) delt, xkzm_m, xkzm_h, xkzm_s + real(kind=kind_phys) dv(im,km), du(im,km), + & tdt(im,km), rtg(im,km,ntrac), + & u1(ix,km), v1(ix,km), + & t1(ix,km), q1(ix,km,ntrac), + & swh(ix,km), hlw(ix,km), + & xmu(im), garea(im), + & psk(ix), rbsoil(im), + & zorl(im), tsea(im), + & u10m(im), v10m(im), + & fm(im), fh(im), + & evap(im), heat(im), + & stress(im), spd1(im), + & prsi(ix,km+1), del(ix,km), + & prsl(ix,km), prslk(ix,km), + & phii(ix,km+1), phil(ix,km), + & dusfc(im), dvsfc(im), + & dtsfc(im), dqsfc(im), + & hpbl(im) +! + logical dspheat +! flag for tke dissipative heating +! +!---------------------------------------------------------------------- +!*** +!*** local variables +!*** + integer i,is,k,kk,n,km1,kmpbl,kmscu,ntrac1 + integer lcld(im),kcld(im),krad(im),mrad(im) + integer kx1(im), kpblx(im) +! + real(kind=kind_phys) tke(im,km), tkeh(im,km-1) +! + real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km), + & qlx(im,km), thetae(im,km), + & slx(im,km), svx(im,km), thlx(im,km), + & tvx(im,km), pix(im,km), + & qtx(im,km), radx(im,km-1), + & dku(im,km-1),dkt(im,km-1), dkq(im,km-1), + & cku(im,km-1),ckt(im,km-1) +! + real(kind=kind_phys) plyr(im,km), rhly(im,km), cfly(im,km), + & qstl(im,km) +! + real(kind=kind_phys) dtdz1(im), gdx(im), + & phih(im), phim(im), prn(im), + & rbdn(im), rbup(im), thermal(im), + & ustar(im), wstar(im), hpblx(im), + & ust3(im), wst3(im), + & z0(im), crb(im), + & hgamt(im), hgamq(im), + & wscale(im),vpert(im), + & zol(im), sflux(im), radj(im), + & tx1(im), tx2(im) +! + real(kind=kind_phys) radmin(im) +! + real(kind=kind_phys) zi(im,km+1), zl(im,km), zm(im,km), + & xkzo(im,km-1),xkzmo(im,km-1), + & xkzm_hx(im), xkzm_mx(im), + & rdzt(im,km-1), + & al(im,km-1), ad(im,km), au(im,km-1), + & f1(im,km), f2(im,km*(ntrac-1)) +! + real(kind=kind_phys) elm(im,km), ele(im,km), rle(im,km-1), + & diss(im,km-1),prod(im,km-1), + & bf(im,km-1), shr2(im,km-1), + & xlamue(im,km-1), xlamde(im,km-1), + & gotvx(im,km), rlam(im,km-1) +! +! variables for updrafts (thermals) +! + real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac), + & ucko(im,km), vcko(im,km), + & buou(im,km), xmf(im,km) +! +! variables for stratocumulus-top induced downdrafts +! + real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac), + & ucdo(im,km), vcdo(im,km), + & buod(im,km), xmfd(im,km) +! + logical pblflg(im), sfcflg(im), flg(im) + logical scuflg(im), pcnvflg(im) + logical mlenflg +! +! pcnvflg: true for unstable pbl +! + real(kind=kind_phys) aphi16, aphi5, + & wfac, cfac, + & gamcrt, gamcrq, sfcfrac, + & conq, cont, conw, + & dsdz2, dsdzt, dkmax, + & dsig, dt2, dtodsd, + & dtodsu, g, factor, dz, + & gocp, gravi, zol1, zolcru, + & buop, shrp, dtn, cdtn, + & prnum, prmax, prmin, prtke, + & prscu, dw2, dw2min, zk, + & elmfac, elefac, dspmax, + & alp, clwt, cql, + & f0, robn, crbmin, crbmax, + & es, qs, value, onemrh, + & cfh, gamma, elocp, el2orc, + & epsi, beta, chx, cqx, + & rdt, rdz, qmin, qlmin, + & ri, rimin, + & rbcr, rbint, tdzmin, + & rlmn, rlmx, elmx, + & ttend, utend, vtend, qtend, + & zfac, zfmin, vk, spdk2, + & tkmin, xkzinv, dspfac, xkgdx, + & zlup, zldn, bsum, + & tem, tem1, tem2, + & ptem, ptem0, ptem1, ptem2 +! + real(kind=kind_phys) ck0, ch0, ch1, ce0, rchck +! + real(kind=kind_phys) qlcr, zstblmax +! + real(kind=kind_phys) h1 +!! + parameter(gravi=1.0/grav) + parameter(g=grav) + parameter(gocp=g/cp) + parameter(cont=cp/g,conq=hvap/g,conw=1.0/g) ! for del in pa +! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa + parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) + parameter(wfac=7.0,cfac=4.0) + parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) + parameter(vk=0.4,rimin=-100.) + parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3) + parameter(rlmn=30.,rlmx=300.,elmx=300.) + parameter(prmin=0.25,prmax=4.0,prtke=1.0,prscu=0.67) + parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35) + parameter(tkmin=5.e-5,dspfac=0.5,dspmax=10.0) + parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8) + parameter(aphi5=5.,aphi16=16.) + parameter(elmfac=1.0,elefac=1.0,cql=100.) + parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=25000.) + parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.15) + parameter(h1=0.33333333) + parameter(ck0=0.4,ch0=0.4,ch1=0.2,ce0=0.7) + parameter(rchck=1.5,cdtn=25.) +! +!************************************************************************ + dt2 = delt + rdt = 1. / dt2 +! +! the code is written assuming ntke=ntrac +! if ntrac > ntke, the code needs to be modified +! + ntrac1 = ntrac - 1 + km1 = km - 1 + kmpbl = km / 2 + kmscu = km / 2 +! + do k=1,km + do i=1,im + zi(i,k) = phii(i,k) * gravi + zl(i,k) = phil(i,k) * gravi + xmf(i,k) = 0. + xmfd(i,k) = 0. + buou(i,k) = 0. + buod(i,k) = 0. + enddo + enddo + do i=1,im + zi(i,km+1) = phii(i,km+1) * gravi + enddo + do k=1,km + do i=1,im + zm(i,k) = zi(i,k+1) + enddo + enddo +! horizontal grid size + do i=1,im + gdx(i) = sqrt(garea(i)) + enddo +! + do k=1,km + do i=1,im + tke(i,k) = max(q1(i,k,ntke), tkmin) + enddo + enddo + do k=1,km1 + do i=1,im + tkeh(i,k) = 0.5 * (tke(i,k) + tke(i,k+1)) + enddo + enddo +! + do k = 1,km1 + do i=1,im + rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k)) + enddo + enddo +! +! set background diffusivities as a function of +! horizontal grid size with xkzm_h & xkzm_m for gdx >= 25km +! and 0.01 for gdx=5m, i.e., +! xkzm_hx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) +! xkzm_mx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) +! + do i=1,im + kx1(i) = 1 + tx1(i) = 1.0 / prsi(i,1) + tx2(i) = tx1(i) + if(gdx(i) >= xkgdx) then + xkzm_hx(i) = xkzm_h + xkzm_mx(i) = xkzm_m + else + tem = 1. / (xkgdx - 5.) + tem1 = (xkzm_h - 0.01) * tem + tem2 = (xkzm_m - 0.01) * tem + ptem = gdx(i) - 5. + xkzm_hx(i) = 0.01 + tem1 * ptem + xkzm_mx(i) = 0.01 + tem2 * ptem + endif + enddo + do k = 1,km1 + do i=1,im + xkzo(i,k) = 0.0 + xkzmo(i,k) = 0.0 + if (k < kinver(i)) then +! vertical background diffusivity + ptem = prsi(i,k+1) * tx1(i) + tem1 = 1.0 - ptem + tem1 = tem1 * tem1 * 10.0 + xkzo(i,k) = xkzm_hx(i) * min(1.0, exp(-tem1)) +! vertical background diffusivity for momentum + if (ptem >= xkzm_s) then + xkzmo(i,k) = xkzm_mx(i) + kx1(i) = k + 1 + else + if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k) + tem1 = 1.0 - prsi(i,k+1) * tx2(i) + tem1 = tem1 * tem1 * 5.0 + xkzmo(i,k) = xkzm_mx(i) * min(1.0, exp(-tem1)) + endif + endif + enddo + enddo +! + do i = 1,im + z0(i) = 0.01 * zorl(i) + dusfc(i) = 0. + dvsfc(i) = 0. + dtsfc(i) = 0. + dqsfc(i) = 0. + kpbl(i) = 1 + hpbl(i) = 0. + kpblx(i) = 1 + hpblx(i) = 0. + pblflg(i)= .true. + sfcflg(i)= .true. + if(rbsoil(i) > 0.) sfcflg(i) = .false. + pcnvflg(i)= .false. + scuflg(i)= .true. + if(scuflg(i)) then + radmin(i)= 0. + mrad(i) = km1 + krad(i) = 1 + lcld(i) = km1 + kcld(i) = km1 + endif + enddo +! + do k=1,km + do i=1,im + pix(i,k) = psk(i) / prslk(i,k) + theta(i,k) = t1(i,k) * pix(i,k) + qlx(i,k) = max(q1(i,k,ntcw),qlmin) + qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k) + ptem = qlx(i,k) + tem = 1.+fv*max(q1(i,k,1),qmin)-ptem + tvx(i,k) = t1(i,k) * tem + slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap * ptem + thlx(i,k) = theta(i,k) - pix(i,k) * elocp * ptem + thlvx(i,k) = thlx(i,k)*(1.+fv*qtx(i,k)) + svx(i,k) = cp * tvx(i,k) + phil(i,k) + ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin) + thetae(i,k)= theta(i,k) + ptem1 + thvx(i,k) = theta(i,k) * tem + gotvx(i,k) = g / tvx(i,k) + enddo + enddo +! +! The background vertical diffusivities in the inversion layers are limited +! to be less than or equal to xkzminv +! + do k = 1,km1 + do i=1,im + tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k) + if(tem1 > 1.e-5) then + xkzo(i,k) = min(xkzo(i,k),xkzinv) + xkzmo(i,k) = min(xkzmo(i,k),xkzinv) + endif + enddo + enddo +! +! compute an empirical cloud fraction based on +! Xu & Randall's (1996,JAS) study +! + do k = 1, km + do i = 1, im + plyr(i,k) = 0.01 * prsl(i,k) ! pa to mb (hpa) +! --- ... compute relative humidity + es = 0.01 * fpvs(t1(i,k)) ! fpvs in pa + qs = max(qmin, eps * es / (plyr(i,k) + epsm1*es)) + rhly(i,k) = max(0.0, min(1.0, max(qmin, q1(i,k,1))/qs)) + qstl(i,k) = qs + enddo + enddo +! + do k = 1, km + do i = 1, im + cfly(i,k) = 0. + clwt = 1.0e-6 * (plyr(i,k)*0.001) + if (qlx(i,k) > clwt) then + onemrh= max(1.e-10, 1.0-rhly(i,k)) + tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) + tem1 = cql / tem1 + value = max(min( tem1*qlx(i,k), 50.0), 0.0) + tem2 = sqrt(sqrt(rhly(i,k))) + cfly(i,k) = max(tem2*(1.0-exp(-value)), 0.0) + endif + enddo + enddo +! +! compute buoyancy modified by clouds +! + do k = 1, km1 + do i = 1, im + tem = 0.5 * (svx(i,k) + svx(i,k+1)) + tem1 = 0.5 * (t1(i,k) + t1(i,k+1)) + tem2 = 0.5 * (qstl(i,k) + qstl(i,k+1)) + cfh = min(cfly(i,k+1),0.5*(cfly(i,k)+cfly(i,k+1))) + alp = g / tem + gamma = el2orc * tem2 / (tem1**2) + epsi = tem1 / elocp + beta = (1. + gamma*epsi*(1.+fv)) / (1. + gamma) + chx = cfh * alp * beta + (1. - cfh) * alp + cqx = cfh * alp * hvap * (beta - epsi) + cqx = cqx + (1. - cfh) * fv * g + ptem1 = (slx(i,k+1)-slx(i,k))*rdzt(i,k) + ptem2 = (qtx(i,k+1)-qtx(i,k))*rdzt(i,k) + bf(i,k) = chx * ptem1 + cqx * ptem2 + enddo + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! + do k=1,km1 + do i=1,im + dku(i,k) = 0. + dkt(i,k) = 0. + dkq(i,k) = 0. + cku(i,k) = 0. + ckt(i,k) = 0. + tem = zi(i,k+1)-zi(i,k) + radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k)) + enddo + enddo +! + do i = 1,im + sflux(i) = heat(i) + evap(i)*fv*theta(i,1) + if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false. + enddo +! +! compute critical bulk richardson number +! + do i = 1,im + if(pblflg(i)) then +! thermal(i) = thvx(i,1) + thermal(i) = thlvx(i,1) + crb(i) = rbcr + else + thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + tem = sqrt(u10m(i)**2+v10m(i)**2) + tem = max(tem, 1.) + robn = tem / (f0 * z0(i)) + tem1 = 1.e-7 * robn + crb(i) = 0.16 * (tem1 ** (-0.18)) + crb(i) = max(min(crb(i), crbmax), crbmin) + endif + enddo +! + do i=1,im + dtdz1(i) = dt2 / (zi(i,2)-zi(i,1)) + enddo +! + do i=1,im + ustar(i) = sqrt(stress(i)) + enddo +! +! compute buoyancy (bf) and winshear square +! + do k = 1, km1 + do i = 1, im + rdz = rdzt(i,k) +! bf(i,k) = gotvx(i,k)*(thvx(i,k+1)-thvx(i,k))*rdz + dw2 = (u1(i,k)-u1(i,k+1))**2 + & + (v1(i,k)-v1(i,k+1))**2 + shr2(i,k) = max(dw2,dw2min)*rdz*rdz + enddo + enddo +! +! find pbl height based on bulk richardson number (mrf pbl scheme) +! and also for diagnostic purpose +! + do i=1,im + flg(i) = .false. + rbup(i) = rbsoil(i) + enddo +! + do k = 1, kmpbl + do i = 1, im + if(.not.flg(i)) then + rbdn(i) = rbup(i) + spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.) +! rbup(i) = (thvx(i,k)-thermal(i))* +! & (g*zl(i,k)/thvx(i,1))/spdk2 + rbup(i) = (thlvx(i,k)-thermal(i))* + & (g*zl(i,k)/thlvx(i,1))/spdk2 + kpblx(i) = k + flg(i) = rbup(i) > crb(i) + endif + enddo + enddo + do i = 1,im + if(kpblx(i) > 1) then + k = kpblx(i) + if(rbdn(i) >= crb(i)) then + rbint = 0. + elseif(rbup(i) <= crb(i)) then + rbint = 1. + else + rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i)) + endif + hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1)) + if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1 + else + hpblx(i) = zl(i,1) + kpblx(i) = 1 + endif + hpbl(i) = hpblx(i) + kpbl(i) = kpblx(i) + if(kpbl(i) <= 1) pblflg(i)=.false. + enddo +! +! compute similarity parameters +! + do i=1,im + zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) + if(sfcflg(i)) then + zol(i) = min(zol(i),-zfmin) + else + zol(i) = max(zol(i),zfmin) + endif +! + zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) + if(sfcflg(i)) then + tem = 1.0 / (1. - aphi16*zol1) + phih(i) = sqrt(tem) + phim(i) = sqrt(phih(i)) + else + phim(i) = 1. + aphi5*zol1 + phih(i) = phim(i) + endif + enddo +! + do i=1,im + if(pblflg(i)) then + if(zol(i) < zolcru) then + pcnvflg(i) = .true. + endif + wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i) + wstar(i)= wst3(i)**h1 + ust3(i) = ustar(i)**3. + wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1 + ptem = ustar(i)/aphi5 + wscale(i) = max(wscale(i),ptem) + endif + enddo +! +! compute a thermal excess +! + do i = 1,im + if(pcnvflg(i)) then + hgamt(i) = heat(i)/wscale(i) + hgamq(i) = evap(i)/wscale(i) + vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1) + vpert(i) = max(vpert(i),0.) + tem = min(cfac*vpert(i),gamcrt) + thermal(i)= thermal(i) + tem + endif + enddo +! +! enhance the pbl height by considering the thermal excess +! (overshoot pbl top) +! + do i=1,im + flg(i) = .true. + if(pcnvflg(i)) then + flg(i) = .false. + rbup(i) = rbsoil(i) + endif + enddo + do k = 2, kmpbl + do i = 1, im + if(.not.flg(i)) then + rbdn(i) = rbup(i) + spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.) + rbup(i) = (thlvx(i,k)-thermal(i))* + & (g*zl(i,k)/thlvx(i,1))/spdk2 + kpbl(i) = k + flg(i) = rbup(i) > crb(i) + endif + enddo + enddo + do i = 1,im + if(pcnvflg(i)) then + k = kpbl(i) + if(rbdn(i) >= crb(i)) then + rbint = 0. + elseif(rbup(i) <= crb(i)) then + rbint = 1. + else + rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i)) + endif + hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1)) + if(hpbl(i) < zi(i,kpbl(i))) then + kpbl(i) = kpbl(i) - 1 + endif + if(kpbl(i) <= 1) then + pcnvflg(i) = .false. + pblflg(i) = .false. + endif + endif + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! look for stratocumulus +! + do i=1,im + flg(i) = scuflg(i) + enddo + do k = 1, km1 + do i=1,im + if(flg(i).and.zl(i,k) >= zstblmax) then + lcld(i)=k + flg(i)=.false. + endif + enddo + enddo + do i = 1, im + flg(i)=scuflg(i) + enddo + do k = kmscu,1,-1 + do i = 1, im + if(flg(i) .and. k <= lcld(i)) then + if(qlx(i,k) >= qlcr) then + kcld(i)=k + flg(i)=.false. + endif + endif + enddo + enddo + do i = 1, im + if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false. + enddo +! + do i = 1, im + flg(i)=scuflg(i) + enddo + do k = kmscu,1,-1 + do i = 1, im + if(flg(i) .and. k <= kcld(i)) then + if(qlx(i,k) >= qlcr) then + if(radx(i,k) < radmin(i)) then + radmin(i)=radx(i,k) + krad(i)=k + endif + else + flg(i)=.false. + endif + endif + enddo + enddo + do i = 1, im + if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false. + if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false. + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute components for mass flux mixing by large thermals +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! + do k = 1, km + do i = 1, im + if(pcnvflg(i)) then + tcko(i,k) = t1(i,k) + ucko(i,k) = u1(i,k) + vcko(i,k) = v1(i,k) + endif + if(scuflg(i)) then + tcdo(i,k) = t1(i,k) + ucdo(i,k) = u1(i,k) + vcdo(i,k) = v1(i,k) + endif + enddo + enddo + do kk = 1, ntrac1 + do k = 1, km + do i = 1, im + if(pcnvflg(i)) then + qcko(i,k,kk) = q1(i,k,kk) + endif + if(scuflg(i)) then + qcdo(i,k,kk) = q1(i,k,kk) + endif + enddo + enddo + enddo +! + call mfpblt(im,ix,km,kmpbl,ntcw,ntrac1,dt2, + & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx, + & gdx,hpbl,kpbl,vpert,buou,xmf, + & tcko,qcko,ucko,vcko,xlamue) +! + call mfscu(im,ix,km,kmscu,ntcw,ntrac1,dt2, + & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix, + & thlx,thvx,thlvx,gdx,thetae,radj, + & krad,mrad,radmin,buod,xmfd, + & tcdo,qcdo,ucdo,vcdo,xlamde) +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute prandtl number +! + do i = 1, im + prn(i) = phih(i)/phim(i) + prn(i) = min(prn(i),prmax) + prn(i) = max(prn(i),prmin) + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute an asymtotic mixing length +! + do k = 1, km1 + do i = 1, im + zlup = 0.0 + bsum = 0.0 + mlenflg = .true. + do n = k, km1 + if(mlenflg) then + dz = zl(i,n+1) - zl(i,n) +! ptem = gotvx(i,n)*(thvx(i,n+1)-thvx(i,k))*dz + ptem = gotvx(i,n)*(thlvx(i,n+1)-thlvx(i,k))*dz + bsum = bsum + ptem + zlup = zlup + dz + if(bsum >= tke(i,k)) then + if(ptem >= 0.) then + tem2 = max(ptem, zfmin) + else + tem2 = min(ptem, -zfmin) + endif + ptem1 = (bsum - tke(i,k)) / tem2 + zlup = zlup - ptem1 * dz + zlup = max(zlup, 0.) + mlenflg = .false. + endif + endif + enddo + zldn = 0.0 + bsum = 0.0 + mlenflg = .true. + do n = k, 1, -1 + if(mlenflg) then + if(n == 1) then + dz = zl(i,1) + tem1 = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + else + dz = zl(i,n) - zl(i,n-1) +! tem1 = thvx(i,n-1) + tem1 = thlvx(i,n-1) + endif +! ptem = gotvx(i,n)*(thvx(i,k)-tem1)*dz + ptem = gotvx(i,n)*(thlvx(i,k)-tem1)*dz + bsum = bsum + ptem + zldn = zldn + dz + if(bsum >= tke(i,k)) then + if(ptem >= 0.) then + tem2 = max(ptem, zfmin) + else + tem2 = min(ptem, -zfmin) + endif + ptem1 = (bsum - tke(i,k)) / tem2 + zldn = zldn - ptem1 * dz + zldn = max(zldn, 0.) + mlenflg = .false. + endif + endif + enddo +! + tem = 0.5 * (zi(i,k+1)-zi(i,k)) + tem1 = min(tem, rlmn) +! + ptem2 = min(zlup,zldn) + rlam(i,k) = elmfac * ptem2 + rlam(i,k) = max(rlam(i,k), tem1) + rlam(i,k) = min(rlam(i,k), rlmx) +! + ptem2 = sqrt(zlup*zldn) + ele(i,k) = elefac * ptem2 + ele(i,k) = max(ele(i,k), tem1) + ele(i,k) = min(ele(i,k), elmx) +! + enddo + enddo +! + do k = 1, km1 + do i = 1, im + tem = vk * zl(i,k) + if (zol(i) < 0.) then + ptem = 1. - 100. * zol(i) + ptem1 = ptem**0.2 + zk = tem * ptem1 + elseif (zol(i) >= 1.) then + zk = tem / 3.7 + else + ptem = 1. + 2.7 * zol(i) + zk = tem / ptem + endif + elm(i,k) = zk*rlam(i,k)/(rlam(i,k)+zk) +! + dz = zi(i,k+1) - zi(i,k) + tem = max(gdx(i),dz) + elm(i,k) = min(elm(i,k), tem) + ele(i,k) = min(ele(i,k), tem) +! + enddo + enddo + do i = 1, im + elm(i,km) = elm(i,km1) + ele(i,km) = ele(i,km1) + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute eddy diffusivities +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! + do k = 1, km1 + do i = 1, im + tem = 0.5 * (elm(i,k) + elm(i,k+1)) + tem = tem * sqrt(tkeh(i,k)) + if(k < kpbl(i)) then + if(pblflg(i)) then + dku(i,k) = ck0 * tem + dkt(i,k) = dku(i,k) / prn(i) + else + dkt(i,k) = ch0 * tem + dku(i,k) = dkt(i,k) * prn(i) + endif + else + ri = max(bf(i,k)/shr2(i,k),rimin) + if(ri < 0.) then ! unstable regime + dku(i,k) = ck0 * tem + dkt(i,k) = rchck * dku(i,k) + else ! stable regime + dkt(i,k) = ch1 * tem + prnum = 1.0 + 2.1*ri + prnum = min(prnum,prmax) + dku(i,k) = dkt(i,k) * prnum + endif + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + tem1 = ck0 * tem + ptem1 = tem1 / prscu + dku(i,k) = max(dku(i,k), tem1) + dkt(i,k) = max(dkt(i,k), ptem1) + endif + endif +! + dkq(i,k) = prtke * dkt(i,k) +! + dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = max(dkt(i,k),xkzo(i,k)) + dkq(i,k) = min(dkq(i,k),dkmax) + dkq(i,k) = max(dkq(i,k),xkzo(i,k)) + dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = max(dku(i,k),xkzmo(i,k)) +! + enddo + enddo +! + do i = 1, im + if(scuflg(i)) then + k = krad(i) + tem = bf(i,k) / gotvx(i,k) + tem1 = max(tem, tdzmin) + ptem = radj(i) / tem1 + dkt(i,k) = dkt(i,k) + ptem + dku(i,k) = dku(i,k) + ptem + dkq(i,k) = dkq(i,k) + ptem + endif + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute buoyancy and shear productions of tke +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! + do k = 1, km1 + do i = 1, im + if (k == 1) then + tem = -dkt(i,1) * bf(i,1) +! if(pcnvflg(i)) then +! ptem1 = xmf(i,1) * buou(i,1) +! else + ptem1 = 0. +! endif + if(scuflg(i) .and. mrad(i) == 1) then + ptem2 = xmfd(i,1) * buod(i,1) + else + ptem2 = 0. + endif + tem = tem + ptem1 + ptem2 + buop = 0.5 * (gotvx(i,1) * sflux(i) + tem) +! + tem1 = dku(i,1) * shr2(i,1) +! + tem = (u1(i,2)-u1(i,1))*rdzt(i,1) +! if(pcnvflg(i)) then +! ptem = xmf(i,1) * tem +! ptem1 = 0.5 * ptem * (u1(i,2)-ucko(i,2)) +! else + ptem1 = 0. +! endif + if(scuflg(i) .and. mrad(i) == 1) then + ptem = ucdo(i,1)+ucdo(i,2)-u1(i,1)-u1(i,2) + ptem = 0.5 * tem * xmfd(i,1) * ptem + else + ptem = 0. + endif + ptem1 = ptem1 + ptem +! + tem = (v1(i,2)-v1(i,1))*rdzt(i,1) +! if(pcnvflg(i)) then +! ptem = xmf(i,1) * tem +! ptem2 = 0.5 * ptem * (v1(i,2)-vcko(i,2)) +! else + ptem2 = 0. +! endif + if(scuflg(i) .and. mrad(i) == 1) then + ptem = vcdo(i,1)+vcdo(i,2)-v1(i,1)-v1(i,2) + ptem = 0.5 * tem * xmfd(i,1) * ptem + else + ptem = 0. + endif + ptem2 = ptem2 + ptem +! +! tem2 = stress(i)*spd1(i)/zl(i,1) + tem2 = stress(i)*ustar(i)*phim(i)/(vk*zl(i,1)) + shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2) + else + tem1 = -dkt(i,k-1) * bf(i,k-1) + tem2 = -dkt(i,k) * bf(i,k) + tem = 0.5 * (tem1 + tem2) + if(pcnvflg(i) .and. k <= kpbl(i)) then + ptem = 0.5 * (xmf(i,k-1) + xmf(i,k)) + ptem1 = ptem * buou(i,k) + else + ptem1 = 0. + endif + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem0 = 0.5 * (xmfd(i,k-1) + xmfd(i,k)) + ptem2 = ptem0 * buod(i,k) + else + ptem2 = 0. + endif + else + ptem2 = 0. + endif + buop = tem + ptem1 + ptem2 +! + tem1 = dku(i,k-1) * shr2(i,k-1) + tem2 = dku(i,k) * shr2(i,k) + tem = 0.5 * (tem1 + tem2) + tem1 = (u1(i,k+1)-u1(i,k))*rdzt(i,k) + tem2 = (u1(i,k)-u1(i,k-1))*rdzt(i,k-1) + if(pcnvflg(i) .and. k <= kpbl(i)) then + ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2 + ptem1 = 0.5 * ptem * (u1(i,k)-ucko(i,k)) + else + ptem1 = 0. + endif + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2 + ptem2 = 0.5 * ptem0 * (ucdo(i,k)-u1(i,k)) + else + ptem2 = 0. + endif + else + ptem2 = 0. + endif + shrp = tem + ptem1 + ptem2 + tem1 = (v1(i,k+1)-v1(i,k))*rdzt(i,k) + tem2 = (v1(i,k)-v1(i,k-1))*rdzt(i,k-1) + if(pcnvflg(i) .and. k <= kpbl(i)) then + ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2 + ptem1 = 0.5 * ptem * (v1(i,k)-vcko(i,k)) + else + ptem1 = 0. + endif + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2 + ptem2 = 0.5 * ptem0 * (vcdo(i,k)-v1(i,k)) + else + ptem2 = 0. + endif + else + ptem2 = 0. + endif + shrp = shrp + ptem1 + ptem2 + endif + prod(i,k) = buop + shrp + enddo + enddo +! +!---------------------------------------------------------------------- +! first predict tke due to tke production & dissipation(diss) +! + do k = 1,km1 + do i=1,im + rle(i,k) = ce0 / ele(i,k) + enddo + enddo + kk = max(nint(dt2/cdtn), 1) + dtn = dt2 / float(kk) + do n = 1, kk + do k = 1,km1 + do i=1,im + tem = sqrt(tke(i,k)) + diss(i,k) = rle(i,k) * tke(i,k) * tem + tem1 = prod(i,k) + tke(i,k) / dtn + diss(i,k)=max(min(diss(i,k), tem1), 0.) + tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k)) + tke(i,k) = max(tke(i,k), tkmin) + enddo + enddo + enddo +! +! compute updraft & downdraft properties for tke +! + do k = 1, km + do i = 1, im + if(pcnvflg(i)) then + qcko(i,k,ntke) = tke(i,k) + endif + if(scuflg(i)) then + qcdo(i,k,ntke) = tke(i,k) + endif + enddo + enddo + do k = 2, kmpbl + do i = 1, im + if (pcnvflg(i) .and. k <= kpbl(i)) then + dz = zl(i,k) - zl(i,k-1) + tem = 0.5 * xlamue(i,k-1) * dz + factor = 1. + tem + qcko(i,k,ntke)=((1.-tem)*qcko(i,k-1,ntke)+tem* + & (tke(i,k)+tke(i,k-1)))/factor + endif + enddo + enddo + do k = kmscu, 1, -1 + do i = 1, im + if (scuflg(i) .and. k < krad(i)) then + if(k >= mrad(i)) then + dz = zl(i,k+1) - zl(i,k) + tem = 0.5 * xlamde(i,k) * dz + factor = 1. + tem + qcdo(i,k,ntke)=((1.-tem)*qcdo(i,k+1,ntke)+tem* + & (tke(i,k)+tke(i,k+1)))/factor + endif + endif + enddo + enddo +! +!---------------------------------------------------------------------- +! compute tridiagonal matrix elements for turbulent kinetic energy +! + do i=1,im + ad(i,1) = 1.0 + f1(i,1) = tke(i,1) + enddo +! + do k = 1,km1 + do i=1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig * dkq(i,k) * rdz + dsdz2 = tem1 * rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1)= 1.-al(i,k) + tem2 = dsig * rdz +! + if(pcnvflg(i) .and. k < kpbl(i)) then + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = tke(i,k) + tke(i,k+1) + ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke) + f1(i,k) = f1(i,k)-(ptem-tem)*ptem1 + f1(i,k+1) = tke(i,k+1)+(ptem-tem)*ptem2 + else + f1(i,k+1) = tke(i,k+1) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem = 0.5 * tem2 * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = tke(i,k) + tke(i,k+1) + ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke) + f1(i,k) = f1(i,k) + (ptem - tem) * ptem1 + f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2 + endif + endif +! + enddo + enddo +c +c solve tridiagonal problem for tke +c + call tridit(im,km,1,al,ad,au,f1,au,f1) +c +c recover tendency of tke +c + do k = 1,km + do i = 1,im +! f1(i,k) = max(f1(i,k), tkmin) + qtend = (f1(i,k)-q1(i,k,ntke))*rdt + rtg(i,k,ntke) = rtg(i,k,ntke)+qtend + enddo + enddo +c +c compute tridiagonal matrix elements for heat and moisture +c + do i=1,im + ad(i,1) = 1. + f1(i,1) = t1(i,1) + dtdz1(i) * heat(i) + f2(i,1) = q1(i,1,1) + dtdz1(i) * evap(i) + enddo + if(ntrac1 >= 2) then + do kk = 2, ntrac1 + is = (kk-1) * km + do i = 1, im + f2(i,1+is) = q1(i,1,kk) + enddo + enddo + endif +c + do k = 1,km1 + do i = 1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig * dkt(i,k) * rdz + dsdzt = tem1 * gocp + dsdz2 = tem1 * rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1)= 1.-al(i,k) + tem2 = dsig * rdz +! + if(pcnvflg(i) .and. k < kpbl(i)) then + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = t1(i,k) + t1(i,k+1) + ptem = tcko(i,k) + tcko(i,k+1) + f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1 + f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2 + tem = q1(i,k,1) + q1(i,k+1,1) + ptem = qcko(i,k,1) + qcko(i,k+1,1) + f2(i,k) = f2(i,k) - (ptem - tem) * ptem1 + f2(i,k+1) = q1(i,k+1,1) + (ptem - tem) * ptem2 + else + f1(i,k) = f1(i,k)+dtodsd*dsdzt + f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt + f2(i,k+1) = q1(i,k+1,1) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem = 0.5 * tem2 * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ptem = tcdo(i,k) + tcdo(i,k+1) + tem = t1(i,k) + t1(i,k+1) + f1(i,k) = f1(i,k) + (ptem - tem) * ptem1 + f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2 + tem = q1(i,k,1) + q1(i,k+1,1) + ptem = qcdo(i,k,1) + qcdo(i,k+1,1) + f2(i,k) = f2(i,k) + (ptem - tem) * ptem1 + f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2 + endif + endif + enddo + enddo +! + if(ntrac1 >= 2) then + do kk = 2, ntrac1 + is = (kk-1) * km + do k = 1, km1 + do i = 1, im + if(pcnvflg(i) .and. k < kpbl(i)) then + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + tem = dsig * rdzt(i,k) + ptem = 0.5 * tem * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem1 = qcko(i,k,kk) + qcko(i,k+1,kk) + tem2 = q1(i,k,kk) + q1(i,k+1,kk) + f2(i,k+is) = f2(i,k+is) - (tem1 - tem2) * ptem1 + f2(i,k+1+is)= q1(i,k+1,kk) + (tem1 - tem2) * ptem2 + else + f2(i,k+1+is) = q1(i,k+1,kk) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + tem = dsig * rdzt(i,k) + ptem = 0.5 * tem * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem1 = qcdo(i,k,kk) + qcdo(i,k+1,kk) + tem2 = q1(i,k,kk) + q1(i,k+1,kk) + f2(i,k+is) = f2(i,k+is) + (tem1 - tem2) * ptem1 + f2(i,k+1+is)= f2(i,k+1+is) - (tem1 - tem2) * ptem2 + endif + endif +! + enddo + enddo + enddo + endif +c +c solve tridiagonal problem for heat and moisture +c + call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2) +c +c recover tendencies of heat and moisture +c + do k = 1,km + do i = 1,im + ttend = (f1(i,k)-t1(i,k))*rdt + qtend = (f2(i,k)-q1(i,k,1))*rdt + tdt(i,k) = tdt(i,k)+ttend + rtg(i,k,1) = rtg(i,k,1)+qtend + dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend + dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend + enddo + enddo +! + if(ntrac1 >= 2) then + do kk = 2, ntrac1 + is = (kk-1) * km + do k = 1, km + do i = 1, im + qtend = (f2(i,k+is)-q1(i,k,kk))*rdt + rtg(i,k,kk) = rtg(i,k,kk)+qtend + enddo + enddo + enddo + endif +! +! add tke dissipative heating to temperature tendency +! + if(dspheat) then + do k = 1,km1 + do i = 1,im +! tem = min(diss(i,k), dspmax) +! ttend = tem / cp + ttend = diss(i,k) / cp + tdt(i,k) = tdt(i,k) + dspfac * ttend + enddo + enddo + endif +c +c compute tridiagonal matrix elements for momentum +c + do i=1,im + ad(i,1) = 1.0 + dtdz1(i) * stress(i) / spd1(i) + f1(i,1) = u1(i,1) + f2(i,1) = v1(i,1) + enddo +c + do k = 1,km1 + do i=1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig * dku(i,k) * rdz + dsdz2 = tem1*rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1)= 1.-al(i,k) + tem2 = dsig * rdz +! + if(pcnvflg(i) .and. k < kpbl(i)) then + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = u1(i,k) + u1(i,k+1) + ptem = ucko(i,k) + ucko(i,k+1) + f1(i,k) = f1(i,k) - (ptem - tem) * ptem1 + f1(i,k+1) = u1(i,k+1) + (ptem - tem) * ptem2 + tem = v1(i,k) + v1(i,k+1) + ptem = vcko(i,k) + vcko(i,k+1) + f2(i,k) = f2(i,k) - (ptem - tem) * ptem1 + f2(i,k+1) = v1(i,k+1) + (ptem - tem) * ptem2 + else + f1(i,k+1) = u1(i,k+1) + f2(i,k+1) = v1(i,k+1) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem = 0.5 * tem2 * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = u1(i,k) + u1(i,k+1) + ptem = ucdo(i,k) + ucdo(i,k+1) + f1(i,k) = f1(i,k) + (ptem - tem) *ptem1 + f1(i,k+1) = f1(i,k+1) - (ptem - tem) *ptem2 + tem = v1(i,k) + v1(i,k+1) + ptem = vcdo(i,k) + vcdo(i,k+1) + f2(i,k) = f2(i,k) + (ptem - tem) * ptem1 + f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2 + endif + endif +! + enddo + enddo +c +c solve tridiagonal problem for momentum +c + call tridi2(im,km,al,ad,au,f1,f2,au,f1,f2) +c +c recover tendencies of momentum +c + do k = 1,km + do i = 1,im + utend = (f1(i,k)-u1(i,k))*rdt + vtend = (f2(i,k)-v1(i,k))*rdt + du(i,k) = du(i,k)+utend + dv(i,k) = dv(i,k)+vtend + dusfc(i) = dusfc(i)+conw*del(i,k)*utend + dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend + enddo + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! pbl height for diagnostic purpose +! + do i = 1, im + hpbl(i) = hpblx(i) + kpbl(i) = kpblx(i) + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + return + end +! +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- + subroutine tridit(l,n,nt,cl,cm,cu,rt,au,at) +!----------------------------------------------------------------------- +cc + use machine , only : kind_phys + implicit none + integer is,k,kk,n,nt,l,i + real(kind=kind_phys) fk(l) +cc + real(kind=kind_phys) cl(l,2:n), cm(l,n), cu(l,n-1), + & rt(l,n*nt), + & au(l,n-1), at(l,n*nt), + & fkk(l,2:n-1) +c----------------------------------------------------------------------- + do i=1,l + fk(i) = 1./cm(i,1) + au(i,1) = fk(i)*cu(i,1) + enddo + do k = 1, nt + is = (k-1) * n + do i = 1, l + at(i,1+is) = fk(i) * rt(i,1+is) + enddo + enddo + do k=2,n-1 + do i=1,l + fkk(i,k) = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) + au(i,k) = fkk(i,k)*cu(i,k) + enddo + enddo + do kk = 1, nt + is = (kk-1) * n + do k=2,n-1 + do i=1,l + at(i,k+is) = fkk(i,k)*(rt(i,k+is)-cl(i,k)*at(i,k+is-1)) + enddo + enddo + enddo + do i=1,l + fk(i) = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) + enddo + do k = 1, nt + is = (k-1) * n + do i = 1, l + at(i,n+is) = fk(i)*(rt(i,n+is)-cl(i,n)*at(i,n+is-1)) + enddo + enddo + do kk = 1, nt + is = (kk-1) * n + do k=n-1,1,-1 + do i=1,l + at(i,k+is) = at(i,k+is) - au(i,k)*at(i,k+is+1) + enddo + enddo + enddo +c----------------------------------------------------------------------- + return + end