From c8b80d84d732444a710c5f23fd1d331ce95dcba0 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Sat, 8 Apr 2023 09:36:56 -0600 Subject: [PATCH 1/3] Rework handling of evaporation constraint in SoilFluxes 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 the deleted lines in this commit have just been moved as-is to earlier in this subroutine. --- src/biogeophys/SoilFluxesMod.F90 | 159 +++++++++++++++++-------------- 1 file changed, 86 insertions(+), 73 deletions(-) diff --git a/src/biogeophys/SoilFluxesMod.F90 b/src/biogeophys/SoilFluxesMod.F90 index 5f4030c6e1..48820df41b 100644 --- a/src/biogeophys/SoilFluxesMod.F90 +++ b/src/biogeophys/SoilFluxesMod.F90 @@ -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) @@ -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) = h2osoi_liq(c,j)/(frac_sno_eff(c)*dtime) + qflx_solidevap_from_top_layer(p) = h2osoi_ice(c,j)/(frac_sno_eff(c)*dtime) ! conserve total energy flux eflx_sh_grnd(p) = eflx_sh_grnd(p) + frac_sno_eff(c)*(evaporation_demand - evaporation_limit)*htvp(c) endif @@ -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) = h2osoi_liq(c,j)/dtime + qflx_solidevap_from_top_layer(p) = h2osoi_ice(c,j)/dtime ! 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') @@ -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) From 1c099f3a752a2d4a07ea0c98ccae76a46ad96ccb Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Mon, 17 Apr 2023 15:14:33 -0600 Subject: [PATCH 2/3] Don't allow qflx_liqevap and qflx_solidevap fluxes to be negative These could sometimes be roundoff-level negative due to h2osoi_liq starting out roundoff-level negative. --- src/biogeophys/SoilFluxesMod.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/biogeophys/SoilFluxesMod.F90 b/src/biogeophys/SoilFluxesMod.F90 index 48820df41b..c316d30fe3 100644 --- a/src/biogeophys/SoilFluxesMod.F90 +++ b/src/biogeophys/SoilFluxesMod.F90 @@ -289,8 +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) = h2osoi_liq(c,j)/(frac_sno_eff(c)*dtime) - qflx_solidevap_from_top_layer(p) = h2osoi_ice(c,j)/(frac_sno_eff(c)*dtime) + 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 @@ -305,8 +305,8 @@ 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) = h2osoi_liq(c,j)/dtime - qflx_solidevap_from_top_layer(p) = h2osoi_ice(c,j)/dtime + 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 From 4e58ce37dd8cce4a0fe10ed96594c606635682be Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Sun, 23 Apr 2023 19:48:52 -0600 Subject: [PATCH 3/3] Update ChangeLog --- doc/ChangeLog | 166 ++++++++++++++++++++++++++++++++++++++++++++++++++ doc/ChangeSum | 1 + 2 files changed, 167 insertions(+) diff --git a/doc/ChangeLog b/doc/ChangeLog index edb6223e2e..f25b413a38 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -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 diff --git a/doc/ChangeSum b/doc/ChangeSum index 24a0b303c6..faadd8cee6 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -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)