Skip to content

Commit

Permalink
Merge tag 'ctsm5.1.dev122'
Browse files Browse the repository at this point in the history 
Rework handling of evaporation constraint in SoilFluxes

Occasionally, h2osoi_ice was going significantly negative in
UpdateState_TopLayerFluxes - see
ESCOMP#1979. As noted in that issue, this
seems to be due to h2osoi_ice having a very different magnitude from
h2osoi_liq, leading to greater-than-roundoff-level differences from zero
final state in a relative sense (i.e., relative to the magnitude of
h2osoi_ice) - I think because of the appearance of the sum (h2osoi_ice +
h2osoi_liq) in the equations that limit fluxes.

To try to deal with this, I have reworked the handling of the
evaporation constraint to directly limit both the liqevap and solidevap,
so that both of them should result in the equivalent liq or ice states
going to 0 within roundoff.

To do that, I needed to move the partitioning of the total flux into
liquid and solid to earlier in the subroutine and then recalculate those
partitioning fluxes in conditions where we're applying an evaporation
constraint.

Note that I applied a max of 0 to the new fluxes because many initial
conditions files have roundoff-level negative H2OSOI_LIQ, so without
this limit, we were getting roundoff-level negative fluxes.

Resolves ESCOMP#1979
  • Loading branch information
@ekluzek
ekluzek committed Apr 27, 2023
2 parents 95c5c87 + cdb6187 commit 798c29c
Show file tree
Hide file tree
Showing 3 changed files with 253 additions and 73 deletions.
166 changes: 166 additions & 0 deletions doc/ChangeLog
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,170 @@
===============================================================
Tag name: ctsm5.1.dev122
Originator(s): sacks (Bill Sacks)
Date: Sun Apr 23 19:36:37 MDT 2023
One-line Summary: Rework handling of evaporation constraint in SoilFluxes

Purpose and description of changes
----------------------------------

Occasionally, h2osoi_ice was going significantly negative in
UpdateState_TopLayerFluxes - see
https://github.com/ESCOMP/CTSM/issues/1979. As noted in that issue, this
seems to be due to h2osoi_ice having a very different magnitude from
h2osoi_liq, leading to greater-than-roundoff-level differences from zero
final state in a relative sense (i.e., relative to the magnitude of
h2osoi_ice) - I think because of the appearance of the sum (h2osoi_ice +
h2osoi_liq) in the equations that limit fluxes.

To try to deal with this, I have reworked the handling of the
evaporation constraint to directly limit both the liqevap and solidevap,
so that both of them should result in the equivalent liq or ice states
going to 0 within roundoff.

To do that, I needed to move the partitioning of the total flux into
liquid and solid to earlier in the subroutine and then recalculate those
partitioning fluxes in conditions where we're applying an evaporation
constraint.

Note that I applied a max of 0 to the new fluxes because many initial
conditions files have roundoff-level negative H2OSOI_LIQ, so without
this limit, we were getting roundoff-level negative fluxes.

Significant changes to scientifically-supported configurations
--------------------------------------------------------------

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.]

[ ] clm5_1

[ ] clm5_0

[ ] ctsm5_0-nwp

[ ] clm4_5


Bugs fixed or introduced
------------------------
CTSM issues fixed (include CTSM Issue #):
- Resolves ESCOMP/CTSM#1979 (Need some changes to avoid negative h2osoi_ice in UpdateState_TopLayerFluxes)


Testing summary:
----------------

regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing):

cheyenne ---- OK
izumi ------- OK

Tests passed, some baseline differences as expected.

Answer changes
--------------

Changes answers relative to baseline: YES

Summarize any changes to answers, i.e.,
- what code configurations: potentially all
- what platforms/compilers: potentially all
- nature of change (roundoff; larger than roundoff/same climate; new climate):
roundoff

Differences were only observed in a few tests:
- ERP_P36x2_Ld30.f45_f45_mg37.I2000Clm51FatesSpCruRsGs.cheyenne_intel.clm-FatesColdSatPhen
- ERI_D_Ld9_P48x1.T31_g37.I2000Clm50Sp.izumi_nag.clm-reduceOutput
- SMS_D_Ln9_P36x3.f19_g17.IHistClm50Sp.cheyenne_intel.clm-waccmx_offline
- SMS_D_Ln9_P36x3_Vmct.f19_g17.IHistClm50Sp.cheyenne_intel.clm-waccmx_offline

If bitwise differences were observed, how did you show they were no worse
than roundoff?

Only two tests had greater-than-roundoff-level differences in the
cprnc output:
SMS_D_Ln9_P36x3.f19_g17.IHistClm50Sp.cheyenne_intel.clm-waccmx_offline
and the mct equivalent,
SMS_D_Ln9_P36x3_Vmct.f19_g17.IHistClm50Sp.cheyenne_intel.clm-waccmx_offline.
To verify that differences were fundamentally no greater than
roundoff-level, I introduced temporary code; this minimal diff
ended up being enough to give just roundoff-level differences from baseline:

diff --git a/src/biogeophys/SoilFluxesMod.F90 b/src/biogeophys/SoilFluxesMod.F90
index c316d30fe..6a958c0ee 100644
--- a/src/biogeophys/SoilFluxesMod.F90
+++ b/src/biogeophys/SoilFluxesMod.F90
@@ -45,7 +45,7 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
! Update surface fluxes based on the new ground temperature
!
! !USES:
- use clm_time_manager , only : get_step_size_real
+ use clm_time_manager , only : get_step_size_real, get_nstep
use clm_varcon , only : hvap, cpair, grav, vkc, tfrz, sb
use landunit_varcon , only : istsoil, istcrop
use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv
@@ -79,7 +79,9 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
real(r8) :: t_grnd0(bounds%begc:bounds%endc) ! t_grnd of previous time step
real(r8) :: lw_grnd
real(r8) :: evaporation_limit ! top layer moisture available for evaporation
- real(r8) :: evaporation_demand ! evaporative demand
+ real(r8) :: evaporation_demand ! evaporative demand
+ real(r8) :: qflx_liqevap_orig
+ real(r8) :: qflx_solidevap_orig
!-----------------------------------------------------------------------

associate( &
@@ -291,6 +293,7 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
qflx_evap_soi(p) = qflx_evap_soi(p) - frac_sno_eff(c)*(evaporation_demand - evaporation_limit)
qflx_liqevap_from_top_layer(p) = max(h2osoi_liq(c,j)/(frac_sno_eff(c)*dtime), 0._r8)
qflx_solidevap_from_top_layer(p) = max(h2osoi_ice(c,j)/(frac_sno_eff(c)*dtime), 0._r8)
+
! conserve total energy flux
eflx_sh_grnd(p) = eflx_sh_grnd(p) + frac_sno_eff(c)*(evaporation_demand - evaporation_limit)*htvp(c)
endif
@@ -307,6 +310,24 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
qflx_ev_snow(p) = qflx_evap_soi(p)
qflx_liqevap_from_top_layer(p) = max(h2osoi_liq(c,j)/dtime, 0._r8)
qflx_solidevap_from_top_layer(p) = max(h2osoi_ice(c,j)/dtime, 0._r8)
+
+ if (h2osoi_liq(c,j) + h2osoi_ice(c,j) > 0._r8) then
+ qflx_liqevap_orig = max(qflx_evap_soi(p)*(h2osoi_liq(c,j)/ &
+ (h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
+ else
+ qflx_liqevap_orig = 0._r8
+ end if
+ qflx_solidevap_orig = qflx_evap_soi(p) - qflx_liqevap_orig
+ if (qflx_solidevap_from_top_layer(p) == 0._r8 .and. &
+ qflx_solidevap_orig < 0._r8 .and. &
+ qflx_solidevap_orig > -1.e-16_r8) then
+ write(iulog,'(a, i0, 1x, i0, 1x, 5e24.17)') &
+ 'WJS adj urb: solid orig le 0, new 0: nstep, p, orig, new, qflx_evap_soi, h2osoi_liq, h2osoi_ice = ', &
+ get_nstep(), p, qflx_solidevap_orig, qflx_solidevap_from_top_layer(p), &
+ qflx_evap_soi(p), h2osoi_liq(c,j), h2osoi_ice(c,j)
+ qflx_solidevap_from_top_layer(p) = qflx_solidevap_orig
+ end if
+
! conserve total energy flux
eflx_sh_grnd(p) = eflx_sh_grnd(p) +(evaporation_demand -evaporation_limit)*htvp(c)
endif

(Note that the diffs in
ERP_P36x2_Ld30.f45_f45_mg37.I2000Clm51FatesSpCruRsGs.cheyenne_intel.clm-FatesColdSatPhen
were ambiguous as to whether they were roundoff-level due to the
single-precision output in that test; I reran with double precision
for the baseline and the branch and was able to verify that the
diffs were only double-precision roundoff-level.)

Other details
-------------

Pull Requests that document the changes (include PR ids):
https://github.com/ESCOMP/CTSM/pull/1987

===============================================================
===============================================================
Tag name: ctsm5.1.dev121
Originator(s): glemieux (Gregory Lemieux,LBL/NGEET,510-486-5049)
Date: Wed Apr 5 13:34:09 MDT 2023
Expand Down
1 change: 1 addition & 0 deletions doc/ChangeSum
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
Tag Who Date Summary
============================================================================================================================
ctsm5.1.dev122 sacks 04/23/2023 Rework handling of evaporation constraint in SoilFluxes
ctsm5.1.dev121 glemieux 04/05/2023 Changes soil moisture initialization logic for FATES
ctsm5.1.dev120 sacks 03/25/2023 Update externals and minor fixes
ctsm5.1.dev119 slevis 03/16/2023 Allow gross unrepresented land use transition (PR #309)
Expand Down
159 changes: 86 additions & 73 deletions src/biogeophys/SoilFluxesMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,74 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
endif
end do

! Partition evaporation into liquid and solid
do fp = 1, num_nolakep
p = filter_nolakep(fp)
c = patch%column(p)
l = patch%landunit(p)
j = col%snl(c)+1

qflx_liqevap_from_top_layer(p) = 0._r8
qflx_solidevap_from_top_layer(p) = 0._r8
qflx_soliddew_to_top_layer(p) = 0._r8
qflx_liqdew_to_top_layer(p) = 0._r8

! Partition the evaporation from snow/soil surface into liquid evaporation,
! solid evaporation (sublimation), liquid dew, or solid dew. Note that the variables
! affected here are all related to the snow subgrid patch only because of the use of qflx_ev_snow.
! In the situations where there are snow layers or there is snow without an explicit snow layer,
! the partitioned variables will represent the components of snow evaporation
! (qflx_ev_snow = qflx_liqevap_from_top_layer + qflx_solidevap_from_top_layer
! - qflx_liqdew_to_top_layer - qflx_soliddew_to_top_layer).
! In the case of no snow, qflx_ev_snow has already been set equal to qflx_ev_soil (the evaporation
! from the subgrid soil patch) and the partitioned variables will then represent evaporation from the
! subgrid soil patch.
! In the case of urban columns (and lake columns - see LakeHydrologyMod), there are no subgrid
! patches and qflx_evap_soi is used. qflx_evap_soi = qflx_liqevap_from_top_layer + qflx_solidevap_from_top_layer
! - qflx_liqdew_to_top_layer - qflx_soliddew_to_top_layer.
if (.not. lun%urbpoi(l)) then
if (qflx_ev_snow(p) >= 0._r8) then
! for evaporation partitioning between liquid evap and ice sublimation,
! use the ratio of liquid to (liquid+ice) in the top layer to determine split
if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
qflx_liqevap_from_top_layer(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/ &
(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
else
qflx_liqevap_from_top_layer(p) = 0._r8
end if
qflx_solidevap_from_top_layer(p) = qflx_ev_snow(p) - qflx_liqevap_from_top_layer(p)
else
if (t_grnd(c) < tfrz) then
qflx_soliddew_to_top_layer(p) = abs(qflx_ev_snow(p))
else
qflx_liqdew_to_top_layer(p) = abs(qflx_ev_snow(p))
end if
end if

else ! Urban columns

if (qflx_evap_soi(p) >= 0._r8) then
! for evaporation partitioning between liquid evap and ice sublimation,
! use the ratio of liquid to (liquid+ice) in the top layer to determine split
if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
qflx_liqevap_from_top_layer(p) = max(qflx_evap_soi(p)*(h2osoi_liq(c,j)/ &
(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
else
qflx_liqevap_from_top_layer(p) = 0._r8
end if
qflx_solidevap_from_top_layer(p) = qflx_evap_soi(p) - qflx_liqevap_from_top_layer(p)
else
if (t_grnd(c) < tfrz) then
qflx_soliddew_to_top_layer(p) = abs(qflx_evap_soi(p))
else
qflx_liqdew_to_top_layer(p) = abs(qflx_evap_soi(p))
end if
end if

end if

end do

! Constrain evaporation from snow to be <= available moisture
do fp = 1,num_nolakep
p = filter_nolakep(fp)
Expand All @@ -221,6 +289,8 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
evaporation_demand = qflx_ev_snow(p)
qflx_ev_snow(p) = evaporation_limit
qflx_evap_soi(p) = qflx_evap_soi(p) - frac_sno_eff(c)*(evaporation_demand - evaporation_limit)
qflx_liqevap_from_top_layer(p) = max(h2osoi_liq(c,j)/(frac_sno_eff(c)*dtime), 0._r8)
qflx_solidevap_from_top_layer(p) = max(h2osoi_ice(c,j)/(frac_sno_eff(c)*dtime), 0._r8)
! conserve total energy flux
eflx_sh_grnd(p) = eflx_sh_grnd(p) + frac_sno_eff(c)*(evaporation_demand - evaporation_limit)*htvp(c)
endif
Expand All @@ -235,11 +305,27 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
evaporation_demand = qflx_evap_soi(p)
qflx_evap_soi(p) = evaporation_limit
qflx_ev_snow(p) = qflx_evap_soi(p)
qflx_liqevap_from_top_layer(p) = max(h2osoi_liq(c,j)/dtime, 0._r8)
qflx_solidevap_from_top_layer(p) = max(h2osoi_ice(c,j)/dtime, 0._r8)
! conserve total energy flux
eflx_sh_grnd(p) = eflx_sh_grnd(p) +(evaporation_demand -evaporation_limit)*htvp(c)
endif
endif

! limit only solid evaporation (sublimation) from top soil layer
! (liquid evaporation from soil should not be limited)
if (j==1 .and. frac_h2osfc(c) < 1._r8) then
evaporation_limit = h2osoi_ice(c,j)/(dtime*(1._r8 - frac_h2osfc(c)))
if (qflx_solidevap_from_top_layer(p) >= evaporation_limit) then
evaporation_demand = qflx_solidevap_from_top_layer(p)
qflx_solidevap_from_top_layer(p) &
= evaporation_limit
qflx_liqevap_from_top_layer(p) &
= qflx_liqevap_from_top_layer(p) &
+ (evaporation_demand - evaporation_limit)
endif
endif

enddo

call t_stopf('bgp2_loop_1')
Expand Down Expand Up @@ -299,79 +385,6 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
eflx_sh_tot_u(p)= eflx_sh_tot(p)
end if

qflx_liqevap_from_top_layer(p) = 0._r8
qflx_solidevap_from_top_layer(p) = 0._r8
qflx_soliddew_to_top_layer(p) = 0._r8
qflx_liqdew_to_top_layer(p) = 0._r8

! Partition the evaporation from snow/soil surface into liquid evaporation,
! solid evaporation (sublimation), liquid dew, or solid dew. Note that the variables
! affected here are all related to the snow subgrid patch only because of the use of qflx_ev_snow.
! In the situations where there are snow layers or there is snow without an explicit snow layer,
! the partitioned variables will represent the components of snow evaporation
! (qflx_ev_snow = qflx_liqevap_from_top_layer + qflx_solidevap_from_top_layer
! - qflx_liqdew_to_top_layer - qflx_soliddew_to_top_layer).
! In the case of no snow, qflx_ev_snow has already been set equal to qflx_ev_soil (the evaporation
! from the subgrid soil patch) and the partitioned variables will then represent evaporation from the
! subgrid soil patch.
! In the case of urban columns (and lake columns - see LakeHydrologyMod), there are no subgrid
! patches and qflx_evap_soi is used. qflx_evap_soi = qflx_liqevap_from_top_layer + qflx_solidevap_from_top_layer
! - qflx_liqdew_to_top_layer - qflx_soliddew_to_top_layer.
if (.not. lun%urbpoi(l)) then
if (qflx_ev_snow(p) >= 0._r8) then
! for evaporation partitioning between liquid evap and ice sublimation,
! use the ratio of liquid to (liquid+ice) in the top layer to determine split
if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
qflx_liqevap_from_top_layer(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/ &
(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
else
qflx_liqevap_from_top_layer(p) = 0._r8
end if
qflx_solidevap_from_top_layer(p) = qflx_ev_snow(p) - qflx_liqevap_from_top_layer(p)
else
if (t_grnd(c) < tfrz) then
qflx_soliddew_to_top_layer(p) = abs(qflx_ev_snow(p))
else
qflx_liqdew_to_top_layer(p) = abs(qflx_ev_snow(p))
end if
end if

else ! Urban columns

if (qflx_evap_soi(p) >= 0._r8) then
! for evaporation partitioning between liquid evap and ice sublimation,
! use the ratio of liquid to (liquid+ice) in the top layer to determine split
if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
qflx_liqevap_from_top_layer(p) = max(qflx_evap_soi(p)*(h2osoi_liq(c,j)/ &
(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
else
qflx_liqevap_from_top_layer(p) = 0._r8
end if
qflx_solidevap_from_top_layer(p) = qflx_evap_soi(p) - qflx_liqevap_from_top_layer(p)
else
if (t_grnd(c) < tfrz) then
qflx_soliddew_to_top_layer(p) = abs(qflx_evap_soi(p))
else
qflx_liqdew_to_top_layer(p) = abs(qflx_evap_soi(p))
end if
end if

end if

! limit only solid evaporation (sublimation) from top soil layer
! (liquid evaporation from soil should not be limited)
if (j==1 .and. frac_h2osfc(c) < 1._r8) then
evaporation_limit = h2osoi_ice(c,j)/(dtime*(1._r8 - frac_h2osfc(c)))
if (qflx_solidevap_from_top_layer(p) >= evaporation_limit) then
evaporation_demand = qflx_solidevap_from_top_layer(p)
qflx_solidevap_from_top_layer(p) &
= evaporation_limit
qflx_liqevap_from_top_layer(p) &
= qflx_liqevap_from_top_layer(p) &
+ (evaporation_demand - evaporation_limit)
endif
endif

! Variables needed by history tape

qflx_evap_can(p) = qflx_evap_veg(p) - qflx_tran_veg(p)
Expand Down

0 comments on commit 798c29c

Please sign in to comment.