Skip to content

Commit

Permalink
Merge remote-tracking branch 'haiqinli/ufs/dev-smoke-dust' into featu…
Browse files Browse the repository at this point in the history
…re/noahmp
  • Loading branch information
uturuncoglu committed Jan 24, 2024
2 parents ca1a6e3 + a0acaed commit 042d156
Show file tree
Hide file tree
Showing 16 changed files with 411 additions and 311 deletions.
5 changes: 3 additions & 2 deletions physics/CONV/Grell_Freitas/cu_gf_deep.F90
Original file line number Diff line number Diff line change
Expand Up @@ -142,13 +142,13 @@ subroutine cu_gf_deep_run( &
!! betwee -1 and +1
,do_capsuppress,cap_suppress_j & !
,k22 & !
,jmin,tropics) !
,jmin,kdt,tropics) !

implicit none

integer &
,intent (in ) :: &
nranflag,itf,ktf,its,ite, kts,kte,ipr,imid
nranflag,itf,ktf,its,ite, kts,kte,ipr,imid,kdt
integer, intent (in ) :: &
ichoice,nchem
real(kind=kind_phys), dimension (its:ite,4) &
Expand Down Expand Up @@ -591,6 +591,7 @@ subroutine cu_gf_deep_run( &
sig(i)=(1.-frh)**2
!frh_out(i) = frh
if(forcing(i,7).eq.0.)sig(i)=1.
if(kdt.le.(3600./dtime))sig(i)=1.
frh_out(i) = frh*sig(i)
enddo
!$acc end kernels
Expand Down
8 changes: 4 additions & 4 deletions physics/CONV/Grell_Freitas/cu_gf_driver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
dfi_radar_max_intervals,ldiag3d,qci_conv,do_cap_suppress, &
maxupmf,maxMF,do_mynnedmf,ichoice_in,ichoicem_in,ichoice_s_in, &
spp_cu_deep,spp_wts_cu_deep,nchem,chem3d,fscav,wetdpc_deep, &
do_smoke_transport,errmsg,errflg)
do_smoke_transport,kdt,errmsg,errflg)
!-------------------------------------------------------------
implicit none
integer, parameter :: maxiens=1
Expand All @@ -95,7 +95,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
integer :: ishallow_g3 ! depend on imfshalcnv
!-------------------------------------------------------------
integer :: its,ite, jts,jte, kts,kte
integer, intent(in ) :: im,km,ntracer, nchem
integer, intent(in ) :: im,km,ntracer,nchem,kdt
integer, intent(in ) :: ichoice_in,ichoicem_in,ichoice_s_in
logical, intent(in ) :: flag_init, flag_restart, do_mynnedmf
logical, intent(in ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend
Expand Down Expand Up @@ -766,7 +766,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
! betwee -1 and +1
,do_cap_suppress_here,cap_suppress_j &
,k22m &
,jminm,tropics)
,jminm,kdt,tropics)
!$acc kernels
do i=its,itf
do k=kts,ktf
Expand Down Expand Up @@ -853,7 +853,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
! betwee -1 and +1
,do_cap_suppress_here,cap_suppress_j &
,k22 &
,jmin,tropics)
,jmin,kdt,tropics)
jpr=0
ipr=0
!$acc kernels
Expand Down
7 changes: 7 additions & 0 deletions physics/CONV/Grell_Freitas/cu_gf_driver.meta
Original file line number Diff line number Diff line change
Expand Up @@ -651,6 +651,13 @@
type = real
kind = kind_phys
intent = inout
[kdt]
standard_name = index_of_timestep
long_name = current forecast iteration
units = index
dimensions = ()
type = integer
intent = in
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
Expand Down
146 changes: 66 additions & 80 deletions physics/PBL/MYNN_EDMF/module_bl_mynn.F90
Original file line number Diff line number Diff line change
Expand Up @@ -2001,9 +2001,9 @@ SUBROUTINE mym_length ( &
uonset= 15.
wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5))
cns = 2.7 !was 3.5
alp1 = 0.22
alp1 = 0.23
alp2 = 0.3
alp3 = 2.0 * wt_u !taper off bouyancy enhancement in shear-driven pbls
alp3 = 2.5 * wt_u !taper off bouyancy enhancement in shear-driven pbls
alp4 = 5.0
alp5 = 0.3
alp6 = 50.
Expand Down Expand Up @@ -2059,12 +2059,12 @@ SUBROUTINE mym_length ( &

! ** Length scale limited by the buoyancy effect **
IF ( dtv(k) .GT. 0.0 ) THEN
bv = max( sqrt( gtr*dtv(k) ), 0.001)
bv = max( sqrt( gtr*dtv(k) ), 0.0001)
elb = MAX(alp2*qkw(k), &
& alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv &
& *( 1.0 + alp3*SQRT( vsc/(bv*elt) ) )
elb = MIN(elb, zwk)
elf = 0.80 * qkw(k)/bv
elf = 1.0 * qkw(k)/bv
elBLavg(k) = MAX(elBLavg(k), alp6*edmf_a1(k-1)*edmf_w1(k-1)/bv)
ELSE
elb = 1.0e10
Expand All @@ -2084,8 +2084,10 @@ SUBROUTINE mym_length ( &
!add blending to use BouLac mixing length in free atmos;
!defined relative to the PBLH (zi) + transition layer (h1)
!el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf)
!try squared-blending
el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2)))
!try squared-blending - but take out elb (makes it underdiffusive)
!el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2)))
el(k) = sqrt( els**2/(1. + (els**2/elt**2)))
el(k) = min(el(k), elb)
el(k) = MIN (el(k), elf)
el(k) = el(k)*(1.-wt) + alp5*elBLavg(k)*wt

Expand Down Expand Up @@ -3633,13 +3635,13 @@ SUBROUTINE mym_condensation (kts,kte, &

real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll, &
&q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb, &
&ls,wt,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, &
&ls,wt,wt2,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, &
&qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc
real(kind_phys), parameter :: qpct_sfc=0.025
real(kind_phys), parameter :: qpct_pbl=0.030
real(kind_phys), parameter :: qpct_trp=0.040
real(kind_phys), parameter :: rhcrit =0.83 !for cloudpdf = 2
real(kind_phys), parameter :: rhmax =1.01 !for cloudpdf = 2
real(kind_phys), parameter :: rhmax =1.02 !for cloudpdf = 2
integer :: i,j,k

real(kind_phys):: erf
Expand Down Expand Up @@ -3864,25 +3866,18 @@ SUBROUTINE mym_condensation (kts,kte, &
!Add condition for falling/settling into low-RH layers, so at least
!some cloud fraction is applied for all qc, qs, and qi.
rh_hack= rh(k)
wt2 = min(max( zagl - pblh2, 0.0 )/300., 1.0)
!ensure adequate RH & q1 when qi is at least 1e-9 (above the PBLH)
if (qi(k)>1.e-9 .and. zagl .gt. pblh2) then
rh_hack =min(rhmax, rhcrit + 0.07*(9.0 + log10(qi(k))))
if ((qi(k)+qs(k))>1.e-9 .and. (zagl .gt. pblh2)) then
rh_hack =min(rhmax, rhcrit + wt2*0.045*(9.0 + log10(qi(k)+qs(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
q1(k) =max(q1_rh, q1(k) )
endif
!ensure adequate RH & q1 when qc is at least 1e-6
if (qc(k)>1.e-6) then
rh_hack =min(rhmax, rhcrit + 0.09*(6.0 + log10(qc(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
q1(k) =max(q1_rh, q1(k) )
endif
!ensure adequate RH & q1 when qs is at least 1e-8 (above the PBLH)
if (qs(k)>1.e-8 .and. zagl .gt. pblh2) then
rh_hack =min(rhmax, rhcrit + 0.07*(8.0 + log10(qs(k))))
!ensure adequate rh & q1 when qc is at least 1e-6 (above the PBLH)
if (qc(k)>1.e-6 .and. (zagl .gt. pblh2)) then
rh_hack =min(rhmax, rhcrit + wt2*0.08*(6.0 + log10(qc(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
Expand Down Expand Up @@ -3994,7 +3989,7 @@ SUBROUTINE mym_condensation (kts,kte, &
fac_damp = min(zagl * 0.0025, 1.0)
!cld_factor = 1.0 + fac_damp*MAX(0.0, ( RH(k) - 0.75 ) / 0.26 )**1.9 !HRRRv4
!cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.25 )**2, 0.3)
cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.35)
cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.37)
cldfra_bl1D(K) = min( 1., cld_factor*cldfra_bl1D(K) )
enddo

Expand Down Expand Up @@ -4181,38 +4176,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &

k=kts

!original approach (drag in b-vector):
! a(1)=0.
! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff
! c(1)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff
! d(1)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff + &
! sub_u(k)*delt + det_u(k)*delt

!rho-weighted (drag in b-vector):
a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff - &
& dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt

!rho-weighted with drag term moved out of b-array
! a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! d(k)=u(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*uoce*ust**2/wspd - &
! !!!d(k)=u(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*uoce*ust**2/wspd - &
! & dtz(k)*rhoinv(k)*s_awu(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt

DO k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + &
& 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff + &
& sub_u(k)*delt + det_u(k)*delt
ENDDO
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*uoce*ust**2/wspd &
& - dtz(k)*rhoinv(k)*s_awu(k+1)*onoff &
& + dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff &
& + sub_u(k)*delt + det_u(k)*delt

do k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
& + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+ dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
& + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff &
& - dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff &
& + sub_u(k)*delt + det_u(k)*delt
enddo

!! no flux at the top
! a(kte)=-1.
Expand Down Expand Up @@ -4247,37 +4237,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &

k=kts

!original approach (drag in b-vector):
! a(1)=0.
! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff
! c(1)= - dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff
! d(1)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff + &
! sub_v(k)*delt + det_v(k)*delt

!rho-weighted (drag in b-vector):
a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + &
& sub_v(k)*delt + det_v(k)*delt

!rho-weighted with drag term moved out of b-array
! a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! d(k)=v(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*voce*ust**2/wspd - &
! !!!d(k)=v(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*voce*ust**2/wspd - &
! & dtz(k)*rhoinv(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + sub_v(k)*delt + det_v(k)*delt

DO k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + &
& 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff + &
& sub_v(k)*delt + det_v(k)*delt
ENDDO
b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*voce*ust**2/wspd &
& - dtz(k)*rhoinv(k)*s_awv(k+1)*onoff &
& + dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff &
& + sub_v(k)*delt + det_v(k)*delt

do k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
& + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
& + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff &
& - dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff &
& + sub_v(k)*delt + det_v(k)*delt
enddo

!! no flux at the top
! a(kte)=-1.
Expand Down
6 changes: 3 additions & 3 deletions physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1687,7 +1687,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia
endif

if(newsn > zero ) then
SNOWFRACnewsn=MIN(one,SNHEI/SNHEI_CRIT_newsn)
SNOWFRACnewsn=MIN(one,snowfallac*1.e-3_kind_phys/SNHEI_CRIT_newsn)
endif

!-- due to steep slopes and blown snow, limit snow fraction in the
Expand All @@ -1700,7 +1700,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia
if(snowfrac < 0.75_kind_phys) snow_mosaic = one

KEEP_SNOW_ALBEDO = zero
IF (NEWSN > zero .and. snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN
IF (snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN
! new snow
KEEP_SNOW_ALBEDO = one
! turn off separate treatment of snow covered and snow-free portions of the grid cell
Expand Down Expand Up @@ -1735,7 +1735,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia
! hwlps with these biases..
if( snow_mosaic == one) then
ALBsn=alb_snow
if(newsn > zero .and. KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then
if(KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then
!-- Albedo correction with fresh snow and deep snow pack
!-- will reduce warm bias in western Canada
!-- and US West coast, where max snow albedo is low (0.3-0.5).
Expand Down
Loading

0 comments on commit 042d156

Please sign in to comment.