Skip to content

Commit

Permalink
Merge branch 'main' of https://github.com/NCAR/ccpp-physics into HEAD
Browse files Browse the repository at this point in the history
  • Loading branch information
dustinswales committed Aug 1, 2023
2 parents 5ab1a5b + 9b77c4e commit 07a25fe
Show file tree
Hide file tree
Showing 8 changed files with 306 additions and 175 deletions.
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

0 comments on commit 07a25fe

Please sign in to comment.