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

UFS-dev PR#65 #1020

Merged
merged 13 commits into from
Jul 7, 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
20 changes: 10 additions & 10 deletions physics/module_mp_thompson.F90
Original file line number Diff line number Diff line change
Expand Up @@ -469,11 +469,11 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
end if
if (mpirank==mpiroot) then
if (is_aerosol_aware) then
write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics'
write (*,'(a)') 'Using aerosol-aware version of Thompson microphysics'
else if(merra2_aerosol_aware) then
write (0,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics'
write (*,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics'
else
write (0,'(a)') 'Using non-aerosol-aware version of Thompson microphysics'
write (*,'(a)') 'Using non-aerosol-aware version of Thompson microphysics'
end if
end if

Expand Down Expand Up @@ -896,22 +896,22 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
!> - Call table_ccnact() to read a static file containing CCN activation of aerosols. The
!! data were created from a parcel model by Feingold & Heymsfield with
!! further changes by Eidhammer and Kriedenweis
if (mpirank==mpiroot) write(0,*) ' calling table_ccnAct routine'
if (mpirank==mpiroot) write(*,*) ' calling table_ccnAct routine'
call table_ccnAct(errmsg,errflg)
if (.not. errflg==0) return

!> - Call table_efrw() and table_efsw() to creat collision efficiency table
!! between rain/snow and cloud water
if (mpirank==mpiroot) write(0,*) ' creating qc collision eff tables'
if (mpirank==mpiroot) write(*,*) ' creating qc collision eff tables'
call table_Efrw
call table_Efsw

!> - Call table_dropevap() to creat rain drop evaporation table
if (mpirank==mpiroot) write(0,*) ' creating rain evap table'
if (mpirank==mpiroot) write(*,*) ' creating rain evap table'
call table_dropEvap

!> - Call qi_aut_qs() to create conversion of some ice mass into snow category
if (mpirank==mpiroot) write(0,*) ' creating ice converting to snow table'
if (mpirank==mpiroot) write(*,*) ' creating ice converting to snow table'
call qi_aut_qs

call cpu_time(etime)
Expand Down Expand Up @@ -942,7 +942,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
call cpu_time(stime)

!> - Call qr_acr_qg() to create rain collecting graupel & graupel collecting rain table
if (mpirank==mpiroot) write(0,*) ' creating rain collecting graupel table'
if (mpirank==mpiroot) write(*,*) ' creating rain collecting graupel table'
call cpu_time(stime)
call qr_acr_qg
call cpu_time(etime)
Expand All @@ -956,7 +956,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
if (mpirank==mpiroot) print '("Computing rain collecting snow table took ",f10.3," seconds.")', etime-stime

!> - Call freezeh2o() to create cloud water and rain freezing (Bigg, 1953) table
if (mpirank==mpiroot) write(0,*) ' creating freezing of water drops table'
if (mpirank==mpiroot) write(*,*) ' creating freezing of water drops table'
call cpu_time(stime)
call freezeH2O(threads)
call cpu_time(etime)
Expand All @@ -969,7 +969,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &

endif if_not_iiwarm

if (mpirank==mpiroot) write(0,*) ' ... DONE microphysical lookup tables'
if (mpirank==mpiroot) write(*,*) ' ... DONE microphysical lookup tables'

endif if_micro_init

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