From eaddc9ba1bd4ae128d5ba6b3ecc17f798344b59d Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Aug 2026 16:37:41 -0600 Subject: [PATCH 01/20] Spread the DYNBAL water and energy fluxes over many years Rather than releasing all of the DYNBAL fluxes in the year in which they're generated, instead put them in storage pools that are released with a 20-year turnover time. This reduces these fluxes, especially in the oscillatory case that is somewhat common with a dynamic ice sheet. This partially addresses https://github.com/ESCOMP/CTSM/issues/4188 (specifically, option (1) in that issue). --- src/biogeophys/BalanceCheckMod.F90 | 50 +++++++-------------- src/biogeophys/EnergyFluxType.F90 | 12 ------ src/biogeophys/TemperatureType.F90 | 16 +++++++ src/biogeophys/WaterFluxType.F90 | 19 -------- src/biogeophys/WaterStateType.F90 | 47 ++++++++++++++++++++ src/dyn_subgrid/dynConsBiogeophysMod.F90 | 55 +++++++++++++++++------- src/dyn_subgrid/dynSubgridDriverMod.F90 | 5 +++ 7 files changed, 123 insertions(+), 81 deletions(-) diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index b79fcee46e..ff2ed56f4d 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -161,7 +161,6 @@ subroutine WaterGridcellBalance(bounds, & water_inst%bulk_and_tracers(i)%waterstate_inst, & water_inst%bulk_and_tracers(i)%waterdiagnostic_inst, & water_inst%bulk_and_tracers(i)%waterbalance_inst, & - water_inst%bulk_and_tracers(i)%waterflux_inst, & use_aquifer_layer = use_aquifer_layer, flag = flag) end do @@ -212,7 +211,7 @@ end subroutine BeginWaterColumnBalance subroutine WaterGridcellBalanceSingle(bounds, & num_nolakec, filter_nolakec, num_lakec, filter_lakec, & lakestate_inst, waterstate_inst, waterdiagnostic_inst, & - waterbalance_inst, waterflux_inst, use_aquifer_layer, flag) + waterbalance_inst, use_aquifer_layer, flag) ! ! !DESCRIPTION: ! Grid cell-level water balance for bulk or a single tracer @@ -232,7 +231,6 @@ subroutine WaterGridcellBalanceSingle(bounds, & class(waterstate_type) , intent(inout) :: waterstate_inst class(waterdiagnostic_type), intent(in) :: waterdiagnostic_inst class(waterbalance_type) , intent(inout) :: waterbalance_inst - class(waterflux_type) , intent(inout) :: waterflux_inst logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run character(len=5) , intent(in) :: flag ! specifies begwb or endwb ! @@ -241,8 +239,6 @@ subroutine WaterGridcellBalanceSingle(bounds, & integer :: begc, endc, begl, endl, begg, endg ! bounds real(r8) :: wb_col(bounds%begc:bounds%endc) ! temporary column-level water mass real(r8) :: wb_grc(bounds%begg:bounds%endg) ! temporary grid cell-level water mass - real(r8) :: qflx_liq_dynbal_left_to_dribble(bounds%begg:bounds%endg) ! grc liq dynamic land cover change conversion runoff flux - real(r8) :: qflx_ice_dynbal_left_to_dribble(bounds%begg:bounds%endg) ! grc ice dynamic land cover change conversion runoff flux real(r8) :: wa_reset_nonconservation_gain_grc(bounds%begg:bounds%endg) ! grc mass gained from resetting water in the unconfined aquifer, wa_col (negative indicates mass lost) (mm) character(len=*), parameter :: subname = 'WaterGridcellBalanceSingle' @@ -283,38 +279,18 @@ subroutine WaterGridcellBalanceSingle(bounds, & enddo endif - ! Call the beginning or ending version of the subroutine according - ! to flag value - if (flag == 'begwb') then - call waterflux_inst%qflx_liq_dynbal_dribbler%get_amount_left_to_dribble_beg( & - bounds, & - qflx_liq_dynbal_left_to_dribble(begg:endg)) - call waterflux_inst%qflx_ice_dynbal_dribbler%get_amount_left_to_dribble_beg( & - bounds, & - qflx_ice_dynbal_left_to_dribble(begg:endg)) - else if (flag == 'endwb') then - call waterflux_inst%qflx_liq_dynbal_dribbler%get_amount_left_to_dribble_end( & - bounds, & - qflx_liq_dynbal_left_to_dribble(begg:endg)) - call waterflux_inst%qflx_ice_dynbal_dribbler%get_amount_left_to_dribble_end( & - bounds, & - qflx_ice_dynbal_left_to_dribble(begg:endg)) - else - write(iulog,*) 'Unknown flag passed into this subroutine.' - write(iulog,*) 'Expecting either begwb or endwb.' - call endrun(msg=errmsg(sourcefile, __LINE__)) - end if - - ! These dynbal dribblers store the delta state, (end - beg). Thus, the - ! amount dribbled out is the negative of the amount stored in the - ! dribblers. Therefore, conservation requires us to subtract the amount - ! remaining to dribble. + ! Account for the dynbal storage pools. These pools hold the delta state, + ! (end - beg), arising from dynamic landunit adjustments; this amount is released + ! gradually to the dynbal fluxes. Thus, the amount released is the negative of the + ! amount held in the storage pools, so conservation requires us to subtract the + ! storage. ! This sign convention is opposite to the convention chosen for the - ! respective dribble terms used in the carbon balance. At some point - ! it may be worth making the two conventions consistent. + ! dribble terms used in the carbon balance. At some point it may be worth making the + ! two conventions consistent. do g = begg, endg - wb_grc(g) = wb_grc(g) - qflx_liq_dynbal_left_to_dribble(g) & - - qflx_ice_dynbal_left_to_dribble(g) + wb_grc(g) = wb_grc(g) & + - waterstate_inst%dynbal_liq_storage_grc(g) & + - waterstate_inst%dynbal_ice_storage_grc(g) end do ! Map wb_grc to beginning/ending water balance according to flag @@ -343,6 +319,10 @@ subroutine WaterGridcellBalanceSingle(bounds, & do g = begg, endg endwb_grc(g) = wb_grc(g) - wa_reset_nonconservation_gain_grc(g) end do + else + write(iulog,*) 'Unknown flag passed into this subroutine.' + write(iulog,*) 'Expecting either begwb or endwb.' + call endrun(msg=errmsg(sourcefile, __LINE__)) end if end associate diff --git a/src/biogeophys/EnergyFluxType.F90 b/src/biogeophys/EnergyFluxType.F90 index df8abb3f01..f05769e05e 100644 --- a/src/biogeophys/EnergyFluxType.F90 +++ b/src/biogeophys/EnergyFluxType.F90 @@ -13,7 +13,6 @@ module EnergyFluxType use LandunitType , only : lun use ColumnType , only : col use PatchType , only : patch - use AnnualFluxDribbler, only : annual_flux_dribbler_type, annual_flux_dribbler_gridcell ! implicit none save @@ -117,10 +116,6 @@ module EnergyFluxType real(r8), pointer :: errlon_patch (:) ! longwave radiation conservation error (W/m**2) real(r8), pointer :: errlon_col (:) ! longwave radiation conservation error (W/m**2) - ! Objects that help convert once-per-year dynamic land cover changes into fluxes - ! that are dribbled throughout the year - type(annual_flux_dribbler_type) :: eflx_dynbal_dribbler - contains procedure, public :: Init ! Public initialization method @@ -272,11 +267,6 @@ subroutine InitAllocate(this, bounds) allocate( this%errlon_patch (begp:endp)) ; this%errlon_patch (:) = nan allocate( this%errlon_col (begc:endc)) ; this%errlon_col (:) = nan - this%eflx_dynbal_dribbler = annual_flux_dribbler_gridcell( & - bounds = bounds, & - name = 'eflx_dynbal', & - units = 'J/m**2') - end subroutine InitAllocate !------------------------------------------------------------------------ @@ -903,8 +893,6 @@ subroutine Restart(this, bounds, ncid, flag, is_simple_buildtemp, is_prog_buildt long_name='instantaneous daily minimum of transpiration wetness factor', units='', & interpinic_flag='interp', readvar=readvar, data=this%btran_min_inst_patch) - call this%eflx_dynbal_dribbler%Restart(bounds, ncid, flag) - end subroutine Restart !----------------------------------------------------------------------- subroutine InitAccBuffer (this, bounds) diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 6fd8faf037..79e90de6e5 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -106,6 +106,7 @@ module TemperatureType real(r8), pointer :: heat2_grc (:) ! grc post land cover change total heat content real(r8), pointer :: liquid_water_temp1_grc (:) ! grc initial weighted average liquid water temperature (K) real(r8), pointer :: liquid_water_temp2_grc (:) ! grc post land cover change weighted average liquid water temperature (K) + real(r8), pointer :: dynbal_heat_storage_grc (:) ! grc heat storage from dynbal adjustments, to be released gradually [J/m^2] ! Flags integer , pointer :: imelt_col (:,:) ! flag for melting (=1), freezing (=2), Not=0 (-nlevsno+1:nlevgrnd) @@ -288,6 +289,7 @@ subroutine InitAllocate(this, bounds) allocate(this%heat2_grc (begg:endg)) ; this%heat2_grc (:) = nan allocate(this%liquid_water_temp1_grc (begg:endg)) ; this%liquid_water_temp1_grc (:) = nan allocate(this%liquid_water_temp2_grc (begg:endg)) ; this%liquid_water_temp2_grc (:) = nan + allocate(this%dynbal_heat_storage_grc (begg:endg)) ; this%dynbal_heat_storage_grc (:) = nan ! flags allocate(this%imelt_col (begc:endc,-nlevsno+1:nlevmaxurbgrnd)) ; this%imelt_col (:,:) = huge(1) @@ -555,6 +557,11 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp ) avgflag='A', long_name='initial gridcell weighted average liquid water temperature', & ptr_lnd=this%liquid_water_temp1_grc, default='inactive') + this%dynbal_heat_storage_grc(begg:endg) = spval + call hist_addfld1d (fname='DYNBAL_HEAT_STORAGE', units='J/m^2', & + avgflag='A', long_name='gridcell heat storage from dynamic landunit adjustments', & + ptr_lnd=this%dynbal_heat_storage_grc, default='inactive') + this%snot_top_col(begc:endc) = spval call hist_addfld1d (fname='SNOTTOPL', units='K', & avgflag='A', long_name='snow temperature (top layer)', & @@ -919,6 +926,9 @@ subroutine InitCold(this, bounds, & ! set. this%dynbal_baseline_heat_col(bounds%begc:bounds%endc) = 0._r8 + ! The dynbal storage pools start out empty + this%dynbal_heat_storage_grc(bounds%begg:bounds%endg) = 0._r8 + end subroutine InitCold !------------------------------------------------------------------------ @@ -1075,6 +1085,12 @@ subroutine Restart(this, bounds, ncid, flag, is_simple_buildtemp, is_prog_buildt units='J/m2', & interpinic_flag='interp', readvar=readvar, data=this%dynbal_baseline_heat_col) + call restartvar(ncid=ncid, flag=flag, varname='DYNBAL_HEAT_STORAGE', xtype=ncd_double, & + dim1name='gridcell', & + long_name="heat storage from dynbal adjustments", & + units='J/m2', & + interpinic_flag='interp', readvar=readvar, data=this%dynbal_heat_storage_grc) + if (use_crop) then call restartvar(ncid=ncid, flag=flag, varname='gdd1020', xtype=ncd_double, & dim1name='pft', long_name='20 year average of growing degree-days base 10C from planting', units='ddays', & diff --git a/src/biogeophys/WaterFluxType.F90 b/src/biogeophys/WaterFluxType.F90 index 88b1b3e5cf..9dd058ef84 100644 --- a/src/biogeophys/WaterFluxType.F90 +++ b/src/biogeophys/WaterFluxType.F90 @@ -13,7 +13,6 @@ module WaterFluxType use decompMod , only : subgrid_level_patch, subgrid_level_column, subgrid_level_landunit, subgrid_level_gridcell use LandunitType , only : lun use ColumnType , only : col - use AnnualFluxDribbler, only : annual_flux_dribbler_type, annual_flux_dribbler_gridcell use WaterInfoBaseType, only : water_info_base_type use WaterTracerContainerType, only : water_tracer_container_type use WaterTracerUtils, only : AllocateVar1d, AllocateVar2d @@ -108,11 +107,6 @@ module WaterFluxType real(r8), pointer :: qflx_irrig_drip_patch (:) ! patch drip irrigation real(r8), pointer :: qflx_irrig_sprinkler_patch(:) ! patch sprinkler irrigation - ! Objects that help convert once-per-year dynamic land cover changes into fluxes - ! that are dribbled throughout the year - type(annual_flux_dribbler_type) :: qflx_liq_dynbal_dribbler - type(annual_flux_dribbler_type) :: qflx_ice_dynbal_dribbler - contains procedure, public :: Init @@ -384,16 +378,6 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%qflx_irrig_sprinkler_patch, name = 'qflx_irrig_sprinkler_patch', & container = tracer_vars, & bounds = bounds, subgrid_level = subgrid_level_patch) - - this%qflx_liq_dynbal_dribbler = annual_flux_dribbler_gridcell( & - bounds = bounds, & - name = this%info%fname('qflx_liq_dynbal'), & - units = 'mm H2O') - - this%qflx_ice_dynbal_dribbler = annual_flux_dribbler_gridcell( & - bounds = bounds, & - name = this%info%fname('qflx_ice_dynbal'), & - units = 'mm H2O') end subroutine InitAllocate @@ -959,9 +943,6 @@ subroutine Restart(this, bounds, ncid, flag) this%qflx_snow_drain_col(bounds%begc:bounds%endc) = 0._r8 endif - call this%qflx_liq_dynbal_dribbler%Restart(bounds, ncid, flag) - call this%qflx_ice_dynbal_dribbler%Restart(bounds, ncid, flag) - end subroutine Restart end module WaterFluxType diff --git a/src/biogeophys/WaterStateType.F90 b/src/biogeophys/WaterStateType.F90 index 35441d65d9..ffe32e9d13 100644 --- a/src/biogeophys/WaterStateType.F90 +++ b/src/biogeophys/WaterStateType.F90 @@ -51,6 +51,11 @@ module WaterStateType real(r8), pointer :: dynbal_baseline_liq_col(:) ! baseline liquid water content subtracted from each column's total liquid water calculation (mm H2O) real(r8), pointer :: dynbal_baseline_ice_col(:) ! baseline ice content subtracted from each column's total ice calculation (mm H2O) + ! Storage pools holding the water added / removed by dynamic landunit adjustments; + ! these are released gradually to the dynbal fluxes + real(r8), pointer :: dynbal_liq_storage_grc(:) ! grc liquid water storage from dynbal adjustments, to be released gradually (mm H2O) + real(r8), pointer :: dynbal_ice_storage_grc(:) ! grc ice storage from dynbal adjustments, to be released gradually (mm H2O) + real(r8) :: aquifer_water_baseline ! baseline value for water in the unconfined aquifer (wa_col) for this bulk / tracer (mm) real(r8), pointer :: excess_ice_col (:,:) ! col excess ice (kg/m2) (new) (-nlevsno+1:nlevgrnd) @@ -160,6 +165,12 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%dynbal_baseline_ice_col, name = 'dynbal_baseline_ice_col', & container = tracer_vars, & bounds = bounds, subgrid_level = subgrid_level_column) + call AllocateVar1d(var = this%dynbal_liq_storage_grc, name = 'dynbal_liq_storage_grc', & + container = tracer_vars, & + bounds = bounds, subgrid_level = subgrid_level_gridcell) + call AllocateVar1d(var = this%dynbal_ice_storage_grc, name = 'dynbal_ice_storage_grc', & + container = tracer_vars, & + bounds = bounds, subgrid_level = subgrid_level_gridcell) call AllocateVar1d(var = this%stream_water_volume_lun, name = 'stream_water_volume_lun', & container = tracer_vars, & bounds = bounds, subgrid_level = subgrid_level_landunit) @@ -292,6 +303,22 @@ subroutine InitHistory(this, bounds, use_aquifer_layer) ptr_col=this%wa_col, l2g_scale_type='veg') end if + this%dynbal_liq_storage_grc(begg:endg) = spval + call hist_addfld1d ( & + fname=this%info%fname('DYNBAL_LIQ_STORAGE'), & + units='mm', & + avgflag='A', & + long_name=this%info%lname('gridcell liquid water storage from dynamic landunit adjustments'), & + ptr_lnd=this%dynbal_liq_storage_grc, default='inactive') + + this%dynbal_ice_storage_grc(begg:endg) = spval + call hist_addfld1d ( & + fname=this%info%fname('DYNBAL_ICE_STORAGE'), & + units='mm', & + avgflag='A', & + long_name=this%info%lname('gridcell ice storage from dynamic landunit adjustments'), & + ptr_lnd=this%dynbal_ice_storage_grc, default='inactive') + if (use_hillslope) then this%stream_water_volume_lun(begl:endl) = spval call hist_addfld1d (fname=this%info%fname('STREAM_WATER_VOLUME'), units='m3', & @@ -549,6 +576,10 @@ subroutine InitCold(this, bounds, & this%dynbal_baseline_liq_col(bounds%begc:bounds%endc) = 0._r8 this%dynbal_baseline_ice_col(bounds%begc:bounds%endc) = 0._r8 + ! The dynbal storage pools start out empty + this%dynbal_liq_storage_grc(bounds%begg:bounds%endg) = 0._r8 + this%dynbal_ice_storage_grc(bounds%begg:bounds%endg) = 0._r8 + !Initialize excess ice this%exice_bulk_init(bounds%begc:bounds%endc) = exice_init_conc_col(bounds%begc:bounds%endc) this%excess_ice_col(bounds%begc:bounds%endc,:) = 0.0_r8 @@ -714,6 +745,22 @@ subroutine Restart(this, bounds, ncid, flag, & units='kg/m2', & interpinic_flag='interp', readvar=readvar, data=this%dynbal_baseline_ice_col) + call restartvar(ncid=ncid, flag=flag, & + varname=this%info%fname('DYNBAL_LIQ_STORAGE'), & + xtype=ncd_double, & + dim1name='gridcell', & + long_name=this%info%lname('liquid water storage from dynbal adjustments'), & + units='kg/m2', & + interpinic_flag='interp', readvar=readvar, data=this%dynbal_liq_storage_grc) + + call restartvar(ncid=ncid, flag=flag, & + varname=this%info%fname('DYNBAL_ICE_STORAGE'), & + xtype=ncd_double, & + dim1name='gridcell', & + long_name=this%info%lname('ice storage from dynbal adjustments'), & + units='kg/m2', & + interpinic_flag='interp', readvar=readvar, data=this%dynbal_ice_storage_grc) + call restartvar(ncid=ncid, flag=flag, & varname=this%info%fname('STREAM_WATER_VOLUME'), & xtype=ncd_double, & diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 9ef474d401..672b66c311 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -11,6 +11,7 @@ module dynConsBiogeophysMod ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_const_mod , only : SHR_CONST_CDAY use decompMod , only : bounds_type use UrbanParamsType , only : urbanparams_type use EnergyFluxType , only : energyflux_type @@ -40,6 +41,7 @@ module dynConsBiogeophysMod use landunit_varcon , only : istsoil, istice use dynSubgridControlMod , only : get_for_testing_zero_dynbal_fluxes use filterColMod , only : filter_col_type, col_filter_from_ltypes + use clm_time_manager , only : get_step_size_real ! ! !PUBLIC MEMBER FUNCTIONS: implicit none @@ -57,6 +59,16 @@ module dynConsBiogeophysMod private :: dyn_water_content ! compute gridcell total liquid and ice water contents private :: dyn_heat_content ! compute gridcell total heat contents + ! + ! !PRIVATE DATA MEMBERS: + + ! Turnover rate of the dynbal storage pools [1/s]. Note that we specify the time in + ! years and then convert it to seconds. (It isn't important to account for leap years + ! here: we just want a value that roughly equates to our desired residence time in + ! years: it doesn't need to be exact.) + real(r8), parameter, private :: dynbal_storage_turnover_rate = & + 1._r8 / (20._r8 * (365._r8 * SHR_CONST_CDAY)) + character(len=*), parameter, private :: sourcefile = & __FILE__ !--------------------------------------------------------------------------- @@ -506,6 +518,7 @@ subroutine dyn_hwcontent_final(bounds, & ! !LOCAL VARIABLES: integer :: i integer :: g ! grid cell index + real(r8) :: dtime ! model time step [s] real(r8) :: this_delta_liq(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk or one tracer real(r8) :: delta_liq_bulk(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk water real(r8) :: delta_heat(bounds%begg:bounds%endg) ! change in gridcell heat content @@ -513,7 +526,10 @@ subroutine dyn_hwcontent_final(bounds, & associate( & begg => bounds%begg, & - endg => bounds%endg) + endg => bounds%endg, & + dynbal_heat_storage => temperature_inst%dynbal_heat_storage_grc, & ! Output: [real(r8) (:)] heat storage from dynbal adjustments, to be released gradually (J/m^2) + eflx_dynbal => energyflux_inst%eflx_dynbal_grc & ! Output: [real(r8) (:)] dynamic land cover change conversion energy flux (W/m^2) + ) do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end associate(bulk_or_tracer => water_inst%bulk_and_tracers(i)) @@ -559,10 +575,12 @@ subroutine dyn_hwcontent_final(bounds, & liquid_water_temp2 = temperature_inst%liquid_water_temp2_grc(begg:endg), & delta_heat = delta_heat(begg:endg)) - call energyflux_inst%eflx_dynbal_dribbler%set_curr_delta(bounds, & - delta_heat(begg:endg)) - call energyflux_inst%eflx_dynbal_dribbler%get_curr_flux(bounds, & - energyflux_inst%eflx_dynbal_grc(begg:endg)) + dtime = get_step_size_real() + do g = begg, endg + dynbal_heat_storage(g) = dynbal_heat_storage(g) + delta_heat(g) + eflx_dynbal(g) = dynbal_heat_storage(g) * dynbal_storage_turnover_rate + dynbal_heat_storage(g) = dynbal_heat_storage(g) - (eflx_dynbal(g) * dtime) + end do end associate @@ -586,7 +604,7 @@ subroutine dyn_water_content_final(bounds, & integer , intent(in) :: filter_nolakec(:) integer , intent(in) :: num_lakec integer , intent(in) :: filter_lakec(:) - class(waterstate_type) , intent(in) :: waterstate_inst + class(waterstate_type) , intent(inout) :: waterstate_inst class(waterdiagnostic_type) , intent(in) :: waterdiagnostic_inst class(waterbalance_type) , intent(inout) :: waterbalance_inst class(waterflux_type) , intent(inout) :: waterflux_inst @@ -595,6 +613,7 @@ subroutine dyn_water_content_final(bounds, & ! ! !LOCAL VARIABLES: integer :: g + real(r8) :: dtime ! model time step [s] real(r8) :: delta_ice(bounds%begg:bounds%endg) ! change in gridcell h2o ice content character(len=*), parameter :: subname = 'dyn_water_content_final' @@ -604,7 +623,12 @@ subroutine dyn_water_content_final(bounds, & associate( & begg => bounds%begg, & - endg => bounds%endg) + endg => bounds%endg, & + dynbal_liq_storage => waterstate_inst%dynbal_liq_storage_grc, & ! Output: [real(r8) (:)] liquid water storage from dynbal adjustments, to be released gradually (mm H2O) + dynbal_ice_storage => waterstate_inst%dynbal_ice_storage_grc, & ! Output: [real(r8) (:)] ice storage from dynbal adjustments, to be released gradually (mm H2O) + qflx_liq_dynbal => waterflux_inst%qflx_liq_dynbal_grc, & ! Output: [real(r8) (:)] liq dynamic land cover change conversion runoff flux (mm H2O/s) + qflx_ice_dynbal => waterflux_inst%qflx_ice_dynbal_grc & ! Output: [real(r8) (:)] ice dynamic land cover change conversion runoff flux (mm H2O/s) + ) call dyn_water_content(bounds, & num_nolakec, filter_nolakec, & @@ -625,15 +649,16 @@ subroutine dyn_water_content_final(bounds, & end do end if - call waterflux_inst%qflx_liq_dynbal_dribbler%set_curr_delta(bounds, & - delta_liq(begg:endg)) - call waterflux_inst%qflx_liq_dynbal_dribbler%get_curr_flux(bounds, & - waterflux_inst%qflx_liq_dynbal_grc(begg:endg)) + dtime = get_step_size_real() + do g = begg, endg + dynbal_liq_storage(g) = dynbal_liq_storage(g) + delta_liq(g) + qflx_liq_dynbal(g) = dynbal_liq_storage(g) * dynbal_storage_turnover_rate + dynbal_liq_storage(g) = dynbal_liq_storage(g) - (qflx_liq_dynbal(g) * dtime) - call waterflux_inst%qflx_ice_dynbal_dribbler%set_curr_delta(bounds, & - delta_ice(begg:endg)) - call waterflux_inst%qflx_ice_dynbal_dribbler%get_curr_flux(bounds, & - waterflux_inst%qflx_ice_dynbal_grc(begg:endg)) + dynbal_ice_storage(g) = dynbal_ice_storage(g) + delta_ice(g) + qflx_ice_dynbal(g) = dynbal_ice_storage(g) * dynbal_storage_turnover_rate + dynbal_ice_storage(g) = dynbal_ice_storage(g) - (qflx_ice_dynbal(g) * dtime) + end do end associate diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index ea1210521d..8110a47974 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -328,6 +328,11 @@ subroutine dynSubgrid_driver(bounds_proc, ! is at least *mostly* true, but there may be some subtleties, like resetting of ! some variables, that are needed even in (some) time steps where we haven't ! changed weights. + ! + ! NOTE(wjs, 2026-08-28) Actually, the above statement is NOT true for the call to + ! dyn_hwcontent_final: that needs to be called every time step to properly do the + ! outflux from the dynbal storage pools. I haven't gone back and checked other + ! pieces to see what else may need to be done every time step. ! ======================================================================== call dynSubgrid_wrapup_weight_changes(bounds_clump, glc_behavior) From 9681949c28467c7241ec8c7d6af0fd2da98a9c61 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Aug 2026 17:36:14 -0600 Subject: [PATCH 02/20] Change sign convention for water and energy dynbal fluxes and storage I found the previous sign convention unintuitive. While I'm changing these fluxes anyway, I'm going to flip the sign convention to be more intuitive. --- src/biogeophys/BalanceCheckMod.F90 | 16 ++++++---------- src/biogeophys/EnergyFluxType.F90 | 4 ++-- src/biogeophys/TemperatureType.F90 | 2 +- src/biogeophys/WaterFluxType.F90 | 10 ++++++---- src/biogeophys/WaterStateType.F90 | 4 +++- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 12 +++++++++--- src/main/lnd2atmMod.F90 | 8 ++++---- 7 files changed, 31 insertions(+), 25 deletions(-) diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index ff2ed56f4d..b6e32c2010 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -279,18 +279,14 @@ subroutine WaterGridcellBalanceSingle(bounds, & enddo endif - ! Account for the dynbal storage pools. These pools hold the delta state, - ! (end - beg), arising from dynamic landunit adjustments; this amount is released - ! gradually to the dynbal fluxes. Thus, the amount released is the negative of the - ! amount held in the storage pools, so conservation requires us to subtract the - ! storage. - ! This sign convention is opposite to the convention chosen for the - ! dribble terms used in the carbon balance. At some point it may be worth making the - ! two conventions consistent. + ! Add the dynbal storage pools. These pools hold the amount of water that must be + ! added to the gridcell's total water content in order for that total to be conserved + ! across dynamic landunit adjustments; this water is released gradually to the dynbal + ! fluxes. do g = begg, endg wb_grc(g) = wb_grc(g) & - - waterstate_inst%dynbal_liq_storage_grc(g) & - - waterstate_inst%dynbal_ice_storage_grc(g) + + waterstate_inst%dynbal_liq_storage_grc(g) & + + waterstate_inst%dynbal_ice_storage_grc(g) end do ! Map wb_grc to beginning/ending water balance according to flag diff --git a/src/biogeophys/EnergyFluxType.F90 b/src/biogeophys/EnergyFluxType.F90 index f05769e05e..4bb95cf177 100644 --- a/src/biogeophys/EnergyFluxType.F90 +++ b/src/biogeophys/EnergyFluxType.F90 @@ -52,7 +52,7 @@ module EnergyFluxType real(r8), pointer :: eflx_snomelt_u_col (:) ! col urban snow melt heat flux (W/m**2) real(r8), pointer :: eflx_gnet_patch (:) ! patch net heat flux into ground (W/m**2) real(r8), pointer :: eflx_grnd_lake_patch (:) ! patch net heat flux into lake / snow surface, excluding light transmission (W/m**2) - real(r8), pointer :: eflx_dynbal_grc (:) ! grc dynamic land cover change conversion energy flux (W/m**2) + real(r8), pointer :: eflx_dynbal_grc (:) ! grc dynamic land cover change conversion energy flux (W/m**2) [+ to atm] real(r8), pointer :: eflx_bot_col (:) ! col heat flux from beneath the soil or ice column (W/m**2) real(r8), pointer :: eflx_fgr12_col (:) ! col ground heat flux between soil layers 1 and 2 (W/m**2) real(r8), pointer :: eflx_fgr_col (:,:) ! col (rural) soil downward heat flux (W/m2) (1:nlevgrnd) (pos upward; usually eflx_bot >= 0) @@ -309,7 +309,7 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp) this%eflx_dynbal_grc(begg:endg) = spval call hist_addfld1d (fname='EFLX_DYNBAL', units='W/m^2', & - avgflag='A', long_name='dynamic land cover change conversion energy flux', & + avgflag='A', long_name='dynamic land cover change conversion energy flux [+ to atm]', & ptr_lnd=this%eflx_dynbal_grc) this%eflx_snomelt_col(begc:endc) = spval diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 79e90de6e5..652a78f9e9 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -106,7 +106,7 @@ module TemperatureType real(r8), pointer :: heat2_grc (:) ! grc post land cover change total heat content real(r8), pointer :: liquid_water_temp1_grc (:) ! grc initial weighted average liquid water temperature (K) real(r8), pointer :: liquid_water_temp2_grc (:) ! grc post land cover change weighted average liquid water temperature (K) - real(r8), pointer :: dynbal_heat_storage_grc (:) ! grc heat storage from dynbal adjustments, to be released gradually [J/m^2] + real(r8), pointer :: dynbal_heat_storage_grc (:) ! grc heat storage from dynbal adjustments, to be released gradually; positive means heat that still needs to be released to the atmosphere [J/m^2] ! Flags integer , pointer :: imelt_col (:,:) ! flag for melting (=1), freezing (=2), Not=0 (-nlevsno+1:nlevgrnd) diff --git a/src/biogeophys/WaterFluxType.F90 b/src/biogeophys/WaterFluxType.F90 index 9dd058ef84..561a658c21 100644 --- a/src/biogeophys/WaterFluxType.F90 +++ b/src/biogeophys/WaterFluxType.F90 @@ -97,8 +97,8 @@ module WaterFluxType real(r8), pointer :: qflx_snow_percolation_col(:,:) ! col liquid percolation out of the bottom of snow layer j (mm H2O /s) ! Dynamic land cover change - real(r8), pointer :: qflx_liq_dynbal_grc (:) ! grc liq dynamic land cover change conversion runoff flux - real(r8), pointer :: qflx_ice_dynbal_grc (:) ! grc ice dynamic land cover change conversion runoff flux + real(r8), pointer :: qflx_liq_dynbal_grc (:) ! grc liq dynamic land cover change conversion runoff flux (mm H2O/s) (positive means addition to runoff) + real(r8), pointer :: qflx_ice_dynbal_grc (:) ! grc ice dynamic land cover change conversion runoff flux (mm H2O/s) (positive means addition to runoff) real(r8), pointer :: qflx_sfc_irrig_col (:) ! col surface irrigation flux (mm H2O/s) [+] real(r8), pointer :: qflx_gw_uncon_irrig_col (:) ! col unconfined groundwater irrigation flux (mm H2O/s) @@ -531,7 +531,8 @@ subroutine InitHistory(this, bounds) fname=this%info%fname('QFLX_LIQ_DYNBAL'), & units='mm/s', & avgflag='A', & - long_name=this%info%lname('liq dynamic land cover change conversion runoff flux'), & + long_name=this%info%lname( & + 'liq dynamic land cover change conversion runoff flux (positive means addition to runoff)'), & ptr_lnd=this%qflx_liq_dynbal_grc) this%qflx_ice_dynbal_grc(begg:endg) = spval @@ -539,7 +540,8 @@ subroutine InitHistory(this, bounds) fname=this%info%fname('QFLX_ICE_DYNBAL'), & units='mm/s', & avgflag='A', & - long_name=this%info%lname('ice dynamic land cover change conversion runoff flux'), & + long_name=this%info%lname( & + 'ice dynamic land cover change conversion runoff flux (positive means addition to runoff)'), & ptr_lnd=this%qflx_ice_dynbal_grc) this%qflx_runoff_col(begc:endc) = spval diff --git a/src/biogeophys/WaterStateType.F90 b/src/biogeophys/WaterStateType.F90 index ffe32e9d13..db43e30929 100644 --- a/src/biogeophys/WaterStateType.F90 +++ b/src/biogeophys/WaterStateType.F90 @@ -52,7 +52,9 @@ module WaterStateType real(r8), pointer :: dynbal_baseline_ice_col(:) ! baseline ice content subtracted from each column's total ice calculation (mm H2O) ! Storage pools holding the water added / removed by dynamic landunit adjustments; - ! these are released gradually to the dynbal fluxes + ! these are released gradually to the dynbal fluxes. For these variables, positive + ! values indicate mass that still needs to be released to runoff; negative values + ! indicate mass that still needs to be drawn back in from elsewhere. real(r8), pointer :: dynbal_liq_storage_grc(:) ! grc liquid water storage from dynbal adjustments, to be released gradually (mm H2O) real(r8), pointer :: dynbal_ice_storage_grc(:) ! grc ice storage from dynbal adjustments, to be released gradually (mm H2O) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 672b66c311..391c733e3a 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -575,9 +575,12 @@ subroutine dyn_hwcontent_final(bounds, & liquid_water_temp2 = temperature_inst%liquid_water_temp2_grc(begg:endg), & delta_heat = delta_heat(begg:endg)) + ! Add this time step's spurious heat gain to the storage pool (with a negative to use + ! the proper sign convention for this pool), then release a portion of the pool as + ! this time step's flux. dtime = get_step_size_real() do g = begg, endg - dynbal_heat_storage(g) = dynbal_heat_storage(g) + delta_heat(g) + dynbal_heat_storage(g) = dynbal_heat_storage(g) - delta_heat(g) eflx_dynbal(g) = dynbal_heat_storage(g) * dynbal_storage_turnover_rate dynbal_heat_storage(g) = dynbal_heat_storage(g) - (eflx_dynbal(g) * dtime) end do @@ -649,13 +652,16 @@ subroutine dyn_water_content_final(bounds, & end do end if + ! Add this time step's spurious water gain to the storage pools (with a negative to + ! use the proper sign convention for these pools), then release a portion of each pool + ! as this time step's flux. dtime = get_step_size_real() do g = begg, endg - dynbal_liq_storage(g) = dynbal_liq_storage(g) + delta_liq(g) + dynbal_liq_storage(g) = dynbal_liq_storage(g) - delta_liq(g) qflx_liq_dynbal(g) = dynbal_liq_storage(g) * dynbal_storage_turnover_rate dynbal_liq_storage(g) = dynbal_liq_storage(g) - (qflx_liq_dynbal(g) * dtime) - dynbal_ice_storage(g) = dynbal_ice_storage(g) + delta_ice(g) + dynbal_ice_storage(g) = dynbal_ice_storage(g) - delta_ice(g) qflx_ice_dynbal(g) = dynbal_ice_storage(g) * dynbal_storage_turnover_rate dynbal_ice_storage(g) = dynbal_ice_storage(g) - (qflx_ice_dynbal(g) * dtime) end do diff --git a/src/main/lnd2atmMod.F90 b/src/main/lnd2atmMod.F90 index 503f4b9585..765ee8fde8 100644 --- a/src/main/lnd2atmMod.F90 +++ b/src/main/lnd2atmMod.F90 @@ -273,7 +273,7 @@ subroutine lnd2atm(bounds, & do g = bounds%begg, bounds%endg lnd2atm_inst%eflx_sh_tot_grc(g) = lnd2atm_inst%eflx_sh_tot_grc(g) + & lnd2atm_inst%eflx_sh_precip_conversion_grc(g) + & - eflx_sh_ice_to_liq_grc(g) - & + eflx_sh_ice_to_liq_grc(g) + & energyflux_inst%eflx_dynbal_grc(g) enddo @@ -445,10 +445,10 @@ subroutine lnd2atm(bounds, & do g = bounds%begg, bounds%endg water_inst%waterlnd2atmbulk_inst%qflx_rofliq_qgwl_grc(g) = & - water_inst%waterlnd2atmbulk_inst%qflx_rofliq_qgwl_grc(g) - & + water_inst%waterlnd2atmbulk_inst%qflx_rofliq_qgwl_grc(g) + & water_inst%waterfluxbulk_inst%qflx_liq_dynbal_grc(g) water_inst%waterlnd2atmbulk_inst%qflx_rofliq_grc(g) = & - water_inst%waterlnd2atmbulk_inst%qflx_rofliq_grc(g) - & + water_inst%waterlnd2atmbulk_inst%qflx_rofliq_grc(g) + & water_inst%waterfluxbulk_inst%qflx_liq_dynbal_grc(g) enddo @@ -463,7 +463,7 @@ subroutine lnd2atm(bounds, & c2l_scale_type= 'urbanf', l2g_scale_type='unity' ) do g = bounds%begg, bounds%endg water_inst%waterlnd2atmbulk_inst%qflx_rofice_grc(g) = & - water_inst%waterlnd2atmbulk_inst%qflx_rofice_grc(g) - & + water_inst%waterlnd2atmbulk_inst%qflx_rofice_grc(g) + & water_inst%waterfluxbulk_inst%qflx_ice_dynbal_grc(g) enddo From 975d4100d6ce056188684835cacb3435cefb6fbb Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Aug 2026 18:24:28 -0600 Subject: [PATCH 03/20] Make dynbal storage residence time namelist controllable --- bld/CLMBuildNamelist.pm | 6 +++ bld/namelist_files/namelist_defaults_ctsm.xml | 1 + .../namelist_definition_ctsm.xml | 26 ++++++++++++ src/dyn_subgrid/dynConsBiogeophysMod.F90 | 19 ++++----- src/dyn_subgrid/dynSubgridControlMod.F90 | 42 ++++++++++++++++++- 5 files changed, 82 insertions(+), 12 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 85c9f472c5..560e1bdca7 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2827,6 +2827,12 @@ sub setup_logic_dynamic_subgrid { &remove_leading_and_trailing_quotes($nl_flags->{'clm_start_type'}) eq "branch") { $log->fatal_error("reset_dynbal_baselines has no effect in a branch run"); } + + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'dynbal_storage_residence_time'); + my $dynbal_storage_residence_time = $nl->get_value('dynbal_storage_residence_time'); + if ( $dynbal_storage_residence_time <= 0.0 ) { + $log->fatal_error("dynbal_storage_residence_time must be greater than 0"); + } } sub setup_logic_do_transient_pfts { diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 190745f54e..524b68bfcb 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2750,6 +2750,7 @@ lnd/clm2/surfdata_esmf/NEON/ctsm5.4.0/surfdata_1x1_NEON_TOOL_hist_2000_78pfts_c2 .false. .false. .false. +20.d00 .true. .true. .true. diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 668b23ea6b..3a962dca09 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2893,6 +2893,32 @@ continue run; it is an error for this to be set in a branch run. Furthermore, this setting has no effect in a cold start run. + +Residence time (in years) of the dynbal storage pools. + +Dynamic landunit adjustments change a grid cell's total water and energy content with no +associated physical flux. The resulting spurious gain or loss is put in a storage pool for +each of heat (DYNBAL_HEAT_STORAGE), liquid water (DYNBAL_LIQ_STORAGE) and ice +(DYNBAL_ICE_STORAGE); each pool is then released gradually to a corresponding conservation +flux (EFLX_DYNBAL, QFLX_LIQ_DYNBAL and QFLX_ICE_DYNBAL). This setting controls how +gradually: each pool is drained with a turnover rate of 1 / dynbal_storage_residence_time, +so a larger value spreads the correction over a longer period, giving smaller +instantaneous fluxes. + +The model's own water and energy budgets close every time step regardless of the value +chosen, because the storage pools are counted as part of each grid cell's total water and +energy. What is deferred is the compensating exchange with real, physical state variables: +at any moment, the coupled system carries an outstanding imbalance equal to the current +contents of the dynbal storage pools. Since the pools only approach 0 asymptotically, this +is only fully settled in the limit as time goes to infinity, and larger values (longer +residence times) leave the coupled system further out of balance for longer. + +Must be greater than 0. This can be set to less than 1 year, but should be much longer +than the model time step: if the residence time approaches the length of a time step, the +pools will overshoot and the fluxes will oscillate. + + If TRUE, allow area changes at times other than the year boundary. This should diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 391c733e3a..f0f183485c 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -11,7 +11,6 @@ module dynConsBiogeophysMod ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg - use shr_const_mod , only : SHR_CONST_CDAY use decompMod , only : bounds_type use UrbanParamsType , only : urbanparams_type use EnergyFluxType , only : energyflux_type @@ -40,6 +39,7 @@ module dynConsBiogeophysMod use clm_varcon , only : tfrz, cpliq, hfus, ispval use landunit_varcon , only : istsoil, istice use dynSubgridControlMod , only : get_for_testing_zero_dynbal_fluxes + use dynSubgridControlMod , only : get_dynbal_storage_turnover_rate use filterColMod , only : filter_col_type, col_filter_from_ltypes use clm_time_manager , only : get_step_size_real ! @@ -62,13 +62,6 @@ module dynConsBiogeophysMod ! ! !PRIVATE DATA MEMBERS: - ! Turnover rate of the dynbal storage pools [1/s]. Note that we specify the time in - ! years and then convert it to seconds. (It isn't important to account for leap years - ! here: we just want a value that roughly equates to our desired residence time in - ! years: it doesn't need to be exact.) - real(r8), parameter, private :: dynbal_storage_turnover_rate = & - 1._r8 / (20._r8 * (365._r8 * SHR_CONST_CDAY)) - character(len=*), parameter, private :: sourcefile = & __FILE__ !--------------------------------------------------------------------------- @@ -519,6 +512,7 @@ subroutine dyn_hwcontent_final(bounds, & integer :: i integer :: g ! grid cell index real(r8) :: dtime ! model time step [s] + real(r8) :: turnover_rate ! turnover rate of the dynbal storage pools [1/s] real(r8) :: this_delta_liq(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk or one tracer real(r8) :: delta_liq_bulk(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk water real(r8) :: delta_heat(bounds%begg:bounds%endg) ! change in gridcell heat content @@ -579,9 +573,10 @@ subroutine dyn_hwcontent_final(bounds, & ! the proper sign convention for this pool), then release a portion of the pool as ! this time step's flux. dtime = get_step_size_real() + turnover_rate = get_dynbal_storage_turnover_rate() do g = begg, endg dynbal_heat_storage(g) = dynbal_heat_storage(g) - delta_heat(g) - eflx_dynbal(g) = dynbal_heat_storage(g) * dynbal_storage_turnover_rate + eflx_dynbal(g) = dynbal_heat_storage(g) * turnover_rate dynbal_heat_storage(g) = dynbal_heat_storage(g) - (eflx_dynbal(g) * dtime) end do @@ -617,6 +612,7 @@ subroutine dyn_water_content_final(bounds, & ! !LOCAL VARIABLES: integer :: g real(r8) :: dtime ! model time step [s] + real(r8) :: turnover_rate ! turnover rate of the dynbal storage pools [1/s] real(r8) :: delta_ice(bounds%begg:bounds%endg) ! change in gridcell h2o ice content character(len=*), parameter :: subname = 'dyn_water_content_final' @@ -656,13 +652,14 @@ subroutine dyn_water_content_final(bounds, & ! use the proper sign convention for these pools), then release a portion of each pool ! as this time step's flux. dtime = get_step_size_real() + turnover_rate = get_dynbal_storage_turnover_rate() do g = begg, endg dynbal_liq_storage(g) = dynbal_liq_storage(g) - delta_liq(g) - qflx_liq_dynbal(g) = dynbal_liq_storage(g) * dynbal_storage_turnover_rate + qflx_liq_dynbal(g) = dynbal_liq_storage(g) * turnover_rate dynbal_liq_storage(g) = dynbal_liq_storage(g) - (qflx_liq_dynbal(g) * dtime) dynbal_ice_storage(g) = dynbal_ice_storage(g) - delta_ice(g) - qflx_ice_dynbal(g) = dynbal_ice_storage(g) * dynbal_storage_turnover_rate + qflx_ice_dynbal(g) = dynbal_ice_storage(g) * turnover_rate dynbal_ice_storage(g) = dynbal_ice_storage(g) - (qflx_ice_dynbal(g) * dtime) end do diff --git a/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index c9a6b77824..cec2251cc5 100644 --- a/src/dyn_subgrid/dynSubgridControlMod.F90 +++ b/src/dyn_subgrid/dynSubgridControlMod.F90 @@ -12,6 +12,7 @@ module dynSubgridControlMod ! ! !USES: #include "shr_assert.h" + use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use clm_varctl , only : fname_len @@ -31,6 +32,7 @@ module dynSubgridControlMod public :: get_do_harvest ! return the value of the do_harvest control flag public :: get_do_grossunrep ! return the value of the do_grossunrep control flag public :: get_reset_dynbal_baselines ! return the value of the reset_dynbal_baselines control flag + public :: get_dynbal_storage_turnover_rate ! return the turnover rate of the dynbal storage pools [1/s] public :: get_for_testing_allow_non_annual_changes ! return true if user has requested to allow area changes at times other than the year boundary, for testing purposes public :: get_for_testing_zero_dynbal_fluxes ! return true if user has requested to set the dynbal water and energy fluxes to zero, for testing purposes ! @@ -52,6 +54,11 @@ module dynSubgridControlMod logical :: reset_dynbal_baselines = .false. ! whether to reset baseline values of total column water and energy in the first step of the run + ! Turnover rate of the dynbal storage pools [1/s]. This is derived from the + ! dynbal_storage_residence_time namelist variable, which gives the residence time in + ! years. + real(r8) :: dynbal_storage_turnover_rate = 0._r8 + ! The following is only meant for testing: Whether area changes are allowed at times ! other than the year boundary. This should only arise in some test configurations ! where we artificially create changes more frequently so that we can run short @@ -110,6 +117,7 @@ subroutine read_namelist( NLFilename ) ! Read dyn_subgrid_control namelist variables ! ! !USES: + use shr_const_mod , only : SHR_CONST_CDAY use fileutils , only : getavu, relavu use clm_nlUtilsMod , only : find_nlgroup_name use clm_varctl , only : iulog @@ -130,11 +138,13 @@ subroutine read_namelist( NLFilename ) logical :: do_grossunrep logical :: vars_1dwt_w_time logical :: reset_dynbal_baselines + real(r8) :: dynbal_storage_residence_time logical :: for_testing_allow_non_annual_changes logical :: for_testing_zero_dynbal_fluxes ! other local variables: integer :: nu_nml ! unit for namelist file integer :: nml_error ! namelist i/o error flag + real(r8) :: dynbal_storage_turnover_rate ! turnover rate derived from dynbal_storage_residence_time [1/s] character(len=*), parameter :: subname = 'read_namelist' !----------------------------------------------------------------------- @@ -149,6 +159,7 @@ subroutine read_namelist( NLFilename ) do_grossunrep, & vars_1dwt_w_time, & reset_dynbal_baselines, & + dynbal_storage_residence_time, & for_testing_allow_non_annual_changes, & for_testing_zero_dynbal_fluxes @@ -162,6 +173,7 @@ subroutine read_namelist( NLFilename ) do_grossunrep = .false. vars_1dwt_w_time = .false. reset_dynbal_baselines = .false. + dynbal_storage_residence_time = 0._r8 for_testing_allow_non_annual_changes = .false. for_testing_zero_dynbal_fluxes = .false. @@ -190,9 +202,23 @@ subroutine read_namelist( NLFilename ) call shr_mpi_bcast (do_grossunrep, mpicom) call shr_mpi_bcast (vars_1dwt_w_time, mpicom) call shr_mpi_bcast (reset_dynbal_baselines, mpicom) + call shr_mpi_bcast (dynbal_storage_residence_time, mpicom) call shr_mpi_bcast (for_testing_allow_non_annual_changes, mpicom) call shr_mpi_bcast (for_testing_zero_dynbal_fluxes, mpicom) + ! Convert the residence time to a turnover rate. + if (dynbal_storage_residence_time <= 0._r8) then + write(iulog,*) 'ERROR: dynbal_storage_residence_time must be greater than 0' + write(iulog,*) 'Value given: ', dynbal_storage_residence_time + call endrun(msg='ERROR: dynbal_storage_residence_time must be greater than 0 '// & + errMsg(sourcefile, __LINE__)) + end if + ! It isn't important to account for leap years in this conversion: we just want a + ! value that roughly equates to the desired residence time in years; it doesn't need + ! to be exact.) + dynbal_storage_turnover_rate = & + 1._r8 / (dynbal_storage_residence_time * (365._r8 * SHR_CONST_CDAY)) + dyn_subgrid_control_inst = dyn_subgrid_control_type( & flanduse_timeseries = flanduse_timeseries, & do_transient_pfts = do_transient_pfts, & @@ -203,6 +229,7 @@ subroutine read_namelist( NLFilename ) do_grossunrep = do_grossunrep, & vars_1dwt_w_time = vars_1dwt_w_time, & reset_dynbal_baselines = reset_dynbal_baselines, & + dynbal_storage_turnover_rate = dynbal_storage_turnover_rate, & for_testing_allow_non_annual_changes = for_testing_allow_non_annual_changes, & for_testing_zero_dynbal_fluxes = for_testing_zero_dynbal_fluxes) @@ -222,7 +249,6 @@ subroutine check_namelist_consistency ! Check consistency of namelist settings ! ! !USES: - use shr_kind_mod, only: r8 => shr_kind_r8 use clm_varctl, only : iulog, use_cndv, use_fates, use_cn, use_crop, & n_dom_pfts, n_dom_landunits, collapse_urban, & toosmall_soil, toosmall_crop, toosmall_glacier, & @@ -465,6 +491,20 @@ logical function get_reset_dynbal_baselines() end function get_reset_dynbal_baselines + !----------------------------------------------------------------------- + real(r8) function get_dynbal_storage_turnover_rate() + ! !DESCRIPTION: + ! Return the turnover rate of the dynbal storage pools [1/s] + ! + ! This is derived from the dynbal_storage_residence_time namelist variable. + !----------------------------------------------------------------------- + + SHR_ASSERT_FL(dyn_subgrid_control_inst%initialized, sourcefile, __LINE__) + + get_dynbal_storage_turnover_rate = dyn_subgrid_control_inst%dynbal_storage_turnover_rate + + end function get_dynbal_storage_turnover_rate + !----------------------------------------------------------------------- logical function get_for_testing_allow_non_annual_changes() ! From 6b44163a4e8a0319307aa9be3ffc87101e9b0dcb Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Aug 2026 19:13:15 -0600 Subject: [PATCH 04/20] Reset dynbal storage pools when reset_dynbal_baselines is true This feels conceptually similar to the other resetting done by reset_dynbal_baselines, so it feels like these should be reset at the same time. (This wasn't an issue before the changes on this branch that change the annual dribblers to pools with a residence time, because the annual dribblers would empty themselves by the end of each year.) --- .../namelist_definition_ctsm.xml | 5 +++++ src/dyn_subgrid/dynConsBiogeophysMod.F90 | 19 +++++++++++++++++-- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 3a962dca09..4688184e0c 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2878,6 +2878,11 @@ starting a startup or hybrid run from a partially or entirely spun-up state. This will reset the baseline values based on the state at the start of this run. +This flag also resets the dynbal storage pools to zero. (These pools hold the +still-to-be-distributed water and energy from earlier dynamic column/landunit area +changes. So resetting these to zero means forgetting about these earlier changes, which is +somewhat analogous to the other resetting done by this flag.) + Note that setting this flag can break water and energy conservation! Specifically, any water and energy that has previously been added to or removed from states that contribute to these baselines (currently, (a) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index f0f183485c..9f36688584 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -96,6 +96,13 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icec, filter_icec, & ! the dynbal fluxes; however, it can break conservation. (So, for example, it can be ! done when transitioning from an offline spinup to a coupled run, but it should not ! be done when transitioning from a coupled historical run to a future scenario.) + ! + ! Setting reset_all_baselines also resets the dynbal storage pools to zero. (These + ! pools hold the still-to-be-distributed water and energy from earlier dynamic + ! column/landunit area changes. So resetting these to zero means forgetting about + ! these earlier changes, which is somewhat analogous to the other resetting done by + ! this flag.) + ! ! Other reset_* flags are described below. ! ! !ARGUMENTS: @@ -203,7 +210,9 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & associate( & dynbal_baseline_liq => waterstate_inst%dynbal_baseline_liq_col, & ! Output: [real(r8) (:) ] baseline liquid water content subtracted from each column's total liquid water calculation (mm H2O) - dynbal_baseline_ice => waterstate_inst%dynbal_baseline_ice_col & ! Output: [real(r8) (:) ] baseline ice content subtracted from each column's total ice calculation (mm H2O) + dynbal_baseline_ice => waterstate_inst%dynbal_baseline_ice_col, & ! Output: [real(r8) (:) ] baseline ice content subtracted from each column's total ice calculation (mm H2O) + dynbal_liq_storage => waterstate_inst%dynbal_liq_storage_grc, & ! Output: [real(r8) (:) ] liquid water storage from dynbal adjustments, to be released gradually (mm H2O) + dynbal_ice_storage => waterstate_inst%dynbal_ice_storage_grc & ! Output: [real(r8) (:) ] ice storage from dynbal adjustments, to be released gradually (mm H2O) ) if (reset_all_baselines) then @@ -229,6 +238,9 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, & call set_glacier_baselines(bounds, num_icec, filter_icec, & vals_col = soil_ice_mass_col(bounds%begc:bounds%endc), & baselines_col = dynbal_baseline_ice(bounds%begc:bounds%endc)) + + dynbal_liq_storage(bounds%begg:bounds%endg) = 0._r8 + dynbal_ice_storage(bounds%begg:bounds%endg) = 0._r8 end if if (reset_all_baselines .or. reset_lake_baselines) then @@ -295,7 +307,8 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & !----------------------------------------------------------------------- associate( & - dynbal_baseline_heat => temperature_inst%dynbal_baseline_heat_col & ! Output: [real(r8) (:) ] baseline heat content subtracted from each column's total heat calculation (J/m2) + dynbal_baseline_heat => temperature_inst%dynbal_baseline_heat_col, & ! Output: [real(r8) (:) ] baseline heat content subtracted from each column's total heat calculation (J/m2) + dynbal_heat_storage => temperature_inst%dynbal_heat_storage_grc & ! Output: [real(r8) (:) ] heat storage from dynbal adjustments, to be released gradually (J/m^2) ) if (reset_all_baselines) then @@ -327,6 +340,8 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, & call set_glacier_baselines(bounds, num_icec, filter_icec, & vals_col = soil_heat_col(bounds%begc:bounds%endc), & baselines_col = dynbal_baseline_heat(bounds%begc:bounds%endc)) + + dynbal_heat_storage(bounds%begg:bounds%endg) = 0._r8 end if if (reset_all_baselines .or. reset_lake_baselines) then From bf025d72abee3f0aa1f77d53f81ba532a7a6df62 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Mon, 31 Aug 2026 16:39:02 -0600 Subject: [PATCH 05/20] Update tech note for changes on this branch --- .../Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst b/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst index ad29c64540..f1e2f0156e 100644 --- a/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst +++ b/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst @@ -47,7 +47,7 @@ When subgrid areas change, the water and energy states remain unchanged on a per For example, consider a gridcell with two columns: column 1 has a water mass of 1 kg m\ :sup:`-2` and column 2 has a water mass of 2 kg m\ :sup:`-2` for a given water state variable, where these are expressed per unit column area. If column 1 increases in area at the expense of column 2, then column 1 will still have a water mass of 1 kg m\ :sup:`-2`, but now expressed over the new column area. This results in a decrease in the total gridcell water content. -Water and energy are conserved by summing up the total water and energy content of each gridcell before and after a change in area. Differences in liquid and ice water content are balanced by liquid and ice runoff terms, which can be either positive or negative. (Negative runoff is effectively a withdrawal of water from the ocean.) Differences in energy content are balanced by a sensible heat flux term, which again can be either positive or negative. These balancing fluxes are spread evenly throughout the following year. +Water and energy are conserved by summing up the total water and energy content of each gridcell before and after a change in area. Differences in liquid and ice water content are balanced by liquid and ice runoff terms, which can be either positive or negative. (Negative runoff is effectively a withdrawal of water from the ocean.) Differences in energy content are balanced by a sensible heat flux term, which again can be either positive or negative. These balancing fluxes are released gradually from temporary storage pools with 20-year residence times. There is a special case when a given crop column type newly comes into existence - for example, when temperate corn first comes into existence in a gridcell. In this case, the column's below-ground temperature and water states are copied from the natural vegetated column in its gridcell, so that these state variables begin in a close-to-spun-up state. Other state variables (most of which spin up relatively quickly) begin at their cold start initialization values. This initialization is not necessary for the two other land unit types that currently can grow - natural vegetation and glacier: Those land unit types are always active, even when they have zero area on the gridcell, so their state variables will be spun up immediately when they come into existence. After this initialization, the conservation code described above takes effect. From 1eee35dda2ecb7f2dcaa8c59de38d14bc05460b7 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 2 Sep 2026 17:20:55 -0600 Subject: [PATCH 06/20] Add a note in namelist definition --- bld/namelist_files/namelist_definition_ctsm.xml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 4688184e0c..74b2b250b1 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2922,6 +2922,12 @@ residence times) leave the coupled system further out of balance for longer. Must be greater than 0. This can be set to less than 1 year, but should be much longer than the model time step: if the residence time approaches the length of a time step, the pools will overshoot and the fluxes will oscillate. + +(Prior to the introduction of this namelist option, the balancing fluxes were spread +evenly throughout the following year. There is no way to exactly recover the old CLM5 +behavior, but it can be approximated by setting this option to something between 0.5 and +1.0, depending on the behaviors you want to recover - peak fluxes, effective residence +time, etc.) Date: Thu, 3 Sep 2026 15:21:58 -0600 Subject: [PATCH 07/20] dynbal_storage_residence_time: maintain close to old behavior for clm5_0 It's impossible to exactly recover the old behavior, but using a residence time somewhere between 0.5 year and 1.0 year gives something that roughly approximates the old behavior. So make that the default for clm5_0. Co-authored-by: Erik Kluzek --- bld/namelist_files/namelist_defaults_ctsm.xml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 524b68bfcb..637710ccdb 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2750,7 +2750,8 @@ lnd/clm2/surfdata_esmf/NEON/ctsm5.4.0/surfdata_1x1_NEON_TOOL_hist_2000_78pfts_c2 .false. .false. .false. -20.d00 +20.d00 +0.75d00 .true. .true. .true. From e5768b97019885aef6ddfca93f373c897d6ab256 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 3 Sep 2026 15:58:41 -0600 Subject: [PATCH 08/20] Don't use hard-coded 365 days per year Instead of a hard-coded 365, use get_average_days_per_year (both to avoid hard-coding 365, and to get a more correct value if we're running with a Gregorian calendar). --- src/dyn_subgrid/dynSubgridControlMod.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index cec2251cc5..a8a81b0dde 100644 --- a/src/dyn_subgrid/dynSubgridControlMod.F90 +++ b/src/dyn_subgrid/dynSubgridControlMod.F90 @@ -16,6 +16,7 @@ module dynSubgridControlMod use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use clm_varctl , only : fname_len + use clm_time_manager , only : get_average_days_per_year ! implicit none private @@ -213,11 +214,8 @@ subroutine read_namelist( NLFilename ) call endrun(msg='ERROR: dynbal_storage_residence_time must be greater than 0 '// & errMsg(sourcefile, __LINE__)) end if - ! It isn't important to account for leap years in this conversion: we just want a - ! value that roughly equates to the desired residence time in years; it doesn't need - ! to be exact.) - dynbal_storage_turnover_rate = & - 1._r8 / (dynbal_storage_residence_time * (365._r8 * SHR_CONST_CDAY)) + dynbal_storage_turnover_rate = 1._r8 / & + (dynbal_storage_residence_time * (get_average_days_per_year() * SHR_CONST_CDAY)) dyn_subgrid_control_inst = dyn_subgrid_control_type( & flanduse_timeseries = flanduse_timeseries, & From 54ff1feab4ba6cac221300240cbb10703f9e17f0 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 3 Sep 2026 16:14:19 -0600 Subject: [PATCH 09/20] Remove a note from namelist definition Now that namelist defaults recovers something close to the old clm5_0 behavior, I feel that this extra text in namelist definition is unnecessary. --- bld/namelist_files/namelist_definition_ctsm.xml | 6 ------ 1 file changed, 6 deletions(-) diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 74b2b250b1..4688184e0c 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2922,12 +2922,6 @@ residence times) leave the coupled system further out of balance for longer. Must be greater than 0. This can be set to less than 1 year, but should be much longer than the model time step: if the residence time approaches the length of a time step, the pools will overshoot and the fluxes will oscillate. - -(Prior to the introduction of this namelist option, the balancing fluxes were spread -evenly throughout the following year. There is no way to exactly recover the old CLM5 -behavior, but it can be approximated by setting this option to something between 0.5 and -1.0, depending on the behaviors you want to recover - peak fluxes, effective residence -time, etc.) Date: Thu, 3 Sep 2026 17:09:02 -0600 Subject: [PATCH 10/20] Add a comment describing clm5 default for dynbal_storage_residence_time --- bld/namelist_files/namelist_defaults_ctsm.xml | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 637710ccdb..ba67f56f96 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2750,8 +2750,19 @@ lnd/clm2/surfdata_esmf/NEON/ctsm5.4.0/surfdata_1x1_NEON_TOOL_hist_2000_78pfts_c2 .false. .false. .false. + + 20.d00 0.75d00 + .true. .true. .true. From 3f516e8d89cad3a2d2e50f52a66b6806059cf034 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 3 Sep 2026 17:15:42 -0600 Subject: [PATCH 11/20] More comments on default values for dynbal_storage_residence_time --- bld/namelist_files/namelist_defaults_ctsm.xml | 26 ++++++++++++------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index ba67f56f96..a40a82a56b 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2751,15 +2751,23 @@ lnd/clm2/surfdata_esmf/NEON/ctsm5.4.0/surfdata_1x1_NEON_TOOL_hist_2000_78pfts_c2 .false. .false. - + 20.d00 0.75d00 From 219555257f5de41b0806376443198f7b37662878 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 4 Sep 2026 11:07:56 -0600 Subject: [PATCH 12/20] Move some description out of namelist definition into tech note --- bld/namelist_files/namelist_definition_ctsm.xml | 11 ++--------- .../CLM50_Tech_Note_Transient_Landcover.rst | 2 +- 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 4688184e0c..b556db97c6 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2909,15 +2909,8 @@ each of heat (DYNBAL_HEAT_STORAGE), liquid water (DYNBAL_LIQ_STORAGE) and ice flux (EFLX_DYNBAL, QFLX_LIQ_DYNBAL and QFLX_ICE_DYNBAL). This setting controls how gradually: each pool is drained with a turnover rate of 1 / dynbal_storage_residence_time, so a larger value spreads the correction over a longer period, giving smaller -instantaneous fluxes. - -The model's own water and energy budgets close every time step regardless of the value -chosen, because the storage pools are counted as part of each grid cell's total water and -energy. What is deferred is the compensating exchange with real, physical state variables: -at any moment, the coupled system carries an outstanding imbalance equal to the current -contents of the dynbal storage pools. Since the pools only approach 0 asymptotically, this -is only fully settled in the limit as time goes to infinity, and larger values (longer -residence times) leave the coupled system further out of balance for longer. +instantaneous fluxes but also leaving the coupled system further out of balance for +longer. Must be greater than 0. This can be set to less than 1 year, but should be much longer than the model time step: if the residence time approaches the length of a time step, the diff --git a/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst b/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst index f1e2f0156e..12f081abbf 100644 --- a/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst +++ b/doc/source/tech_note/Transient_Landcover/CLM50_Tech_Note_Transient_Landcover.rst @@ -47,7 +47,7 @@ When subgrid areas change, the water and energy states remain unchanged on a per For example, consider a gridcell with two columns: column 1 has a water mass of 1 kg m\ :sup:`-2` and column 2 has a water mass of 2 kg m\ :sup:`-2` for a given water state variable, where these are expressed per unit column area. If column 1 increases in area at the expense of column 2, then column 1 will still have a water mass of 1 kg m\ :sup:`-2`, but now expressed over the new column area. This results in a decrease in the total gridcell water content. -Water and energy are conserved by summing up the total water and energy content of each gridcell before and after a change in area. Differences in liquid and ice water content are balanced by liquid and ice runoff terms, which can be either positive or negative. (Negative runoff is effectively a withdrawal of water from the ocean.) Differences in energy content are balanced by a sensible heat flux term, which again can be either positive or negative. These balancing fluxes are released gradually from temporary storage pools with 20-year residence times. +Water and energy are conserved by summing up the total water and energy content of each gridcell before and after a change in area. Differences in liquid and ice water content are balanced by liquid and ice runoff terms, which can be either positive or negative. (Negative runoff is effectively a withdrawal of water from the ocean.) Differences in energy content are balanced by a sensible heat flux term, which again can be either positive or negative. These balancing fluxes are released gradually from temporary storage pools with 20-year residence times (with separate storage pools for each of heat, liquid water and ice). CLM's water and energy budgets close every time step because the storage pools are counted as part of each grid cell's total water and energy. What is deferred is the compensating exchange with real, physical state variables: at any moment, the coupled system carries an outstanding imbalance equal to the current contents of the storage pools. Since the pools only approach 0 asymptotically, this is only fully settled in the limit as time goes to infinity. There is a special case when a given crop column type newly comes into existence - for example, when temperate corn first comes into existence in a gridcell. In this case, the column's below-ground temperature and water states are copied from the natural vegetated column in its gridcell, so that these state variables begin in a close-to-spun-up state. Other state variables (most of which spin up relatively quickly) begin at their cold start initialization values. This initialization is not necessary for the two other land unit types that currently can grow - natural vegetation and glacier: Those land unit types are always active, even when they have zero area on the gridcell, so their state variables will be spun up immediately when they come into existence. After this initialization, the conservation code described above takes effect. From 3fe863a21317093786883d37d5e3e98fb6e6e5a9 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 4 Sep 2026 12:28:06 -0600 Subject: [PATCH 13/20] Move dynbal_storage_residence_time NL read into dynConsBiogeophysMod The problem with the previous placement is that dynSubgridControlMod is initialized before the calendar is initialized. This new placement allows us to initialize this variable after the calendar is initialized. It also keeps the namelist variable local to where it's needed. --- bld/CLMBuildNamelist.pm | 17 ++- .../namelist_definition_ctsm.xml | 42 +++---- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 103 ++++++++++++++++-- src/dyn_subgrid/dynSubgridControlMod.F90 | 40 +------ src/main/clm_initializeMod.F90 | 6 + 5 files changed, 140 insertions(+), 68 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 560e1bdca7..3ce5fe24a5 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1931,6 +1931,11 @@ sub process_namelist_inline_logic { # namelist group: clm_temperature_inparm # ########################################## setup_logic_coldstart_temp($opts,$nl_flags, $definition, $defaults, $nl); + + ############################################## + # namelist group: dyn_cons_biogeophys_inparm # + ############################################## + setup_logic_dyn_cons_biogeophys($opts, $nl_flags, $definition, $defaults, $nl); } #------------------------------------------------------------------------------- @@ -2827,6 +2832,15 @@ sub setup_logic_dynamic_subgrid { &remove_leading_and_trailing_quotes($nl_flags->{'clm_start_type'}) eq "branch") { $log->fatal_error("reset_dynbal_baselines has no effect in a branch run"); } +} + +#------------------------------------------------------------------------------- + +sub setup_logic_dyn_cons_biogeophys { + # + # Options controlling conservation of water and energy with dynamic land cover + # + my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'dynbal_storage_residence_time'); my $dynbal_storage_residence_time = $nl->get_value('dynbal_storage_residence_time'); @@ -5419,7 +5433,8 @@ sub write_output_files { soil_resis_inparm bgc_shared canopyfluxes_inparm aerosol clmu_inparm clm_soilstate_inparm clm_nitrogen clm_snowhydrology_inparm hillslope_hydrology_inparm hillslope_properties_inparm cnprecision_inparm clm_glacier_behavior crop_inparm irrigation_inparm - surfacealbedo_inparm water_tracers_inparm tillage_inparm); + surfacealbedo_inparm water_tracers_inparm tillage_inparm + dyn_cons_biogeophys_inparm); #@groups = qw(clm_inparm clm_canopyhydrology_inparm clm_soilhydrology_inparm # finidat_consistency_checks dynpft_consistency_checks); diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index b556db97c6..cc20e00fc2 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2898,25 +2898,6 @@ continue run; it is an error for this to be set in a branch run. Furthermore, this setting has no effect in a cold start run. - -Residence time (in years) of the dynbal storage pools. - -Dynamic landunit adjustments change a grid cell's total water and energy content with no -associated physical flux. The resulting spurious gain or loss is put in a storage pool for -each of heat (DYNBAL_HEAT_STORAGE), liquid water (DYNBAL_LIQ_STORAGE) and ice -(DYNBAL_ICE_STORAGE); each pool is then released gradually to a corresponding conservation -flux (EFLX_DYNBAL, QFLX_LIQ_DYNBAL and QFLX_ICE_DYNBAL). This setting controls how -gradually: each pool is drained with a turnover rate of 1 / dynbal_storage_residence_time, -so a larger value spreads the correction over a longer period, giving smaller -instantaneous fluxes but also leaving the coupled system further out of balance for -longer. - -Must be greater than 0. This can be set to less than 1 year, but should be much longer -than the model time step: if the residence time approaches the length of a time step, the -pools will overshoot and the fluxes will oscillate. - - If TRUE, allow area changes at times other than the year boundary. This should @@ -2945,6 +2926,29 @@ outcome as vars_1dwt_w_time = .true.. Use this flag if you wish to change FALSE to TRUE when run_has_transient_landcover = .false.. + + + + + +Residence time (in years) of the dynbal storage pools. + +Dynamic landunit adjustments change a grid cell's total water and energy content with no +associated physical flux. The resulting spurious gain or loss is put in a storage pool for +each of heat (DYNBAL_HEAT_STORAGE), liquid water (DYNBAL_LIQ_STORAGE) and ice +(DYNBAL_ICE_STORAGE); each pool is then released gradually to a corresponding conservation +flux (EFLX_DYNBAL, QFLX_LIQ_DYNBAL and QFLX_ICE_DYNBAL). This setting controls how +gradually: each pool is drained with a turnover rate of 1 / dynbal_storage_residence_time, +so a larger value spreads the correction over a longer period, giving smaller +instantaneous fluxes but also leaving the coupled system further out of balance for +longer. + +Must be greater than 0. This can be set to less than 1 year, but should be much longer +than the model time step: if the residence time approaches the length of a time step, the +pools will overshoot and the fluxes will oscillate. + + diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 9f36688584..2f34d2b459 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -39,14 +39,14 @@ module dynConsBiogeophysMod use clm_varcon , only : tfrz, cpliq, hfus, ispval use landunit_varcon , only : istsoil, istice use dynSubgridControlMod , only : get_for_testing_zero_dynbal_fluxes - use dynSubgridControlMod , only : get_dynbal_storage_turnover_rate use filterColMod , only : filter_col_type, col_filter_from_ltypes - use clm_time_manager , only : get_step_size_real + use clm_time_manager , only : get_step_size_real, get_average_days_per_year ! ! !PUBLIC MEMBER FUNCTIONS: implicit none private ! + public :: dynConsBiogeophys_readnl ! read namelist variables public :: dyn_hwcontent_set_baselines ! set start-of-run baseline values for heat and water content in some columns public :: dyn_hwcontent_init ! compute grid-level heat and water content, before land cover change public :: dyn_hwcontent_final ! compute grid-level heat and water content and dynbal fluxes after land cover change @@ -62,12 +62,99 @@ module dynConsBiogeophysMod ! ! !PRIVATE DATA MEMBERS: + ! Turnover rate of the dynbal storage pools [1/s]. This is derived from the + ! dynbal_storage_residence_time namelist variable, which gives the residence time in + ! years. + real(r8), private :: dynbal_storage_turnover_rate + + ! Whether the namelist has been read (and thus whether dynbal_storage_turnover_rate has + ! been set) + logical, private :: namelist_read = .false. + character(len=*), parameter, private :: sourcefile = & __FILE__ !--------------------------------------------------------------------------- contains + !----------------------------------------------------------------------- + subroutine dynConsBiogeophys_readnl(NLFilename) + ! + ! !DESCRIPTION: + ! Read the dyn_cons_biogeophys_inparm namelist. + ! + ! Note that this must be called AFTER the time manager has been initialized, because + ! converting dynbal_storage_residence_time to a turnover rate depends on the model + ! calendar. + ! + ! !USES: + use shr_const_mod , only : SHR_CONST_CDAY + use abortutils , only : endrun + use fileutils , only : getavu, relavu + use clm_nlUtilsMod , only : find_nlgroup_name + use clm_varctl , only : iulog + use spmdMod , only : masterproc, mpicom + use shr_mpi_mod , only : shr_mpi_bcast + ! + ! !ARGUMENTS: + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! !LOCAL VARIABLES: + ! temporary variable corresponding to the namelist variable: + real(r8) :: dynbal_storage_residence_time ! residence time of the dynbal storage pools [years] + ! other local variables: + integer :: nu_nml ! unit for namelist file + integer :: nml_error ! namelist i/o error flag + + character(len=*), parameter :: subname = 'dynConsBiogeophys_readnl' + !----------------------------------------------------------------------- + + namelist /dyn_cons_biogeophys_inparm/ & + dynbal_storage_residence_time + + dynbal_storage_residence_time = 0._r8 + + if (masterproc) then + nu_nml = getavu() + open( nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) + call find_nlgroup_name(nu_nml, 'dyn_cons_biogeophys_inparm', status=nml_error) + if (nml_error == 0) then + read(nu_nml, nml=dyn_cons_biogeophys_inparm, iostat=nml_error) + if (nml_error /= 0) then + call endrun(msg='ERROR reading dyn_cons_biogeophys_inparm namelist'//errMsg(sourcefile, __LINE__)) + end if + else + call endrun(msg='ERROR finding dyn_cons_biogeophys_inparm namelist'//errMsg(sourcefile, __LINE__)) + end if + close(nu_nml) + call relavu( nu_nml ) + endif + + call shr_mpi_bcast (dynbal_storage_residence_time, mpicom) + + ! Convert the residence time to a turnover rate. Note that we check for a valid + ! residence time before doing this conversion, since the conversion divides by the + ! residence time. + if (dynbal_storage_residence_time <= 0._r8) then + write(iulog,*) 'ERROR: dynbal_storage_residence_time must be greater than 0' + write(iulog,*) 'Value given: ', dynbal_storage_residence_time + call endrun(msg='ERROR: dynbal_storage_residence_time must be greater than 0 '// & + errMsg(sourcefile, __LINE__)) + end if + dynbal_storage_turnover_rate = 1._r8 / & + (dynbal_storage_residence_time * (get_average_days_per_year() * SHR_CONST_CDAY)) + + namelist_read = .true. + + if (masterproc) then + write(iulog,*) ' ' + write(iulog,*) 'dyn_cons_biogeophys_inparm settings:' + write(iulog,nml=dyn_cons_biogeophys_inparm) + write(iulog,*) ' ' + end if + + end subroutine dynConsBiogeophys_readnl + !----------------------------------------------------------------------- subroutine dyn_hwcontent_set_baselines(bounds, num_icec, filter_icec, & num_lakec, filter_lakec, & @@ -527,12 +614,13 @@ subroutine dyn_hwcontent_final(bounds, & integer :: i integer :: g ! grid cell index real(r8) :: dtime ! model time step [s] - real(r8) :: turnover_rate ! turnover rate of the dynbal storage pools [1/s] real(r8) :: this_delta_liq(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk or one tracer real(r8) :: delta_liq_bulk(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk water real(r8) :: delta_heat(bounds%begg:bounds%endg) ! change in gridcell heat content !--------------------------------------------------------------------------- + SHR_ASSERT_FL(namelist_read, sourcefile, __LINE__) + associate( & begg => bounds%begg, & endg => bounds%endg, & @@ -588,10 +676,9 @@ subroutine dyn_hwcontent_final(bounds, & ! the proper sign convention for this pool), then release a portion of the pool as ! this time step's flux. dtime = get_step_size_real() - turnover_rate = get_dynbal_storage_turnover_rate() do g = begg, endg dynbal_heat_storage(g) = dynbal_heat_storage(g) - delta_heat(g) - eflx_dynbal(g) = dynbal_heat_storage(g) * turnover_rate + eflx_dynbal(g) = dynbal_heat_storage(g) * dynbal_storage_turnover_rate dynbal_heat_storage(g) = dynbal_heat_storage(g) - (eflx_dynbal(g) * dtime) end do @@ -627,7 +714,6 @@ subroutine dyn_water_content_final(bounds, & ! !LOCAL VARIABLES: integer :: g real(r8) :: dtime ! model time step [s] - real(r8) :: turnover_rate ! turnover rate of the dynbal storage pools [1/s] real(r8) :: delta_ice(bounds%begg:bounds%endg) ! change in gridcell h2o ice content character(len=*), parameter :: subname = 'dyn_water_content_final' @@ -667,14 +753,13 @@ subroutine dyn_water_content_final(bounds, & ! use the proper sign convention for these pools), then release a portion of each pool ! as this time step's flux. dtime = get_step_size_real() - turnover_rate = get_dynbal_storage_turnover_rate() do g = begg, endg dynbal_liq_storage(g) = dynbal_liq_storage(g) - delta_liq(g) - qflx_liq_dynbal(g) = dynbal_liq_storage(g) * turnover_rate + qflx_liq_dynbal(g) = dynbal_liq_storage(g) * dynbal_storage_turnover_rate dynbal_liq_storage(g) = dynbal_liq_storage(g) - (qflx_liq_dynbal(g) * dtime) dynbal_ice_storage(g) = dynbal_ice_storage(g) - delta_ice(g) - qflx_ice_dynbal(g) = dynbal_ice_storage(g) * turnover_rate + qflx_ice_dynbal(g) = dynbal_ice_storage(g) * dynbal_storage_turnover_rate dynbal_ice_storage(g) = dynbal_ice_storage(g) - (qflx_ice_dynbal(g) * dtime) end do diff --git a/src/dyn_subgrid/dynSubgridControlMod.F90 b/src/dyn_subgrid/dynSubgridControlMod.F90 index a8a81b0dde..f8e59d4762 100644 --- a/src/dyn_subgrid/dynSubgridControlMod.F90 +++ b/src/dyn_subgrid/dynSubgridControlMod.F90 @@ -16,7 +16,6 @@ module dynSubgridControlMod use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use clm_varctl , only : fname_len - use clm_time_manager , only : get_average_days_per_year ! implicit none private @@ -33,7 +32,6 @@ module dynSubgridControlMod public :: get_do_harvest ! return the value of the do_harvest control flag public :: get_do_grossunrep ! return the value of the do_grossunrep control flag public :: get_reset_dynbal_baselines ! return the value of the reset_dynbal_baselines control flag - public :: get_dynbal_storage_turnover_rate ! return the turnover rate of the dynbal storage pools [1/s] public :: get_for_testing_allow_non_annual_changes ! return true if user has requested to allow area changes at times other than the year boundary, for testing purposes public :: get_for_testing_zero_dynbal_fluxes ! return true if user has requested to set the dynbal water and energy fluxes to zero, for testing purposes ! @@ -55,11 +53,6 @@ module dynSubgridControlMod logical :: reset_dynbal_baselines = .false. ! whether to reset baseline values of total column water and energy in the first step of the run - ! Turnover rate of the dynbal storage pools [1/s]. This is derived from the - ! dynbal_storage_residence_time namelist variable, which gives the residence time in - ! years. - real(r8) :: dynbal_storage_turnover_rate = 0._r8 - ! The following is only meant for testing: Whether area changes are allowed at times ! other than the year boundary. This should only arise in some test configurations ! where we artificially create changes more frequently so that we can run short @@ -118,7 +111,6 @@ subroutine read_namelist( NLFilename ) ! Read dyn_subgrid_control namelist variables ! ! !USES: - use shr_const_mod , only : SHR_CONST_CDAY use fileutils , only : getavu, relavu use clm_nlUtilsMod , only : find_nlgroup_name use clm_varctl , only : iulog @@ -139,14 +131,12 @@ subroutine read_namelist( NLFilename ) logical :: do_grossunrep logical :: vars_1dwt_w_time logical :: reset_dynbal_baselines - real(r8) :: dynbal_storage_residence_time logical :: for_testing_allow_non_annual_changes logical :: for_testing_zero_dynbal_fluxes ! other local variables: integer :: nu_nml ! unit for namelist file integer :: nml_error ! namelist i/o error flag - real(r8) :: dynbal_storage_turnover_rate ! turnover rate derived from dynbal_storage_residence_time [1/s] - + character(len=*), parameter :: subname = 'read_namelist' !----------------------------------------------------------------------- @@ -160,7 +150,6 @@ subroutine read_namelist( NLFilename ) do_grossunrep, & vars_1dwt_w_time, & reset_dynbal_baselines, & - dynbal_storage_residence_time, & for_testing_allow_non_annual_changes, & for_testing_zero_dynbal_fluxes @@ -174,7 +163,6 @@ subroutine read_namelist( NLFilename ) do_grossunrep = .false. vars_1dwt_w_time = .false. reset_dynbal_baselines = .false. - dynbal_storage_residence_time = 0._r8 for_testing_allow_non_annual_changes = .false. for_testing_zero_dynbal_fluxes = .false. @@ -203,20 +191,9 @@ subroutine read_namelist( NLFilename ) call shr_mpi_bcast (do_grossunrep, mpicom) call shr_mpi_bcast (vars_1dwt_w_time, mpicom) call shr_mpi_bcast (reset_dynbal_baselines, mpicom) - call shr_mpi_bcast (dynbal_storage_residence_time, mpicom) call shr_mpi_bcast (for_testing_allow_non_annual_changes, mpicom) call shr_mpi_bcast (for_testing_zero_dynbal_fluxes, mpicom) - ! Convert the residence time to a turnover rate. - if (dynbal_storage_residence_time <= 0._r8) then - write(iulog,*) 'ERROR: dynbal_storage_residence_time must be greater than 0' - write(iulog,*) 'Value given: ', dynbal_storage_residence_time - call endrun(msg='ERROR: dynbal_storage_residence_time must be greater than 0 '// & - errMsg(sourcefile, __LINE__)) - end if - dynbal_storage_turnover_rate = 1._r8 / & - (dynbal_storage_residence_time * (get_average_days_per_year() * SHR_CONST_CDAY)) - dyn_subgrid_control_inst = dyn_subgrid_control_type( & flanduse_timeseries = flanduse_timeseries, & do_transient_pfts = do_transient_pfts, & @@ -227,7 +204,6 @@ subroutine read_namelist( NLFilename ) do_grossunrep = do_grossunrep, & vars_1dwt_w_time = vars_1dwt_w_time, & reset_dynbal_baselines = reset_dynbal_baselines, & - dynbal_storage_turnover_rate = dynbal_storage_turnover_rate, & for_testing_allow_non_annual_changes = for_testing_allow_non_annual_changes, & for_testing_zero_dynbal_fluxes = for_testing_zero_dynbal_fluxes) @@ -489,20 +465,6 @@ logical function get_reset_dynbal_baselines() end function get_reset_dynbal_baselines - !----------------------------------------------------------------------- - real(r8) function get_dynbal_storage_turnover_rate() - ! !DESCRIPTION: - ! Return the turnover rate of the dynbal storage pools [1/s] - ! - ! This is derived from the dynbal_storage_residence_time namelist variable. - !----------------------------------------------------------------------- - - SHR_ASSERT_FL(dyn_subgrid_control_inst%initialized, sourcefile, __LINE__) - - get_dynbal_storage_turnover_rate = dyn_subgrid_control_inst%dynbal_storage_turnover_rate - - end function get_dynbal_storage_turnover_rate - !----------------------------------------------------------------------- logical function get_for_testing_allow_non_annual_changes() ! diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index 8505276174..ae16e81616 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -160,6 +160,7 @@ subroutine initialize2(ni,nj, currtime) use CIsoAtmTimeseriesMod , only : CIsoAtmReadNML use DaylengthMod , only : InitDaylength use dynSubgridDriverMod , only : dynSubgrid_init + use dynConsBiogeophysMod , only : dynConsBiogeophys_readnl use dynConsBiogeophysMod , only : dyn_hwcontent_set_baselines use fileutils , only : getfil use initInterpMod , only : initInterp @@ -380,6 +381,11 @@ subroutine initialize2(ni,nj, currtime) ! Initialize Balance checking (after time-manager) call BalanceCheckInit() + ! Read the dyn_cons_biogeophys_inparm namelist. This must be done after the time + ! manager is initialized, because converting dynbal_storage_residence_time to a + ! turnover rate depends on the model calendar. + call dynConsBiogeophys_readnl(NLFilename) + ! History file variables if (use_cn .and. .not. use_noio ) then call hist_addfld1d (fname='DAYL', units='s', & From 1e56fd67705cf209780397ed37996b57bff0d311 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 4 Sep 2026 15:37:44 -0600 Subject: [PATCH 14/20] Fix test_dyn_cons_biogeophys unit tests Also, test the new zeroing of the dynbal storage terms in this unit test. --- .../test_dyn_cons_biogeophys.pf | 23 ++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf b/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf index d991f0cc54..98e678fb2c 100644 --- a/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf +++ b/src/dyn_subgrid/test/dynConsBiogeophys_test/test_dyn_cons_biogeophys.pf @@ -6,7 +6,7 @@ module test_dyn_cons_biogeophys use dynConsBiogeophysMod use shr_kind_mod , only : r8 => shr_kind_r8 use unittestSubgridMod - use unittestArrayMod, only : col_array + use unittestArrayMod, only : col_array, grc_array use unittestFilterBuilderMod, only : filter_from_range, filter_empty use unittestWaterTypeFactory, only : unittest_water_type_factory_type use clm_varpar, only : nlevgrnd, nlevsno, maxpatch_glc @@ -73,6 +73,9 @@ contains if (associated(this%temperature_inst%dynbal_baseline_heat_col)) then deallocate(this%temperature_inst%dynbal_baseline_heat_col) end if + if (associated(this%temperature_inst%dynbal_heat_storage_grc)) then + deallocate(this%temperature_inst%dynbal_heat_storage_grc) + end if if (associated(this%urbanparams_inst%nlev_improad)) then deallocate(this%urbanparams_inst%nlev_improad) @@ -109,6 +112,7 @@ contains real(r8), allocatable :: expected_vals_heat_col(:) real(r8), allocatable :: ignored_heatliq_col(:) real(r8), allocatable :: ignored_cvliq_col(:) + real(r8), allocatable :: expected_zero_grc(:) ! ------------------------------------------------------------------------ ! Create subgrid structure. @@ -178,6 +182,10 @@ contains ! Initialize necessary variables in temperature_inst allocate(this%temperature_inst%t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) allocate(this%temperature_inst%dynbal_baseline_heat_col(bounds%begc:bounds%endc)) + allocate(this%temperature_inst%dynbal_heat_storage_grc(bounds%begg:bounds%endg)) + ! Give the storage pool an arbitrary non-zero value, so that we can verify below that + ! reset_all_baselines resets it to 0 + this%temperature_inst%dynbal_heat_storage_grc(bounds%begg:bounds%endg) = 1000._r8 do c = bounds%begc, bounds%endc ! Use a different value of temperature for each column if (c == glc_col) then @@ -203,6 +211,11 @@ contains t_soisno_col = this%temperature_inst%t_soisno_col(bounds%begc:bounds%endc,:), & watsat_col = this%soilstate_inst%watsat_col(bounds%begc:bounds%endc,:)) + ! Give the water storage pools arbitrary non-zero values, so that we can verify below + ! that reset_all_baselines resets them to 0 + this%water_inst%waterstatebulk_inst%dynbal_liq_storage_grc(bounds%begg:bounds%endg) = 100._r8 + this%water_inst%waterstatebulk_inst%dynbal_ice_storage_grc(bounds%begg:bounds%endg) = 200._r8 + ! ------------------------------------------------------------------------ ! Call the routine we're testing ! ------------------------------------------------------------------------ @@ -258,6 +271,14 @@ contains baselines_col = this%temperature_inst%dynbal_baseline_heat_col, & msg = 'heat') + ! Setting reset_all_baselines should also have reset the dynbal storage pools to 0 + expected_zero_grc = grc_array(0._r8) + associate(waterstatebulk_inst => this%water_inst%waterstatebulk_inst) + @assertEqual(expected_zero_grc, waterstatebulk_inst%dynbal_liq_storage_grc, message='liq storage') + @assertEqual(expected_zero_grc, waterstatebulk_inst%dynbal_ice_storage_grc, message='ice storage') + end associate + @assertEqual(expected_zero_grc, this%temperature_inst%dynbal_heat_storage_grc, message='heat storage') + contains subroutine assertBaselines(expected_vals_col, baselines_col, msg) real(r8), intent(in) :: expected_vals_col(bounds%begc:) From cf37d1a25ec4633b2e8baffcb00e361d532ad2b2 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 4 Sep 2026 18:55:04 -0600 Subject: [PATCH 15/20] Add some tests to expected fails When the RUN fails, the COMPARE_base_rest is left in a PEND state --- cime_config/testdefs/ExpectedTestFails.xml | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index 6eb0b2fcdb..407e9c0aa8 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -140,6 +140,11 @@ FAIL #3788 + + + PEND + #3788 + @@ -255,6 +260,11 @@ FAIL #2261 + + + PEND + #2261 + From e10c2e0eb1e0966222fe15a34ef3780b17f691f8 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Sat, 5 Sep 2026 09:57:12 -0600 Subject: [PATCH 16/20] Refactor namelist handling Follow more recent patterns, as exemplified by CIsoAtmReadNML Changes made by Claude, reviewed carefully by myself --- src/dyn_subgrid/dynConsBiogeophysMod.F90 | 68 +++++++++++------------- 1 file changed, 32 insertions(+), 36 deletions(-) diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 2f34d2b459..12d451e037 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -88,25 +88,23 @@ subroutine dynConsBiogeophys_readnl(NLFilename) ! calendar. ! ! !USES: - use shr_const_mod , only : SHR_CONST_CDAY - use abortutils , only : endrun - use fileutils , only : getavu, relavu - use clm_nlUtilsMod , only : find_nlgroup_name - use clm_varctl , only : iulog - use spmdMod , only : masterproc, mpicom - use shr_mpi_mod , only : shr_mpi_bcast - ! + use shr_const_mod, only : SHR_CONST_CDAY + use shr_nl_mod , only : shr_nl_find_group_name + use abortutils , only : endrun + use clm_varctl , only : iulog + use spmdMod , only : masterproc, mpicom + use shr_mpi_mod , only : shr_mpi_bcast + ! !ARGUMENTS: - character(len=*), intent(in) :: NLFilename ! Namelist filename - ! + character(len=*), intent(in) :: NLFilename ! Namelist filename to read + ! !LOCAL VARIABLES: ! temporary variable corresponding to the namelist variable: real(r8) :: dynbal_storage_residence_time ! residence time of the dynbal storage pools [years] ! other local variables: - integer :: nu_nml ! unit for namelist file - integer :: nml_error ! namelist i/o error flag - - character(len=*), parameter :: subname = 'dynConsBiogeophys_readnl' + integer :: ierr ! error code + integer :: unitn ! unit for namelist file + character(len=*), parameter :: nml_name = 'dyn_cons_biogeophys_inparm' ! MUST agree with name in namelist and read !----------------------------------------------------------------------- namelist /dyn_cons_biogeophys_inparm/ & @@ -114,32 +112,30 @@ subroutine dynConsBiogeophys_readnl(NLFilename) dynbal_storage_residence_time = 0._r8 + ! Read in the namelist on the main task if (masterproc) then - nu_nml = getavu() - open( nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) - call find_nlgroup_name(nu_nml, 'dyn_cons_biogeophys_inparm', status=nml_error) - if (nml_error == 0) then - read(nu_nml, nml=dyn_cons_biogeophys_inparm, iostat=nml_error) - if (nml_error /= 0) then - call endrun(msg='ERROR reading dyn_cons_biogeophys_inparm namelist'//errMsg(sourcefile, __LINE__)) + open( newunit=unitn, file=trim(NLFilename), status='old', iostat=ierr ) + call shr_nl_find_group_name(unitn, nml_name, status=ierr) + if (ierr == 0) then + read(unitn, nml=dyn_cons_biogeophys_inparm, iostat=ierr) + if (ierr /= 0) then + call endrun(msg="ERROR reading "//nml_name//" namelist", file=sourcefile, line=__LINE__) end if else - call endrun(msg='ERROR finding dyn_cons_biogeophys_inparm namelist'//errMsg(sourcefile, __LINE__)) + call endrun(msg="ERROR could NOT find "//nml_name//" namelist", file=sourcefile, line=__LINE__) end if - close(nu_nml) - call relavu( nu_nml ) - endif + close( unitn ) + end if - call shr_mpi_bcast (dynbal_storage_residence_time, mpicom) + ! Broadcast namelist values to all tasks + call shr_mpi_bcast( dynbal_storage_residence_time, mpicom ) - ! Convert the residence time to a turnover rate. Note that we check for a valid - ! residence time before doing this conversion, since the conversion divides by the - ! residence time. + ! Convert the residence time to a turnover rate. if (dynbal_storage_residence_time <= 0._r8) then - write(iulog,*) 'ERROR: dynbal_storage_residence_time must be greater than 0' - write(iulog,*) 'Value given: ', dynbal_storage_residence_time - call endrun(msg='ERROR: dynbal_storage_residence_time must be greater than 0 '// & - errMsg(sourcefile, __LINE__)) + write(iulog,*) "ERROR: dynbal_storage_residence_time must be greater than 0" + write(iulog,*) "Value given: ", dynbal_storage_residence_time + call endrun(msg="ERROR: dynbal_storage_residence_time must be greater than 0", & + file=sourcefile, line=__LINE__) end if dynbal_storage_turnover_rate = 1._r8 / & (dynbal_storage_residence_time * (get_average_days_per_year() * SHR_CONST_CDAY)) @@ -147,10 +143,10 @@ subroutine dynConsBiogeophys_readnl(NLFilename) namelist_read = .true. if (masterproc) then - write(iulog,*) ' ' - write(iulog,*) 'dyn_cons_biogeophys_inparm settings:' + write(iulog,*) " " + write(iulog,*) "dyn_cons_biogeophys_inparm settings:" write(iulog,nml=dyn_cons_biogeophys_inparm) - write(iulog,*) ' ' + write(iulog,*) " " end if end subroutine dynConsBiogeophys_readnl From 4e3a4e16a871cc6cd9e923d3f1012aad2a19677a Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 8 Sep 2026 10:13:41 -0600 Subject: [PATCH 17/20] Add test back to expected fails list This DAE test was removed from the expected fails list in the last tag, but it's still failing. (It gets further than before, but still fails in case2.) --- cime_config/testdefs/ExpectedTestFails.xml | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index b677b1526f..b20dbaa8ed 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -147,6 +147,18 @@ + + + FAIL + #4202 + + + + PEND + #4202 + + + FAIL From fe49a7e395da55c3b8ca5e134a976baf63eff709 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 8 Sep 2026 10:32:44 -0600 Subject: [PATCH 18/20] Remove some tests from expected fails list These tests appear to have been passing for a while. See recent discussion in https://github.com/escomp/ctsm/issues/2619 --- cime_config/testdefs/ExpectedTestFails.xml | 24 ---------------------- 1 file changed, 24 deletions(-) diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index b20dbaa8ed..b06311cbd6 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -89,30 +89,6 @@ This failure relates to the following REP failure. - - - FAIL - #2619 - This failure relates to the following REP failure. - - - FAIL - #2619 - This failure relates to the following REP failure. - - - - - FAIL - #2619 - This failure relates to the following REP failure. - - - FAIL - #2619 - This failure relates to the following REP failure. - - FAIL From 820370a025d5cb2467e329e25339090ba0a09cbb Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 8 Sep 2026 10:34:51 -0600 Subject: [PATCH 19/20] Add a test to the expected fails list This test seems to have been failing for a while but wasn't in the expected fails list. See recent discussion in https://github.com/escomp/ctsm/issues/2619 --- cime_config/testdefs/ExpectedTestFails.xml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index b06311cbd6..72a641a901 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -83,6 +83,11 @@ + + FAIL + #2619 + This failure relates to the following REP failure. + FAIL #2619 From a5503cf7739c53e0a1d664b3d6bb978ea487d623 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 8 Sep 2026 16:02:05 -0600 Subject: [PATCH 20/20] Update ChangeLog --- doc/ChangeLog | 125 ++++++++++++++++++++++++++++++++++++++++++++++++++ doc/ChangeSum | 1 + 2 files changed, 126 insertions(+) diff --git a/doc/ChangeLog b/doc/ChangeLog index 6f3e60d07b..6504fe8948 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,129 @@ =============================================================== +Tag name: ctsm5.4.055 +Originator(s): sacks (Bill Sacks) +Date: Tue Sep 8 16:01:56 MDT 2026 +One-line Summary: Add temporal smoothing of DYNBAL conservation correction fluxes + +Purpose and description of changes +---------------------------------- + +Spread the DYNBAL water and energy fluxes over many years: Rather than releasing all of +the DYNBAL fluxes in the year in which they're generated, instead put them in storage +pools that are released with a 20-year residence time (controllable via a new namelist +parameter; the default for CLM50 is set to 0.75 to maintain something closer to the old +behavior). This reduces these fluxes, especially in the oscillatory case that is somewhat +common with a dynamic ice sheet. + +The new storage pools are reset when reset_dynbal_baselines is set to .true., since this +feels conceptually similar to the other resetting done by that flag, so it feels +appropriate to do this resetting at the same time. + +This tag also changes the sign convention for the DYNBAL fluxes: the previous sign +convention was unintuitively the opposite of the sign convention of other fluxes to +river/atm; while I'm changing these fluxes anyway, I'm also flipping the sign convention +to be more intuitive and hopefully less error-prone. + +The motivation for this change is: When subgrid areas change due to dynamic landunits / +columns, fluxes are generated to compensate for changes in gridcell-average water and +energy (QFLX_LIQ_DYNBAL, QFLX_ICE_DYNBAL, EFLX_DYNBAL). These fluxes can (and should) be +reduced via careful setting of reset_dynbal_baselines, but that still doesn't reduce these +fluxes to zero. The previous implementation seems good enough in situations where there is +a gradual change in a given grid cell's subgrid areas over many years. Unfortunately, +though, that is not always the situation when coupling to CISM, as has been seen in recent +CESM3 runs: In grid cells around the ice sheet margin, large fractions of a CTSM grid cell +can blink between mostly-vegetated/bare and mostly-glacier on a 2-year cycle, and this +oscillatory behavior can go on for many years, possibly indefinitely. This can lead to +large dynbal fluxes in each year, even when dribbled throughout the year (as was +previously done). It is problematic for the ocean to get negative dynbal fluxes, but it is +also problematic to remove these negative dynbal fluxes by spreading them globally, as is +currently done, because then regional fluxes that should average to 0 spatially and +temporally (e.g., into the Arctic ocean) end up being biased positive. This PR implements +a relatively easy partial fix, which particularly addresses the oscillatory case that +seems to be a common cause of the observed ice runoff. See +https://github.com/ESCOMP/CTSM/issues/4188 for more details. + + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- +[Remove entire section if none of the boxes are checked.] + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + + [Put an [X] in the box for any configuration with significant answer changes.] + +[X] clm6_0 + +[X] clm5_0 + +[X] ctsm5_0-nwp + +Bugs fixed +---------- + +List of CTSM issues fixed (include CTSM Issue # and description) [one per line]: +- Partially addresses #4188 - Consider one or more changes to reduce DYNBAL fluxes (dynamic landunit conservation correction fluxes) + +Notes of particular relevance for users +--------------------------------------- + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): +- New namelist variable: dynbal_storage_residence_time + +Changes to documentation: +- Updates CLM50_Tech_Note_Transient_Landcover.rst to reflect these changes + +Contributors: @billsacks + +Notes of particular relevance for developers: +--------------------------------------------- + +Changes to tests or testing: +- Fixed some expected test fails based on discussions with Erik Kluzek and Sam Levis + +Testing summary: +---------------- + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + derecho - OK + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + derecho ----- OK + izumi ------- OK + +Answer changes +-------------- + +Changes answers relative to baseline: + + Summarize any changes to answers, i.e., + - what code configurations: + - Transient cases (Hist, SSP) + - Cases with an evolving ice sheet (active CISM) + - FATES cases + - Cases that start mid-year from a restart file that had non-zero values in the old dribblers + - what platforms/compilers: all + - nature of change (roundoff; larger than roundoff/same climate; new climate): + - larger than roundoff + - same vs. new climate depends on the configuration and what variables are of interest + +Non-BGC I compset cases have diffs only in limited variables: This PR only directly +affects sensible heat, liquid and ice runoff fluxes, so the only changes are in these +fluxes and related ROF variables. However, some BGC cases have more extensive diffs, I +think due to the feedback of ROF on CH4. For a few long BGC tests with extensive diffs, I +verified that, by changing the compset to one using SROF, the set of variables with diffs +became small (similar to the set of variables with diffs in non-BGC cases). + + +Other details +------------- + +Pull Requests that document the changes (include PR ids): https://github.com/ESCOMP/CTSM/pull/4189 + +=============================================================== +=============================================================== Tag name: ctsm5.4.054 Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) Date: Mon Aug 31 05:35:20 AM MDT 2026 diff --git a/doc/ChangeSum b/doc/ChangeSum index e9c62924ba..520580bd66 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.4.055 sacks 09/08/2026 Add temporal smoothing of DYNBAL conservation correction fluxes ctsm5.4.054 erik 08/31/2026 Update submodules to the latest needed for cesm3_0_alpha10a ctsm5.4.053 oleson 08/18/2026 Resolve problems with running PLUMBER2 sites ctsm5.4.052 slevis 08/17/2026 Merge b4b-dev to master