Skip to content

Commit

Permalink
Reset lake dynbal baselines if using an old initial conditions file
Browse files Browse the repository at this point in the history 
We have changed the definition of total column water for lake columns,
so the baseline values for lakes are incorrect on old initial conditions
files. This commit adds some code to check if we're using an old initial
conditions file, and if so, resets the dynbal baseline values for lakes
to use the new definition.
  • Loading branch information
@billsacks
billsacks committed Sep 3, 2020
1 parent 52105c4 commit de3e12c
Show file tree
Hide file tree
Showing 3 changed files with 269 additions and 126 deletions.
251 changes: 152 additions & 99 deletions src/dyn_subgrid/dynConsBiogeophysMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ module dynConsBiogeophysMod
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_log_mod , only : errMsg => shr_log_errMsg
use decompMod , only : bounds_type
use spmdMod , only : masterproc
use clm_varctl , only : iulog
use UrbanParamsType , only : urbanparams_type
use EnergyFluxType , only : energyflux_type
use SoilHydrologyType , only : soilhydrology_type
Expand Down Expand Up @@ -66,7 +68,8 @@ module dynConsBiogeophysMod
!-----------------------------------------------------------------------
subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
num_lakec, filter_lakec, &
urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst)
urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, &
reset_all_baselines, reset_lake_baselines, print_resets)
!
! !DESCRIPTION:
! Set start-of-run baseline values for heat and water content in some columns.
Expand All @@ -79,13 +82,18 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
! represented). These corrections will typically reduce the fictitious dynbal
! conservation fluxes.
!
! At a minimum, this should be called during cold start initialization, to initialize
! the baseline values based on cold start states. In addition, this can be called when
! starting a new run from existing initial conditions, if the user desires. This
! optional resetting can further reduce 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.)
! At a minimum, this should be called during cold start initialization with
! reset_all_baselines set to true, to initialize the baseline values based on cold
! start states. This should also be called after reading initial conditions, but the
! various reset_* flags should be set appropriately; setting all of these flags to
! false will result in no baselines being reset in this call.

! Setting reset_all_baselines can be done when starting a new run from existing
! initial conditions if the user desires. This optional resetting can further reduce
! 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.)
! Other reset_* flags are described below.
!
! !ARGUMENTS:
type(bounds_type), intent(in) :: bounds
Expand All @@ -98,12 +106,28 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
integer, intent(in) :: num_lakec ! number of points in filter_lakec
integer, intent(in) :: filter_lakec(:) ! filter for lake columns


type(urbanparams_type), intent(in) :: urbanparams_inst
type(soilstate_type), intent(in) :: soilstate_inst
type(lakestate_type), intent(in) :: lakestate_inst
type(water_type), intent(inout) :: water_inst
type(temperature_type), intent(inout) :: temperature_inst

! Whether to reset baselines for all columns
logical, intent(in) :: reset_all_baselines

! Whether to reset baselines for lake columns. Ignored if reset_all_baselines is
! true.
!
! BACKWARDS_COMPATIBILITY(wjs, 2020-09-02) This is needed when reading old initial
! conditions files created before https://github.com/ESCOMP/CTSM/issues/1140 was
! resolved via https://github.com/ESCOMP/CTSM/pull/1109: The definition of total
! column water content has been changed for lakes, so we need to reset baseline values
! for lakes on older initial conditions files.
logical, intent(in) :: reset_lake_baselines

! Whether to print information about the reset flags
logical, intent(in) :: print_resets

!
! !LOCAL VARIABLES:
integer :: i
Expand All @@ -112,6 +136,18 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &
character(len=*), parameter :: subname = 'dyn_hwcontent_set_baselines'
!-----------------------------------------------------------------------

if (masterproc .and. print_resets) then
if (reset_all_baselines) then
write(iulog,*) ' '
write(iulog,*) 'Resetting dynbal baselines for all columns'
write(iulog,*) ' '
else if (reset_lake_baselines) then
write(iulog,*) ' '
write(iulog,*) 'Resetting dynbal baselines for lake columns'
write(iulog,*) ' '
end if
end if

! Note that we include inactive points in the following. This could be important if
! an inactive point later becomes active, so that we have an appropriate baseline
! value for that point.
Expand All @@ -125,22 +161,26 @@ subroutine dyn_hwcontent_set_baselines(bounds, num_icemecc, filter_icemecc, &

call dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
bulk_or_tracer%waterstate_inst, lakestate_inst)
bulk_or_tracer%waterstate_inst, lakestate_inst, &
reset_all_baselines = reset_all_baselines, &
reset_lake_baselines = reset_lake_baselines)
end associate
end do

call dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
urbanparams_inst, soilstate_inst, lakestate_inst, water_inst%waterstatebulk_inst, &
temperature_inst)

temperature_inst, &
reset_all_baselines = reset_all_baselines, &
reset_lake_baselines = reset_lake_baselines)

end subroutine dyn_hwcontent_set_baselines

!-----------------------------------------------------------------------
subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
waterstate_inst, lakestate_inst)
waterstate_inst, lakestate_inst, &
reset_all_baselines, reset_lake_baselines)
!
! !DESCRIPTION:
! Set start-of-run baseline values for water content, for a single water tracer or
Expand All @@ -155,9 +195,13 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
integer, intent(in) :: filter_icemecc(:) ! filter for icemec (i.e., glacier) columns
integer, intent(in) :: num_lakec ! number of points in filter_lakec
integer, intent(in) :: filter_lakec(:) ! filter for lake columns
type(lakestate_type), intent(in) :: lakestate_inst

type(lakestate_type), intent(in) :: lakestate_inst
class(waterstate_type), intent(inout) :: waterstate_inst

! See dyn_hwcontent_set_baselines for documentation of these arguments
logical, intent(in) :: reset_all_baselines
logical, intent(in) :: reset_lake_baselines
!
! !LOCAL VARIABLES:
integer :: c, fc ! indices
Expand All @@ -174,47 +218,50 @@ subroutine dyn_water_content_set_baselines(bounds, natveg_and_glc_filterc, &
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)
)

soil_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
soil_ice_mass_col (bounds%begc:bounds%endc) = 0._r8
lake_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
lake_ice_mass_col (bounds%begc:bounds%endc) = 0._r8


call AccumulateSoilLiqIceMassNonLake(bounds, &
natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
waterstate_inst, &
liquid_mass = soil_liquid_mass_col(bounds%begc:bounds%endc), &
ice_mass = soil_ice_mass_col(bounds%begc:bounds%endc))

! For glacier columns, the liquid and ice in the "soil" (i.e., in the glacial ice) is
! a virtual pool. So we'll subtract this amount when determining the dynbal
! fluxes. (Note that a positive value in these baseline variables indicates an extra
! stock that will be subtracted later.) But glacier columns do not represent the soil
! under the glacial ice. Let's assume that the soil state under each glacier column is
! the same as the soil state in the natural vegetation landunit on that grid cell. We
! subtract this from the dynbal baseline variables to indicate a missing stock.
call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
vals_col = soil_liquid_mass_col(bounds%begc:bounds%endc), &
baselines_col = dynbal_baseline_liq(bounds%begc:bounds%endc))
call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
vals_col = soil_ice_mass_col(bounds%begc:bounds%endc), &
baselines_col = dynbal_baseline_ice(bounds%begc:bounds%endc))


! set baselines for lake columns

! Calculate the total water volume of the lake column
call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, &
tracer_ratio = waterstate_inst%info%get_ratio(), &
liquid_mass = lake_liquid_mass_col(bounds%begc:bounds%endc), &
ice_mass = lake_ice_mass_col(bounds%begc:bounds%endc))

do fc = 1, num_lakec
c = filter_lakec(fc)
dynbal_baseline_liq(c) = lake_liquid_mass_col(c)
dynbal_baseline_ice(c) = lake_ice_mass_col(c)
end do
if (reset_all_baselines) then
soil_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
soil_ice_mass_col (bounds%begc:bounds%endc) = 0._r8

call AccumulateSoilLiqIceMassNonLake(bounds, &
natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
waterstate_inst, &
liquid_mass = soil_liquid_mass_col(bounds%begc:bounds%endc), &
ice_mass = soil_ice_mass_col(bounds%begc:bounds%endc))

! For glacier columns, the liquid and ice in the "soil" (i.e., in the glacial ice) is
! a virtual pool. So we'll subtract this amount when determining the dynbal
! fluxes. (Note that a positive value in these baseline variables indicates an extra
! stock that will be subtracted later.) But glacier columns do not represent the soil
! under the glacial ice. Let's assume that the soil state under each glacier column is
! the same as the soil state in the natural vegetation landunit on that grid cell. We
! subtract this from the dynbal baseline variables to indicate a missing stock.
call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
vals_col = soil_liquid_mass_col(bounds%begc:bounds%endc), &
baselines_col = dynbal_baseline_liq(bounds%begc:bounds%endc))
call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
vals_col = soil_ice_mass_col(bounds%begc:bounds%endc), &
baselines_col = dynbal_baseline_ice(bounds%begc:bounds%endc))
end if

if (reset_all_baselines .or. reset_lake_baselines) then
! set baselines for lake columns

lake_liquid_mass_col(bounds%begc:bounds%endc) = 0._r8
lake_ice_mass_col (bounds%begc:bounds%endc) = 0._r8

! Calculate the total water volume of the lake column
call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, lakestate_inst, &
tracer_ratio = waterstate_inst%info%get_ratio(), &
liquid_mass = lake_liquid_mass_col(bounds%begc:bounds%endc), &
ice_mass = lake_ice_mass_col(bounds%begc:bounds%endc))

do fc = 1, num_lakec
c = filter_lakec(fc)
dynbal_baseline_liq(c) = lake_liquid_mass_col(c)
dynbal_baseline_ice(c) = lake_ice_mass_col(c)
end do
end if

end associate

end subroutine dyn_water_content_set_baselines
Expand All @@ -223,7 +270,8 @@ end subroutine dyn_water_content_set_baselines
subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
num_icemecc, filter_icemecc, num_lakec, filter_lakec, &
urbanparams_inst, soilstate_inst, lakestate_inst, waterstatebulk_inst, &
temperature_inst)
temperature_inst, &
reset_all_baselines, reset_lake_baselines)
!
! !DESCRIPTION:
! Set start-of-run baseline values for heat content.
Expand All @@ -243,6 +291,10 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
type(lakestate_type) , intent(in) :: lakestate_inst
type(waterstatebulk_type), intent(in) :: waterstatebulk_inst
type(temperature_type), intent(inout) :: temperature_inst

! See dyn_hwcontent_set_baselines for documentation of these arguments
logical, intent(in) :: reset_all_baselines
logical, intent(in) :: reset_lake_baselines
!
! !LOCAL VARIABLES:
integer :: c, fc ! indices
Expand All @@ -258,51 +310,52 @@ subroutine dyn_heat_content_set_baselines(bounds, natveg_and_glc_filterc, &
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)
)

soil_heat_col(bounds%begc:bounds%endc) = 0._r8
soil_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8
soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8

lake_heat_col(bounds%begc:bounds%endc) = 0._r8

call AccumulateSoilHeatNonLake(bounds, &
natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
urbanparams_inst, soilstate_inst, temperature_inst, waterstatebulk_inst, &
heat = soil_heat_col(bounds%begc:bounds%endc), &
heat_liquid = soil_heat_liquid_col(bounds%begc:bounds%endc), &
cv_liquid = soil_cv_liquid_col(bounds%begc:bounds%endc))


! See comments in dyn_water_content_set_baselines for rationale for these glacier
! baselines. Even though the heat in glacier ice can interact with the rest of the
! system (e.g., giving off heat to the atmosphere or receiving heat from the
! atmosphere), it is still a virtual state in the sense that the glacier ice column
! is virtual. And, as for water, we subtract the heat of the soil in the associated
! natural vegetated landunit to account for the fact that we don't explicitly model
! the soil under glacial ice.
!
! Aside from these considerations of what is virtual vs. real, another rationale
! driving the use of these baselines is the desire to minimize the dynbal fluxes. For
! the sake of conservation, it seems okay to pick any fixed baseline when summing the
! energy (or water) content of a given column (as long as that baseline doesn't
! change over time). By using the baselines computed here, we reduce the dynbal
! fluxes to more reasonable values.
call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
vals_col = soil_heat_col(bounds%begc:bounds%endc), &
baselines_col = dynbal_baseline_heat(bounds%begc:bounds%endc))


! Set baselines for lake columns
call AccumulateHeatLake(bounds, num_lakec, filter_lakec, &
temperature_inst, lakestate_inst, &
heat = lake_heat_col)

do fc = 1, num_lakec
c = filter_lakec(fc)
if (reset_all_baselines) then
soil_heat_col(bounds%begc:bounds%endc) = 0._r8
soil_heat_liquid_col(bounds%begc:bounds%endc) = 0._r8
soil_cv_liquid_col(bounds%begc:bounds%endc) = 0._r8

call AccumulateSoilHeatNonLake(bounds, &
natveg_and_glc_filterc%num, natveg_and_glc_filterc%indices, &
urbanparams_inst, soilstate_inst, temperature_inst, waterstatebulk_inst, &
heat = soil_heat_col(bounds%begc:bounds%endc), &
heat_liquid = soil_heat_liquid_col(bounds%begc:bounds%endc), &
cv_liquid = soil_cv_liquid_col(bounds%begc:bounds%endc))

! See comments in dyn_water_content_set_baselines for rationale for these glacier
! baselines. Even though the heat in glacier ice can interact with the rest of the
! system (e.g., giving off heat to the atmosphere or receiving heat from the
! atmosphere), it is still a virtual state in the sense that the glacier ice column
! is virtual. And, as for water, we subtract the heat of the soil in the associated
! natural vegetated landunit to account for the fact that we don't explicitly model
! the soil under glacial ice.
!
! Aside from these considerations of what is virtual vs. real, another rationale
! driving the use of these baselines is the desire to minimize the dynbal fluxes. For
! the sake of conservation, it seems okay to pick any fixed baseline when summing the
! energy (or water) content of a given column (as long as that baseline doesn't
! change over time). By using the baselines computed here, we reduce the dynbal
! fluxes to more reasonable values.
call set_glacier_baselines(bounds, num_icemecc, filter_icemecc, &
vals_col = soil_heat_col(bounds%begc:bounds%endc), &
baselines_col = dynbal_baseline_heat(bounds%begc:bounds%endc))
end if

if (reset_all_baselines .or. reset_lake_baselines) then
! Set baselines for lake columns

lake_heat_col(bounds%begc:bounds%endc) = 0._r8

call AccumulateHeatLake(bounds, num_lakec, filter_lakec, &
temperature_inst, lakestate_inst, &
heat = lake_heat_col)

do fc = 1, num_lakec
c = filter_lakec(fc)
dynbal_baseline_heat(c) = lake_heat_col(c)
end do
end if

dynbal_baseline_heat(c) = lake_heat_col(c)

end do

end associate

end subroutine dyn_heat_content_set_baselines
Expand Down
Loading

0 comments on commit de3e12c

Please sign in to comment.