Skip to content

Commit

Permalink
Merge pull request #3 from billsacks/dynlakes_avoid_tws_changes
Browse files Browse the repository at this point in the history 
Dynlakes: fix some subtle issues

### Description of changes

Fix a variety of subtle issues with dynamic lakes - particularly the accounting of total water and energy.

### Specific notes

This branch contains the following commits:

- a9fa875: This is needed to avoid counting lake water in the begwb and endwb terms, which is needed because these are used to calculate gridcell total water store (TWS), which in turn influences the methane code. Because the methane code was tuned around old values of TWS, changing TWS would lead to unintentional – and potentially large – changes in methane terms. Eventually we'd like to remove methane's dependence on TWS, but for now this workaround is needed to avoid changing behavior too much. See ESCOMP#659 (comment) for more details.

- 52105c4: This minor fix is needed for the sake of water tracers / water isotopes. It shouldn't have any impact outside of that (because the tracer_ratio of bulk water is 1)

- de3e12c and acf0984: This one is especially subtle; it is needed for backwards compatibility with old restart files. The main changes are in de3e12c; acf0984 is just a minor tweak on top of that. The problem is that, on existing initial conditions files, there can be already-existing DYNBAL_BASELINE variables (for LIQUID, ICE & HEAT). But these pre-existing variables will have baseline values of 0 for lake. Before this commit, when you started up from an old initial conditions file, the code would use these 0 values for lake baselines (because baselines are only reset if the user explicitly asks them to be reset with a namelist flag). This commit adds some code to detect if the initial conditions file is old, and if so, recomputes dynbal baselines for lake using the new definition. Note that some even older initial conditions files didn't have the DYNBAL_BASELINE variables at all; those would have been okay before this change: the problem is with initial conditions files that are somewhat but not very old - so have DYNBAL_BASELINE variables on them that use the old definition (where lake baseline values were 0).

- 8088c3c: Minor fix for a pre-existing issue

- a31875d: I'm not sure if this is actually needed, but I thought it would be good to group together the lake water content and the roughly equal-but-opposite baselines, so that these can cancel to near zero before adding the smaller terms. In principle, this should help maintain precision in these smaller terms. I thought this might help resolve some of the larger-than-expected answer changes I was seeing in testing, but I don't think it actually does... but I still thought this would be good to keep in place. **I have double and triple checked these changes, but it would be good to have an extra set of eyes on them to make sure I did this reordering correctly. In particular, I think there were some subtleties about when a term should accumulate on top of an existing value vs. setting the initial value of a variable.**
  • Loading branch information
@billsacks
billsacks authored Sep 28, 2020
2 parents 8ff6ca3 + a31875d commit 2cb54b5
Show file tree
Hide file tree
Showing 6 changed files with 357 additions and 185 deletions.
2 changes: 1 addition & 1 deletion src/biogeophys/BalanceCheckMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@ subroutine BeginWaterBalanceSingle(bounds, &

call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
waterstate_inst, lakestate_inst, &
subtract_dynbal_baselines = .false., &
add_lake_water_and_subtract_dynbal_baselines = .false., &
water_mass = begwb(bounds%begc:bounds%endc))

call waterstate_inst%CalculateTotalH2osno(bounds, num_nolakec, filter_nolakec, &
Expand Down
2 changes: 1 addition & 1 deletion src/biogeophys/LakeHydrologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -650,7 +650,7 @@ subroutine LakeHydrology(bounds, &

call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
b_waterstate_inst, lakestate_inst, &
subtract_dynbal_baselines = .false., &
add_lake_water_and_subtract_dynbal_baselines = .false., &
water_mass = endwb(bounds%begc:bounds%endc))

do j = 1, nlevgrnd
Expand Down
144 changes: 89 additions & 55 deletions src/biogeophys/TotalWaterAndHeatMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@ end subroutine ComputeWaterMassNonLake
!-----------------------------------------------------------------------
subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
waterstate_inst, lakestate_inst, &
subtract_dynbal_baselines, &
add_lake_water_and_subtract_dynbal_baselines, &
water_mass)
!
! !DESCRIPTION:
Expand All @@ -159,9 +159,12 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
class(waterstate_type) , intent(in) :: waterstate_inst
type(lakestate_type) , intent(in) :: lakestate_inst

! Whether to (1) add lake water/ice to total accounting, and (2) subtract
! dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
!
! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658
! is resolved, remove this argument, always assuming it's true.
logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
logical, intent(in) :: add_lake_water_and_subtract_dynbal_baselines

real(r8) , intent(inout) :: water_mass( bounds%begc: ) ! computed water mass (kg m-2)
!
Expand All @@ -181,7 +184,7 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
filter_lakec = filter_lakec, &
waterstate_inst = waterstate_inst, &
lakestate_inst = lakestate_inst, &
subtract_dynbal_baselines = subtract_dynbal_baselines, &
add_lake_water_and_subtract_dynbal_baselines = add_lake_water_and_subtract_dynbal_baselines, &
liquid_mass = liquid_mass(bounds%begc:bounds%endc), &
ice_mass = ice_mass(bounds%begc:bounds%endc))

Expand Down Expand Up @@ -387,8 +390,8 @@ end subroutine AccumulateSoilLiqIceMassNonLake

!-----------------------------------------------------------------------
subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
waterstate_inst, lakestate_inst,&
subtract_dynbal_baselines, &
waterstate_inst, lakestate_inst, &
add_lake_water_and_subtract_dynbal_baselines, &
liquid_mass, ice_mass)
!
! !DESCRIPTION:
Expand All @@ -407,10 +410,12 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
class(waterstate_type), intent(in) :: waterstate_inst
type(lakestate_type) , intent(in) :: lakestate_inst


! Whether to (1) add lake water/ice to total accounting, and (2) subtract
! dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
!
! BUG(wjs, 2019-03-12, ESCOMP/ctsm#659) When https://github.com/ESCOMP/CTSM/issues/658
! is resolved, remove this argument, always assuming it's true.
logical, intent(in) :: subtract_dynbal_baselines ! whether to subtract dynbal_baseline_liq and dynbal_baseline_ice from liquid_mass and ice_mass
logical, intent(in) :: add_lake_water_and_subtract_dynbal_baselines

real(r8) , intent(inout) :: liquid_mass( bounds%begc: ) ! computed liquid water mass (kg m-2)
real(r8) , intent(inout) :: ice_mass( bounds%begc: ) ! computed ice mass (kg m-2)
Expand Down Expand Up @@ -439,6 +444,34 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
ice_mass(c) = 0._r8
end do

! ------------------------------------------------------------------------
! Start with some large terms that will often cancel (the negative baselines and the
! lake water content): In order to maintain precision in the other terms, it can help if
! we deal with these large, often-canceling terms first. (If we accumulated some
! small terms, then added a big term and then subtracted a big term, we would have
! lost some precision in the small terms.)
! ------------------------------------------------------------------------

if (add_lake_water_and_subtract_dynbal_baselines) then
! Subtract baselines set in initialization
do fc = 1, num_lakec
c = filter_lakec(fc)
liquid_mass(c) = liquid_mass(c) - dynbal_baseline_liq(c)
ice_mass(c) = ice_mass(c) - dynbal_baseline_ice(c)
end do

! Lake water content
call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, &
lakestate_inst, &
tracer_ratio = waterstate_inst%info%get_ratio(), &
liquid_mass = liquid_mass(bounds%begc:bounds%endc), &
ice_mass = ice_mass(bounds%begc:bounds%endc))
end if

! ------------------------------------------------------------------------
! Now add some other terms
! ------------------------------------------------------------------------

! Snow water content
do fc = 1, num_lakec
c = filter_lakec(fc)
Expand All @@ -449,11 +482,7 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j)
end do
end do

! Lake water content
call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, &
lakestate_inst, liquid_mass, ice_mass)


! Soil water content of the soil under the lake
do j = 1, nlevgrnd
do fc = 1, num_lakec
Expand All @@ -463,22 +492,13 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
end do
end do

if (subtract_dynbal_baselines) then
! Subtract baselines set in initialization
do fc = 1, num_lakec
c = filter_lakec(fc)
liquid_mass(c) = liquid_mass(c) - dynbal_baseline_liq(c)
ice_mass(c) = ice_mass(c) - dynbal_baseline_ice(c)
end do
end if

end associate

end subroutine ComputeLiqIceMassLake

!-----------------------------------------------------------------------
subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, &
lakestate_inst, liquid_mass, ice_mass)
lakestate_inst, tracer_ratio, liquid_mass, ice_mass)
!
! !DESCRIPTION:
! Accumulate lake water mass of lake columns, separated into liquid and ice.
Expand All @@ -493,6 +513,7 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, &
integer , intent(in) :: num_c ! number of column points in column filter (should not include lake points; can be a subset of the no-lake filter)
integer , intent(in) :: filter_c(:) ! column filter (should not include lake points; can be a subset of the no-lake filter)
type(lakestate_type) , intent(in) :: lakestate_inst
real(r8) , intent(in) :: tracer_ratio ! for water tracers, standard ratio of this tracer to bulk water (1 for bulk water)
real(r8) , intent(inout) :: liquid_mass( bounds%begc: ) ! accumulated liquid mass (kg m-2)
real(r8) , intent(inout) :: ice_mass( bounds%begc: ) ! accumulated ice mass (kg m-2)
!
Expand All @@ -516,12 +537,14 @@ subroutine AccumulateLiqIceMassLake(bounds, num_c, filter_c, &
do j = 1, nlevlak
do fc = 1, num_c
c = filter_c(fc)
! calculate lake liq and ice content per lake layer first
h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j))
! Lake water volume isn't tracked explicitly, so we calculate it from lake
! depth. Because it isn't tracked explicitly, we also don't have any water
! tracer information, so we assume a fixed, standard tracer ratio.
h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j)) * tracer_ratio

! use water density rather than ice density because lake layer depths are not
! adjusted when the layer freezes
h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j)
h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) * tracer_ratio

liquid_mass(c) = liquid_mass(c) + h2olak_liq
ice_mass(c) = ice_mass(c) + h2olak_ice
Expand Down Expand Up @@ -936,6 +959,42 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, &
latent_heat_liquid(c) = 0._r8
end do

! ------------------------------------------------------------------------
! Start with some large terms that will often cancel (the negative baselines and the
! lake water content): In order to maintain precision in the other terms, it can help if
! we deal with these large, often-canceling terms first. (If we accumulated some
! small terms, then added a big term and then subtracted a big term, we would have
! lost some precision in the small terms.)
! ------------------------------------------------------------------------

! Subtract baselines set in initialization
!
! NOTE(wjs, 2019-03-01) I haven't given enough thought to how (if at all) we should
! correct for heat_liquid and cv_liquid, which are used to determine the weighted
! average liquid water temperature. For example, if we're subtracting out a baseline
! water amount because a particular water state is fictitious, we probably shouldn't
! include that particular state when determining the weighted average temperature of
! liquid water. And conversely, if we're adding a state via these baselines, should
! we also add some water temperature of that state? The tricky thing here is what to
! do when we end up subtracting water due to the baselines, so for now I'm simply not
! trying to account for the temperature of these baselines at all. This amounts to
! assuming that the baselines that we add / subtract are at the average temperature
! of the real liquid water in the column.
do fc = 1, num_lakec
c = filter_lakec(fc)
heat(c) = -dynbal_baseline_heat(c)
end do

! Lake water heat content
! Note that we do NOT accumulate heat_liquid and cv_liquid in this call. See the
! comments at the top of AccumulateHeatLake for rationale.
call AccumulateHeatLake(bounds, num_lakec, filter_lakec, temperature_inst, lakestate_inst, &
heat)

! ------------------------------------------------------------------------
! Now add some other terms
! ------------------------------------------------------------------------

! Snow heat content
do fc = 1, num_lakec
c = filter_lakec(fc)
Expand Down Expand Up @@ -976,41 +1035,16 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, &

do fc = 1, num_lakec
c = filter_lakec(fc)
heat(c) = heat_dry_mass(c) + heat_ice(c) + heat_liquid(c) + latent_heat_liquid(c)
heat(c) = heat(c) + heat_dry_mass(c) + heat_ice(c) + heat_liquid(c) + latent_heat_liquid(c)
end do

! Lake water heat content
! Note that we do NOT accumulate heat_liquid and cv_liquid in this call. See the
! comments at the top of AccumulateHeatLake for rationale.
call AccumulateHeatLake(bounds, num_lakec, filter_lakec, temperature_inst, lakestate_inst, &
heat)

! Subtract baselines set in initialization
!
!
! NOTE(wjs, 2019-03-01) I haven't given enough thought to how (if at all) we should
! correct for heat_liquid and cv_liquid, which are used to determine the weighted
! average liquid water temperature. For example, if we're subtracting out a baseline
! water amount because a particular water state is fictitious, we probably shouldn't
! include that particular state when determining the weighted average temperature of
! liquid water. And conversely, if we're adding a state via these baselines, should
! we also add some water temperature of that state? The tricky thing here is what to
! do when we end up subtracting water due to the baselines, so for now I'm simply not
! trying to account for the temperature of these baselines at all. This amounts to
! assuming that the baselines that we add / subtract are at the average temperature
! of the real liquid water in the column.
do fc = 1, num_lakec
c = filter_lakec(fc)
heat(c) = heat(c) - dynbal_baseline_heat(c)
end do

end associate

end subroutine ComputeHeatLake

!-----------------------------------------------------------------------
subroutine AccumulateHeatLake(bounds, num_c, filter_c, &
temperature_inst, lakestate_inst, &
temperature_inst, lakestate_inst, &
heat)
!
! !DESCRIPTION:
Expand Down Expand Up @@ -1097,11 +1131,11 @@ subroutine AccumulateHeatLake(bounds, num_c, filter_c, &
end do
end do

! add ice heat and liquid heat together
! add ice heat and liquid heat together
do fc = 1, num_c
c = filter_c(fc)
heat(c) = heat(c) + lake_heat_ice(c) + &
lake_heat_liquid(c) + lake_latent_heat_liquid(c)
heat(c) = heat(c) + (lake_heat_ice(c) + &
lake_heat_liquid(c) + lake_latent_heat_liquid(c))
end do

end associate
Expand Down
Loading

0 comments on commit 2cb54b5

Please sign in to comment.