Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PBL and Convection and Microphysics update for HR2 #65

Merged
merged 10 commits into from
May 9, 2023
43 changes: 30 additions & 13 deletions physics/mfpbltq.f
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ module mfpbltq_mod
!> @{
subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
& cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
& gdx,hpbl,kpbl,vpert,buo,xmf,
& gdx,hpbl,kpbl,vpert,buo,wush,tkemean,vez0fun,xmf,
& tcko,qcko,ucko,vcko,xlamueq,a1)
!
use machine , only : kind_phys
Expand All @@ -33,7 +33,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
& 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),
& buo(im,km), wush(im,km),
& tkemean(im),vez0fun(im),xmf(im,km),
& tcko(im,km),qcko(im,km,ntrac1),
& ucko(im,km),vcko(im,km),
& xlamueq(im,km-1)
Expand All @@ -44,8 +45,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
integer kpblx(im), kpbly(im)
!
real(kind=kind_phys) dt2, dz, ce0,
& cm, cq,
& factor, gocp,
& cm, cq, tkcrt,
& factor, gocp, cmxfac,
& g, b1, f1,
& bb1, bb2,
& alp, vpertmax,a1, pgcon,
Expand All @@ -59,7 +60,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
!
real(kind=kind_phys) rbdn(im), rbup(im), hpblx(im),
& xlamue(im,km-1), xlamuem(im,km-1)
real(kind=kind_phys) delz(im), xlamax(im)
real(kind=kind_phys) delz(im), xlamax(im), ce0t(im)
!
real(kind=kind_phys) wu2(im,km), thlu(im,km),
& qtx(im,km), qtu(im,km)
Expand All @@ -73,7 +74,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
parameter(g=grav)
parameter(gocp=g/cp)
parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
parameter(ce0=0.4,cm=1.0,cq=1.3)
parameter(ce0=0.4,cm=1.0,cq=1.0,tkcrt=2.,cmxfac=5.)
parameter(qmin=1.e-8,qlmin=1.e-12)
parameter(alp=1.5,vpertmax=3.0,pgcon=0.55)
parameter(b1=0.5,f1=0.15)
Expand Down Expand Up @@ -112,13 +113,27 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
enddo
!
!> - Compute entrainment rate
!
! if tkemean>tkcrt, ce0t=sqrt(tkemean/tkcrt)*ce0
!
do i=1,im
if(cnvflg(i)) then
ce0t(i) = ce0 * vez0fun(i)
if(tkemean(i) > tkcrt) then
tem = sqrt(tkemean(i)/tkcrt)
tem1 = min(tem, cmxfac)
tem2 = tem1 * ce0
ce0t(i) = max(ce0t(i), tem2)
endif
endif
enddo
!
do i=1,im
if(cnvflg(i)) then
k = kpbl(i) / 2
k = max(k, 1)
delz(i) = zl(i,k+1) - zl(i,k)
xlamax(i) = ce0 / delz(i)
xlamax(i) = ce0t(i) / delz(i)
endif
enddo
!
Expand All @@ -129,7 +144,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
ptem = 1./(zm(i,k)+delz(i))
tem = max((hpbl(i)-zm(i,k)+delz(i)) ,delz(i))
ptem1 = 1./tem
xlamue(i,k) = ce0 * (ptem+ptem1)
xlamue(i,k) = ce0t(i) * (ptem+ptem1)
else
xlamue(i,k) = xlamax(i)
endif
Expand Down Expand Up @@ -210,11 +225,13 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
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
tem = 0.25*bb1*(xlamue(i,k-1)+xlamue(i,k))*dz
tem1 = max(wu2(i,k-1), 0.)
tem1 = bb2 * buo(i,k) - wush(i,k) * sqrt(tem1)
tem2 = tem1 * dz
ptem = (1. - tem) * wu2(i,k-1)
ptem1 = 1. + tem
wu2(i,k) = (ptem + tem1) / ptem1
wu2(i,k) = (ptem + tem2) / ptem1
endif
enddo
enddo
Expand Down Expand Up @@ -271,7 +288,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
k = kpbl(i) / 2
k = max(k, 1)
delz(i) = zl(i,k+1) - zl(i,k)
xlamax(i) = ce0 / delz(i)
xlamax(i) = ce0t(i) / delz(i)
endif
enddo
!
Expand All @@ -283,7 +300,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
ptem = 1./(zm(i,k)+delz(i))
tem = max((hpbl(i)-zm(i,k)+delz(i)) ,delz(i))
ptem1 = 1./tem
xlamue(i,k) = ce0 * (ptem+ptem1)
xlamue(i,k) = ce0t(i) * (ptem+ptem1)
else
xlamue(i,k) = xlamax(i)
endif
Expand Down
43 changes: 31 additions & 12 deletions physics/mfscuq.f
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ module mfscuq_mod
subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
& cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,
& thlx,thvx,thlvx,gdx,thetae,
& krad,mrad,radmin,buo,xmfd,
& krad,mrad,radmin,buo,wush,tkemean,vez0fun,xmfd,
& tcdo,qcdo,ucdo,vcdo,xlamdeq,a1)
!
use machine , only : kind_phys
Expand All @@ -38,9 +38,10 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
& gdx(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),
& buo(im,km), wush(im,km),
& tkemean(im),vez0fun(im),xmfd(im,km),
& tcdo(im,km),qcdo(im,km,ntrac1),
& ucdo(im,km),vcdo(im,km),
& xlamdeq(im,km-1)
!
! local variables and arrays
Expand All @@ -51,6 +52,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
!
real(kind=kind_phys) dt2, dz, ce0,
& cm, cq,
& tkcrt, cmxfac,
& gocp, factor, g, tau,
& b1, f1, bb1, bb2,
& a1, a2,
Expand All @@ -67,7 +69,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
& qtx(im,km), qtd(im,km),
& thlvd(im), hrad(im), xlamde(im,km-1),
& xlamdem(im,km-1), ra1(im)
real(kind=kind_phys) delz(im), xlamax(im)
real(kind=kind_phys) delz(im), xlamax(im), ce0t(im)
!
real(kind=kind_phys) xlamavg(im), sigma(im),
& scaldfunc(im), sumx(im)
Expand All @@ -80,7 +82,8 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
parameter(g=grav)
parameter(gocp=g/cp)
parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
parameter(ce0=0.4,cm=1.0,cq=1.3,pgcon=0.55)
parameter(ce0=0.4,cm=1.0,cq=1.0,pgcon=0.55)
parameter(tkcrt=2.,cmxfac=5.)
parameter(qmin=1.e-8,qlmin=1.e-12)
parameter(b1=0.45,f1=0.15)
parameter(a2=0.5)
Expand Down Expand Up @@ -185,13 +188,27 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
!!
!
!> - Compute entrainment rate
!
! if tkemean>tkcrt, ce0t=sqrt(tkemean/tkcrt)*ce0
!
do i=1,im
if(cnvflg(i)) then
ce0t(i) = ce0 * vez0fun(i)
if(tkemean(i) > tkcrt) then
tem = sqrt(tkemean(i)/tkcrt)
tem1 = min(tem, cmxfac)
tem2 = tem1 * ce0
ce0t(i) = max(ce0t(i), tem2)
endif
endif
enddo
!
do i=1,im
if(cnvflg(i)) then
k = mrad(i) + (krad(i)-mrad(i)) / 2
k = max(k, mrad(i))
delz(i) = zl(i,k+1) - zl(i,k)
xlamax(i) = ce0 / delz(i)
xlamax(i) = ce0t(i) / delz(i)
endif
enddo
!
Expand All @@ -206,7 +223,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
endif
tem = max((hrad(i)-zm(i,k)+delz(i)) ,delz(i))
ptem1 = 1./tem
xlamde(i,k) = ce0 * (ptem+ptem1)
xlamde(i,k) = ce0t(i) * (ptem+ptem1)
else
xlamde(i,k) = xlamax(i)
endif
Expand Down Expand Up @@ -289,10 +306,12 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
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
tem1 = max(wd2(i,k+1), 0.)
tem1 = bb2*buo(i,k+1) - wush(i,k+1)*sqrt(tem1)
tem2 = tem1 * dz
ptem = (1. - tem) * wd2(i,k+1)
ptem1 = 1. + tem
wd2(i,k) = (ptem + tem1) / ptem1
wd2(i,k) = (ptem + tem2) / ptem1
endif
enddo
enddo
Expand Down Expand Up @@ -334,7 +353,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
k = mrad(i) + (krad(i)-mrad(i)) / 2
k = max(k, mrad(i))
delz(i) = zl(i,k+1) - zl(i,k)
xlamax(i) = ce0 / delz(i)
xlamax(i) = ce0t(i) / delz(i)
endif
enddo
!
Expand All @@ -349,7 +368,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
endif
tem = max((hrad(i)-zm(i,k)+delz(i)) ,delz(i))
ptem1 = 1./tem
xlamde(i,k) = ce0 * (ptem+ptem1)
xlamde(i,k) = ce0t(i) * (ptem+ptem1)
else
xlamde(i,k) = xlamax(i)
endif
Expand Down
13 changes: 8 additions & 5 deletions physics/radiation_clouds.f
Original file line number Diff line number Diff line change
Expand Up @@ -2081,7 +2081,8 @@ subroutine progcld_thompson_wsm6 &

! --- constant values
real (kind=kind_phys), parameter :: xrc3 = 100.

real (kind=kind_phys), parameter :: snow2ice = 0.25
real (kind=kind_phys), parameter :: coef_t = 0.025
!
!===> ... begin here

Expand All @@ -2097,7 +2098,7 @@ subroutine progcld_thompson_wsm6 &
rei (i,k) = re_ice(i,k)
rer (i,k) = rrain_def ! default rain radius to 1000 micron
res (i,k) = re_snow(i,K)
tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) )
tem2d (i,k) = min(1.0, max( 0.0, (con_ttp-tlyr(i,k))*coef_t))
clwf(i,k) = 0.0
enddo
enddo
Expand Down Expand Up @@ -2138,12 +2139,14 @@ subroutine progcld_thompson_wsm6 &
if(tem1 > 1.e-12 .and. clw(i,k,ntcw) < 1.e-12)
& rew(i,k)=reliq_def
tem2 = cnvw(i,k)*tem2d(i,k)
cip(i,k) = max(0.0, (clw(i,k,ntiw) + tem2 )
& *gfac * delp(i,k))
cip(i,k) = max(0.0, (clw(i,k,ntiw) +
& snow2ice*clw(i,k,ntsw) + tem2) *
& gfac * delp(i,k))
if(tem2 > 1.e-12 .and. clw(i,k,ntiw) < 1.e-12)
& rei(i,k)=reice_def
crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k))
csp(i,k) = max(0.0, clw(i,k,ntsw) * gfac * delp(i,k))
csp(i,k) = max(0.0, (1.-snow2ice)*clw(i,k,ntsw) *
& gfac * delp(i,k))
enddo
enddo

Expand Down
Loading