Skip to content

Robustness guards for LSTDA #50

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

Merged
jjokella merged 7 commits into from
Apr 2, 2026
Merged

Robustness guards for LSTDA #50

Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
7d06515
Added cap for LSTADA from Claude
visakhraja Mar 2, 2026
f4ee368
Merge branch 'tsmp-pdaf-patched-lstda' into tsmp-pdaf-patched-lstda-cap
jjokella Mar 3, 2026
e581021
Merge branch 'tsmp-pdaf-patched-lstda' into tsmp-pdaf-patched-lstda-cap
jjokella Mar 4, 2026
b31838a
Merge branch 'tsmp-pdaf-patched-lstda' into tsmp-pdaf-patched-lstda-cap
jjokella Mar 26, 2026
9fa68be
non-vegetated gridcells: document choice of safe fallback
jjokella Apr 2, 2026
dff2654
warning for NaN increment factors
jjokella Apr 2, 2026
8958241
apply guards to clmupdate_T==3
jjokella Apr 2, 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
100 changes: 81 additions & 19 deletions interface/model/eclm/enkf_clm_mod_5.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -946,16 +946,35 @@ subroutine set_clm_statevec_T(tstartcycle,mype)
clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = 0.0
n_p = 0
do p=clm_begp,clm_endp
if(patch%gridcell(p)==g) then
! Skip bare-ground/lake patches where t_veg retains spval (1e36).
if(patch%gridcell(p)==g .and. t_veg(p) < 1.0e20_r8) then
clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) + t_veg(p)
n_p = n_p + 1
end if
end do
if(n_p > 0) then
clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) / real(n_p, r8)
else
write(*,*) "ERROR: Gridcell g=", g, " has no patches for TVEG averaging"
error stop "Gridcell without patches in set_clm_statevec_T (TVEG)"
! No vegetated patches in this gridcell (bare ground, glacier, lake):
! t_veg retains spval (~1e36) for all patches, so no meaningful TVEG
! average can be formed. TSKIN is used as a fallback to keep the value
! physically plausible within the DA step.
!
! Note: this slot does not affect the model state. The distribution
! step guards against writing back to patches where t_veg >= 1e20
! (see update_clm_statevec_T), so the fallback value is never applied.
!
! The fallback does, however, participate in ensemble covariance
! estimation during DA, introducing an artificial TSKIN-TVEG
! correlation for these cells. A conceptually cleaner solution would
! be to exclude non-vegetated gridcells from the TVEG part of the
! state vector entirely (detectable at initialisation: all patches in
! the gridcell have t_veg >= 1e20). This would require a mask array
! and non-uniform state vector layout, breaking the current assumption
! that every gridcell contributes the same variable block. The mapping
! arrays state_pdaf2clm_p_p / state_clm2pdaf_p already handle
! non-trivial patch mappings and could serve as a template for this.
clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc)
end if

end do
Expand Down Expand Up @@ -1355,6 +1374,10 @@ subroutine update_clm_T(tstartcycle,mype)
integer :: incr_warn_count_skin, incr_warn_count_soisno, &
incr_warn_count_veg, incr_warn_count_grnd

! Guard against applying t_soisno increment multiple times to the same
! column when several patches share a column (T==2 loop is over patches).
logical, allocatable :: col_updated(:)

! LST
t_grnd => temperature_inst%t_grnd_col
t_soisno => temperature_inst%t_soisno_col
Expand All @@ -1369,6 +1392,9 @@ subroutine update_clm_T(tstartcycle,mype)
incr_warn_count_veg = 0
incr_warn_count_grnd = 0

allocate(col_updated(clm_begc:clm_endc))
col_updated = .false.

!hcp: TG, TV
if(clmupdate_T==1) then
do p = clm_begp, clm_endp
Expand Down Expand Up @@ -1401,7 +1427,10 @@ subroutine update_clm_T(tstartcycle,mype)
t_update = t_skin(p) * increment_factor
else
increment_factor = clm_statevec(cc) - clm_statevec_orig(cc)
if (abs(increment_factor) < clmT_max_increment) then
if (ieee_is_nan(increment_factor)) then
print *, "WARNING: t_skin increment_factor is NaN at p=", p, " - leaving t_skin unchanged"
t_update = t_skin(p)
else if (abs(increment_factor) < clmT_max_increment) then
t_update = t_skin(p) + increment_factor
else
t_update = t_skin(p) + sign(clmT_max_increment, increment_factor)
Expand All @@ -1411,10 +1440,12 @@ subroutine update_clm_T(tstartcycle,mype)
if (ieee_is_nan(t_update)) then
print *, "WARNING: t_skin update is NaN at p=", p
else
t_skin(p) = t_update
t_skin(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update))
end if
end if

! Guard: only update each column once even if multiple patches share it.
if (.not. col_updated(c)) then
! --- TSOIL: update with increment factor for each layer ---
do lev=1,nlevgrnd
cc = state_clm2pdaf_p(p, 1+lev)
Expand All @@ -1425,7 +1456,10 @@ subroutine update_clm_T(tstartcycle,mype)
t_update = t_soisno(c,lev) * increment_factor
else
increment_factor = clm_statevec(cc) - clm_statevec_orig(cc)
if (abs(increment_factor) < clmT_max_increment) then
if (ieee_is_nan(increment_factor)) then
print *, "WARNING: t_soisno increment_factor is NaN at c=", c, " lev=", lev, " - leaving t_soisno unchanged"
t_update = t_soisno(c,lev)
else if (abs(increment_factor) < clmT_max_increment) then
t_update = t_soisno(c,lev) + increment_factor
else
t_update = t_soisno(c,lev) + sign(clmT_max_increment, increment_factor)
Expand All @@ -1435,21 +1469,27 @@ subroutine update_clm_T(tstartcycle,mype)
if (ieee_is_nan(t_update)) then
print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev
else
t_soisno(c,lev) = t_update
t_soisno(c,lev) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update))
end if
end if
end do

col_updated(c) = .true.
end if ! col_updated

! --- TVEG: update with increment factor ---
cc = state_clm2pdaf_p(p, 2+nlevgrnd)
! Skip if no significant change in gridcell mean
if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then
if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7 .and. t_veg(p) < 1.0e20_r8) then
if( (clmincrement_type == 0)) then
increment_factor = clm_statevec(cc) / clm_statevec_orig(cc)
t_update = t_veg(p) * increment_factor
else
increment_factor = clm_statevec(cc) - clm_statevec_orig(cc)
if (abs(increment_factor) < clmT_max_increment) then
if (ieee_is_nan(increment_factor)) then
print *, "WARNING: t_veg increment_factor is NaN at p=", p, " - leaving t_veg unchanged"
t_update = t_veg(p)
else if (abs(increment_factor) < clmT_max_increment) then
t_update = t_veg(p) + increment_factor
else
t_update = t_veg(p) + sign(clmT_max_increment, increment_factor)
Expand All @@ -1459,7 +1499,7 @@ subroutine update_clm_T(tstartcycle,mype)
if (ieee_is_nan(t_update)) then
print *, "WARNING: t_veg update is NaN at p=", p
else
t_veg(p) = t_update
t_veg(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update))
end if
end if

Expand All @@ -1473,6 +1513,7 @@ subroutine update_clm_T(tstartcycle,mype)
incr_warn_count_soisno
if (incr_warn_count_veg > 0) print *, "WARNING: t_veg total increments exceeding T_max_increment:", &
incr_warn_count_veg
deallocate(col_updated)
endif

! Skin temperature updating skin, soil, vegetation and ground temperature.
Expand All @@ -1497,7 +1538,10 @@ subroutine update_clm_T(tstartcycle,mype)
t_update = t_skin(p) * increment_factor
else
increment_factor = clm_statevec(cc) - clm_statevec_orig(cc)
if (abs(increment_factor) < clmT_max_increment) then
if (ieee_is_nan(increment_factor)) then
print *, "WARNING: t_skin increment_factor is NaN at p=", p, " - leaving t_skin unchanged"
t_update = t_skin(p)
else if (abs(increment_factor) < clmT_max_increment) then
t_update = t_skin(p) + increment_factor
else
t_update = t_skin(p) + sign(clmT_max_increment, increment_factor)
Expand All @@ -1507,10 +1551,12 @@ subroutine update_clm_T(tstartcycle,mype)
if (ieee_is_nan(t_update)) then
print *, "WARNING: t_skin update is NaN at p=", p
else
t_skin(p) = t_update
t_skin(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update))
end if
end if

! Guard: only update each column once even if multiple patches share it.
if (.not. col_updated(c)) then
! --- TSOIL: update with increment factor for each layer ---
do lev=1,nlevgrnd
cc = state_clm2pdaf_p(p, 1+lev)
Expand All @@ -1521,7 +1567,10 @@ subroutine update_clm_T(tstartcycle,mype)
t_update = t_soisno(c,lev) * increment_factor
else
increment_factor = clm_statevec(cc) - clm_statevec_orig(cc)
if (abs(increment_factor) < clmT_max_increment) then
if (ieee_is_nan(increment_factor)) then
print *, "WARNING: t_soisno increment_factor is NaN at c=", c, " lev=", lev, " - leaving t_soisno unchanged"
t_update = t_soisno(c,lev)
else if (abs(increment_factor) < clmT_max_increment) then
t_update = t_soisno(c,lev) + increment_factor
else
t_update = t_soisno(c,lev) + sign(clmT_max_increment, increment_factor)
Expand All @@ -1531,21 +1580,25 @@ subroutine update_clm_T(tstartcycle,mype)
if (ieee_is_nan(t_update)) then
print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev
else
t_soisno(c,lev) = t_update
t_soisno(c,lev) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update))
end if
end if
end do
end if ! col_updated

! --- TVEG: update with increment factor ---
cc = state_clm2pdaf_p(p, 2+nlevgrnd)
! Skip if no significant change in gridcell mean
if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then
if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7 .and. t_veg(p) < 1.0e20_r8) then
if( (clmincrement_type == 0)) then
increment_factor = clm_statevec(cc) / clm_statevec_orig(cc)
t_update = t_veg(p) * increment_factor
else
increment_factor = clm_statevec(cc) - clm_statevec_orig(cc)
if (abs(increment_factor) < clmT_max_increment) then
if (ieee_is_nan(increment_factor)) then
print *, "WARNING: t_veg increment_factor is NaN at p=", p, " - leaving t_veg unchanged"
t_update = t_veg(p)
else if (abs(increment_factor) < clmT_max_increment) then
t_update = t_veg(p) + increment_factor
else
t_update = t_veg(p) + sign(clmT_max_increment, increment_factor)
Expand All @@ -1555,10 +1608,12 @@ subroutine update_clm_T(tstartcycle,mype)
if (ieee_is_nan(t_update)) then
print *, "WARNING: t_veg update is NaN at p=", p
else
t_veg(p) = t_update
t_veg(p) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update))
end if
end if

! Guard: only update each column once even if multiple patches share it.
if (.not. col_updated(c)) then
! --- TGRND: update with increment factor ---
cc = state_clm2pdaf_p(p, 3+nlevgrnd)
! Skip if no significant change in gridcell mean
Expand All @@ -1568,7 +1623,10 @@ subroutine update_clm_T(tstartcycle,mype)
t_update = t_grnd(c) * increment_factor
else
increment_factor = clm_statevec(cc) - clm_statevec_orig(cc)
if (abs(increment_factor) < clmT_max_increment) then
if (ieee_is_nan(increment_factor)) then
print *, "WARNING: t_grnd increment_factor is NaN at c=", c, " - leaving t_grnd unchanged"
t_update = t_grnd(c)
else if (abs(increment_factor) < clmT_max_increment) then
t_update = t_grnd(c) + increment_factor
else
t_update = t_grnd(c) + sign(clmT_max_increment, increment_factor)
Expand All @@ -1578,10 +1636,13 @@ subroutine update_clm_T(tstartcycle,mype)
if (ieee_is_nan(t_update)) then
print *, "WARNING: t_grnd update is NaN at c=", c
else
t_grnd(c) = t_update
t_grnd(c) = max(SHR_CONST_TKFRZ - 130.0_r8, min(SHR_CONST_TKFRZ + 100.0_r8, t_update))
end if
end if

col_updated(c) = .true.
end if ! col_updated

end if mask_freeze_2
end if mask_snow_2

Expand All @@ -1594,6 +1655,7 @@ subroutine update_clm_T(tstartcycle,mype)
incr_warn_count_veg
if (incr_warn_count_grnd > 0) print *, "WARNING: t_grnd total increments exceeding T_max_increment:", &
incr_warn_count_grnd
deallocate(col_updated)
endif

#ifdef PDAF_DEBUG
Expand Down