Skip to content

Bugfixes to SLC checks, post-LIAU #339

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.

Sign up for GitHub

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

Open
ClaraDraper-NOAA wants to merge 4 commits into
base: ufs/dev
Choose a base branch
from
Open

Bugfixes to SLC checks, post-LIAU #339

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3a66ae6
Fixed too-high lower bound for SMC
ClaraDraper-NOAA Jan 7, 2026
66fedab
Merge branch 'ufs/dev' into fix/stc_check
ClaraDraper-NOAA Jan 7, 2026
bbf9cff
Tidy-up for PR.
ClaraDraper-NOAA Jan 7, 2026
8415ad5
Changes to the post L-IAU checks:
ClaraDraper-NOAA Jan 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
97 changes: 25 additions & 72 deletions physics/SFC_Models/Land/Noahmp/noahmpdrv.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,20 +187,17 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, &

integer :: lsoil_incr
integer, allocatable :: mask_tile(:)
integer,allocatable :: stc_updated(:), slc_updated(:)
integer,allocatable :: stc_updated(:)
logical :: soil_freeze, soil_ice
integer :: soiltype, n_stc, n_slc
integer :: soiltype
real(kind=kind_phys) :: slc_new

integer :: i, j, ij, l, k, ib
integer :: lensfc

real(kind=kind_phys) :: smp !< for computing supercooled water
real(kind=kind_phys) :: hc_incr

integer :: nother, nsnowupd
integer :: nstcupd, nslcupd, nfrozen, nfrozen_upd
logical :: print_update_stats = .False.
real(kind=kind_phys), parameter :: tfreez_noahmp=273.16 ! tfreez used in NoahMP to determine frozen ground
Copy link
Collaborator

@grantfirl grantfirl Jan 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason why noahmp needs to redefine this constant? The host provides this value as a physical constant and it should be used to maintain consistency with the rest of physics.

Copy link
Collaborator

@grantfirl grantfirl Jan 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like con_t0c (temperature_at_zero_celsius) is set to 273.15 and con_ttp (triple_point_temperature_of_water) is set to 273.16.

Copy link
Collaborator Author

@ClaraDraper-NOAA ClaraDraper-NOAA Jan 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are really questions for the land modeling folks. @barlage do you know why Noah-MP uses the slightly higher tfreeze? I'm not keen on using con_ttp here unless Noah-MP is using 273.16 specifically because it is the triple point (as in, this would be the right answer for the wrong reason, which could create problems later on). Also, do you know whyN Noah-MP is not taking the tfreeze constant from the ufs model? If it's just a matter of it hasn't been done yet, we could look into what would be involved to update that now.

Copy link
Collaborator

@grantfirl grantfirl Jan 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a little strange that this PR is going from using con_t0c that is already consistent with the rest of physics and redefining it locally. I guess it's OK to define tfreez_noahmp, but it should be set from the host, either from con_t0c or con_ttp, IMO.

Copy link
Collaborator Author

@ClaraDraper-NOAA ClaraDraper-NOAA Jan 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The problem is that the freezing temp that is used by the NoahMP model (in the noahmp_sflx routine) is set to 273.16 in the model code locally here. con_t0c is set in the driver, but is not passed through to the NoahMP model call. The change in this PR brings the same value into the DA since the model and DA values must be the same.

An alternative to my fix is to pass the con_t0c value through to the NoahMP routine (and remove both local tfreez definitions), but this is a science change to the model. @barlage what do you think?


! --- Initialize CCPP error handling variables
errmsg = ''
Expand All @@ -223,7 +220,7 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, &
endif

if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) then
print*, "adding land iau increments "
print*, "adding land iau increments"
endif

if (Land_IAU_Control%lsoil .ne. km) then
Expand All @@ -236,11 +233,9 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, &
allocate(stc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols
allocate(slc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols
allocate(stc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny))
allocate(slc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny))

!copy background stc
stc_updated = 0
slc_updated = 0
ib = 1
do j = 1, Land_IAU_Control%ny
do k = 1, km
Expand All @@ -260,53 +255,37 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, &
lensfc = Land_IAU_Control%nx * Land_IAU_Control%ny

if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) print*,' adjusting first ', lsoil_incr, ' surface layers only, delt ', delt
! initialize variables for counts statitics to be zeros
nother = 0 ! grid cells not land
nsnowupd = 0 ! grid cells with snow (temperature not yet updated)
nstcupd = 0 ! grid cells that are updated stc
nslcupd = 0 ! grid cells that are updated slc
nfrozen = 0 ! not update as frozen soil
nfrozen_upd = 0 ! not update as frozen soil

!TODO---if only fv3 increment files are used, this can be read from file

allocate(mask_tile(lensfc))
call calculate_landinc_mask(weasd, vegtype, soiltyp, lensfc, isice_table, mask_tile)

!IAU increments are in units of 1/sec !Land_IAU_Control%dtp
!* only updating soil temp for now
dz(1) = -zsoil(1)
do k = 2, km
dz(k) = -zsoil(k) + zsoil(k-1)
enddo
!IAU increments are in units of 1/sec
ij_loop : do ij = 1, lensfc
! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land
if (mask_tile(ij) == 1) then

soil_freeze=.false.
soil_ice=.false.
do k = 1, lsoil_incr ! k = 1, km
if ( stc(ij,k) < con_t0c) soil_freeze=.true.
do k = 1, lsoil_incr
if ( stc(ij,k) < tfreez_noahmp ) soil_freeze=.true.
if ( smc(ij,k) - slc(ij,k) > 0.001 ) soil_ice=.true.

if (Land_IAU_Control%upd_stc) then
stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt !Land_IAU_Control%dtp
if (k==1) then
stc_updated(ij) = 1
nstcupd = nstcupd + 1
endif
if (k==1) stc_updated(ij) = 1
stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt
endif

if ( (stc(ij,k) < con_t0c) .and. (.not. soil_freeze) .and. (k==1) ) nfrozen_upd = nfrozen_upd + 1

! do not do updates if this layer or any above is frozen
! do not do SLC updates if this layer or any above is frozen
if ( (.not. soil_freeze ) .and. (.not. soil_ice ) ) then
if (Land_IAU_Control%upd_slc) then
if (k==1) then
nslcupd = nslcupd + 1
slc_updated(ij) = 1
endif
! apply zero limit here (higher, model-specific limits are later)
slc(ij,k) = max(slc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0)
smc(ij,k) = max(smc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0)
! if soil moisture is <0.1 mm in layer, prevent DA from further reducing it
slc(ij,k) = max(slc(ij,k) + slc_inc_flat(ij,k)*delt, min(0.0001/dz(k), slc(ij,k)))
smc(ij,k) = max(smc(ij,k) + slc_inc_flat(ij,k)*delt, min(0.0001/dz(k), smc(ij,k)))
endif
else
if (k==1) nfrozen = nfrozen+1
endif
enddo
endif ! if soil/snow point
Expand All @@ -318,27 +297,25 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, &
errmsg = 'FATAL ERROR in noahmpdrv_timestep_init: problem in set_soilveg_noahmp'
return
endif
n_stc = 0
n_slc = 0
if (Land_IAU_Control%do_stcsmc_adjustment) then
if (Land_IAU_Control%upd_stc) then
do i=1,lensfc
if (stc_updated(i) == 1 ) then ! soil-only location
n_stc = n_stc+1
soiltype = soiltyp(i)
do l = 1, lsoil_incr
if (abs(stc_inc_flat(i,l)) > Land_IAU_Control%min_T_increment) then
!the following if case applies when updated stc > melting point, it handles both
!case 1: frz ==> frz, recalculate slc, smc remains
!case 2: unfrz ==> frz, recalculate slc, smc remains
if (stc(i,l) .LT. con_t0c )then
if (stc(i,l) .LE. tfreez_noahmp )then
!recompute supercool liquid water,smc_anl remain unchanged
smp = con_hfus*(con_t0c-stc(i,l))/(con_g*stc(i,l)) !(m)
smp = con_hfus*(tfreez_noahmp-stc(i,l))/(con_g*stc(i,l)) !(m)
slc_new=maxsmc(soiltype)*(smp/satpsi(soiltype))**(-1./bb(soiltype))
slc(i,l) = max( min( slc_new, smc(i,l)), 0.0 )
endif
!case 3: frz ==> unfrz (or unfrz ==> unfrz), melt all soil ice (if any)
if (stc(i,l) .GT. con_t0c )then !do not rely on stc_bck
!case 3: frz ==> unfrz melt all soil ice (if any)
if ( stc(i,l) .GT. tfreez_noahmp .and. &
stc(i,l) - stc_inc_flat(i,l)*delt .LE. tfreez_noahmp ) then
slc(i,l)=smc(i,l)
endif
endif
Expand All @@ -347,34 +324,10 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, &
enddo
endif

if (Land_IAU_Control%upd_slc) then
dz(1) = -zsoil(1)
do l = 2, km
dz(l) = -zsoil(l) + zsoil(l-1)
enddo
do i=1,lensfc
if (slc_updated(i) == 1 ) then
n_slc = n_slc+1
! apply SM bounds (later: add upper SMC limit)
do l = 1, lsoil_incr
if (abs(slc_inc_flat(i, l)) > Land_IAU_Control%min_SLC_increment) then
! noah-mp minimum is 1 mm per layer (in SMC)
! no need to maintain frozen amount, would be v. small.
slc(i,l) = max( 0.001/dz(l), slc(i,l) )
smc(i,l) = max( 0.001/dz(l), smc(i,l) )
endif
enddo
endif
enddo
endif
endif

deallocate(stc_inc_flat, slc_inc_flat, stc_updated, slc_updated)
deallocate(mask_tile)
deallocate(stc_inc_flat, slc_inc_flat, stc_updated)
deallocate(mask_tile)

!Remove non-warning, non-error log write
!write(*,'(a,i4,a,i8)') 'noahmpdrv_timestep_init rank ', Land_IAU_Control%me, ' # of cells with stc update ', nstcupd


end subroutine noahmpdrv_timestep_init

Expand Down