diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90 index 98a327a8a4..4886d5e23e 100644 --- a/src/biogeophys/SoilHydrologyMod.F90 +++ b/src/biogeophys/SoilHydrologyMod.F90 @@ -17,6 +17,7 @@ module SoilHydrologyMod use column_varcon , only : icol_road_imperv use landunit_varcon , only : istsoil, istcrop use clm_time_manager , only : get_step_size + use NumericsMod , only : truncate_small_values use EnergyFluxType , only : energyflux_type use InfiltrationExcessRunoffMod, only : infiltration_excess_runoff_type use SoilHydrologyType , only : soilhydrology_type @@ -214,12 +215,12 @@ subroutine SetQflxInputs(bounds, num_hydrologyc, filter_hydrologyc, & qflx_top_soil => waterflux_inst%qflx_top_soil_col , & ! Output: [real(r8) (:)] net water input into soil from top (mm/s) qflx_in_soil => waterflux_inst%qflx_in_soil_col , & ! Output: [real(r8) (:)] surface input to soil (mm/s) qflx_top_soil_to_h2osfc => waterflux_inst%qflx_top_soil_to_h2osfc_col , & ! Output: [real(r8) (:)] portion of qflx_top_soil going to h2osfc, minus evaporation (mm/s) - qflx_rain_plus_snomelt => waterflux_inst%qflx_rain_plus_snomelt_col , & ! Input: [real(r8) (:)] rain plus snow melt falling on the soil (mm/s) + qflx_rain_plus_snomelt => waterflux_inst%qflx_rain_plus_snomelt_col , & ! Input: [real(r8) (:)] rain plus snow melt falling on the soil (mm/s) qflx_snow_h2osfc => waterflux_inst%qflx_snow_h2osfc_col , & ! Input: [real(r8) (:)] snow falling on surface water (mm/s) qflx_floodc => waterflux_inst%qflx_floodc_col , & ! Input: [real(r8) (:)] column flux of flood water from RTM - qflx_ev_soil => waterflux_inst%qflx_ev_soil_col , & ! Input: [real(r8) (:) ] evaporation flux from soil (W/m**2) [+ to atm] - qflx_evap_grnd => waterflux_inst%qflx_evap_grnd_col , & ! Input: [real(r8) (:) ] ground surface evaporation rate (mm H2O/s) [+] - qflx_ev_h2osfc => waterflux_inst%qflx_ev_h2osfc_col , & ! Input: [real(r8) (:) ] evaporation flux from h2osfc (W/m**2) [+ to atm] + qflx_ev_soil => waterflux_inst%qflx_ev_soil_col , & ! Input: [real(r8) (:)] evaporation flux from soil (W/m**2) [+ to atm] + qflx_evap_grnd => waterflux_inst%qflx_evap_grnd_col , & ! Input: [real(r8) (:)] ground surface evaporation rate (mm H2O/s) [+] + qflx_ev_h2osfc => waterflux_inst%qflx_ev_h2osfc_col , & ! Input: [real(r8) (:)] evaporation flux from h2osfc (W/m**2) [+ to atm] qflx_sat_excess_surf => saturated_excess_runoff_inst%qflx_sat_excess_surf_col, & ! Input: [real(r8) (:) ] surface runoff due to saturated surface (mm H2O /s) @@ -345,7 +346,6 @@ subroutine UpdateH2osfc(bounds, num_hydrologyc, filter_hydrologyc, & integer :: c,l,fc ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: h2osfc_partial(bounds%begc:bounds%endc) ! partially-updated h2osfc - logical :: truncate_h2osfc_to_zero(bounds%begc:bounds%endc) !----------------------------------------------------------------------- associate( & @@ -383,36 +383,32 @@ subroutine UpdateH2osfc(bounds, num_hydrologyc, filter_hydrologyc, & h2osfc_partial(c) = h2osfc(c) + (qflx_in_h2osfc(c) - qflx_h2osfc_surf(c)) * dtime end do + call truncate_small_values(num_f = num_hydrologyc, filter_f = filter_hydrologyc, & + lb = bounds%begc, ub = bounds%endc, & + data_baseline = h2osfc(bounds%begc:bounds%endc), & + data = h2osfc_partial(bounds%begc:bounds%endc)) + call QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & h2osfcflag = h2osfcflag, & h2osfc = h2osfc_partial(bounds%begc:bounds%endc), & frac_h2osfc = frac_h2osfc(bounds%begc:bounds%endc), & qinmax = qinmax(bounds%begc:bounds%endc), & - qflx_h2osfc_drain = qflx_h2osfc_drain(bounds%begc:bounds%endc), & - truncate_h2osfc_to_zero = truncate_h2osfc_to_zero(bounds%begc:bounds%endc)) + qflx_h2osfc_drain = qflx_h2osfc_drain(bounds%begc:bounds%endc)) + ! Update h2osfc based on fluxes do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) - - ! The parenthesization of this expression is just needed to maintain bfb answers - ! in the major refactor. Note that the first parenthesized expression is - ! h2osfc_partial(c), but I'm writing it out explicitly to facilitate a possible - ! future removal of h2osfc_partial. - h2osfc(c) = (h2osfc(c) + (qflx_in_h2osfc(c) - qflx_h2osfc_surf(c)) * dtime) & - - qflx_h2osfc_drain(c) * dtime - - ! TODO(wjs, 2017-08-22) This is here to maintain bit-for-bit answers with the old - ! code. If we're okay changing answers, we could remove it. Or maybe we want to - ! put a more general truncation in here, like: - ! if (abs(h2osfc(c)) < 1.e-14_r8 * abs(h2osfc_orig)) h2osfc(c) = 0._r8 - ! But I'm not sure what the general philosophy is regarding whether and when we - ! want to do truncations like this. - if (truncate_h2osfc_to_zero(c)) then - h2osfc(c) = 0._r8 - end if - + h2osfc(c) = h2osfc_partial(c) - qflx_h2osfc_drain(c) * dtime end do + ! Due to rounding errors, fluxes that should have brought h2osfc to exactly 0 may + ! have instead left it slightly less than or slightly greater than 0. Correct for + ! that here. + call truncate_small_values(num_f = num_hydrologyc, filter_f = filter_hydrologyc, & + lb = bounds%begc, ub = bounds%endc, & + data_baseline = h2osfc_partial(bounds%begc:bounds%endc), & + data = h2osfc(bounds%begc:bounds%endc)) + end associate end subroutine UpdateH2osfc @@ -524,16 +520,13 @@ end subroutine QflxH2osfcSurf !----------------------------------------------------------------------- subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & h2osfcflag, h2osfc, frac_h2osfc, qinmax, & - qflx_h2osfc_drain, truncate_h2osfc_to_zero) + qflx_h2osfc_drain) ! ! !DESCRIPTION: ! Compute qflx_h2osfc_drain ! ! Note that, if h2osfc is negative, then qflx_h2osfc_drain will be negative - acting - ! to exactly restore h2osfc to 0. In this case, truncate_h2osfc_to_zero will be set - ! to true; in all other cases, truncate_h2osfc_to_zero will be set to false. This - ! particular behavior of truncate_h2osfc_to_zero is just like this to maintain - ! bit-for-bit answers with the old code (see also the TODO note in UpdateH2osfc). + ! to exactly restore h2osfc to 0. ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -544,7 +537,6 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & real(r8) , intent(in) :: frac_h2osfc( bounds%begc: ) ! fraction of ground covered by surface water (0 to 1) real(r8) , intent(in) :: qinmax( bounds%begc: ) ! maximum infiltration rate (mm H2O /s) real(r8) , intent(inout) :: qflx_h2osfc_drain( bounds%begc: ) ! bottom drainage from h2osfc (mm H2O /s) - logical , intent(inout) :: truncate_h2osfc_to_zero( bounds%begc: ) ! whether h2osfc should be truncated to 0 to correct for roundoff errors, in order to maintain bit-for-bit the same answers as the old code ! ! !LOCAL VARIABLES: integer :: fc, c @@ -557,7 +549,6 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & SHR_ASSERT_ALL((ubound(frac_h2osfc) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) SHR_ASSERT_ALL((ubound(qinmax) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) SHR_ASSERT_ALL((ubound(qflx_h2osfc_drain) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) - SHR_ASSERT_ALL((ubound(truncate_h2osfc_to_zero) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) dtime = get_step_size() @@ -566,14 +557,12 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & if (h2osfc(c) < 0.0) then qflx_h2osfc_drain(c) = h2osfc(c)/dtime - truncate_h2osfc_to_zero(c) = .true. else qflx_h2osfc_drain(c)=min(frac_h2osfc(c)*qinmax(c),h2osfc(c)/dtime) if(h2osfcflag==0) then ! ensure no h2osfc qflx_h2osfc_drain(c)= max(0._r8,h2osfc(c)/dtime) end if - truncate_h2osfc_to_zero(c) = .false. end if end do diff --git a/src/unit_test_shr/unittestFilterBuilderMod.F90 b/src/unit_test_shr/unittestFilterBuilderMod.F90 index 3ed61eb9a2..6999b63d72 100644 --- a/src/unit_test_shr/unittestFilterBuilderMod.F90 +++ b/src/unit_test_shr/unittestFilterBuilderMod.F90 @@ -55,7 +55,7 @@ subroutine filter_from_range(start, end, numf, filter) ! !LOCAL VARIABLES: integer :: i - character(len=*), parameter :: subname = 'filter_from_list' + character(len=*), parameter :: subname = 'filter_from_range' !----------------------------------------------------------------------- numf = end - start + 1 diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index 31a298e28e..4c4c498734 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -17,6 +17,7 @@ list(APPEND clm_sources clm_nlUtilsMod.F90 clm_time_manager.F90 fileutils.F90 + NumericsMod.F90 ) sourcelist_to_parent(clm_sources) diff --git a/src/utils/NumericsMod.F90 b/src/utils/NumericsMod.F90 new file mode 100644 index 0000000000..8d75753bb8 --- /dev/null +++ b/src/utils/NumericsMod.F90 @@ -0,0 +1,89 @@ +module NumericsMod + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Utility routines for assisting with model numerics + ! + ! !USES: +#include "shr_assert.h" + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_log_mod , only : errMsg => shr_log_errMsg + + implicit none + save + private + + ! !PUBLIC MEMBER FUNCTIONS: + public :: truncate_small_values ! Truncate relatively small values to 0 + + ! !PUBLIC MEMBER DATA: + + ! Relative differences below rel_epsilon are considered to be zero. + ! + ! Note that double precision machine epsilon is approximately 1e-16, so this value of + ! 1e-13 allows for 3 orders of magnitude of "slop". + ! + ! Examples of how to use this: + ! + ! (1) Rather than checking + ! if (x == y) + ! instead check + ! if (abs(x - y) < rel_epsilon * x) + ! or + ! if (abs(x - y) < rel_epsilon * y) + ! + ! (2) After a state update, you can truncate the state to 0 based on this condition: + ! if (abs(some_state) < rel_epsilon * abs(some_state_orig)) then + ! some_state = 0._r8 + ! end if + ! where some_state_orig is the value of the state variable before the update + real(r8), public, parameter :: rel_epsilon = 1.e-13_r8 ! Relative differences below this are considered to be zero + + ! !PRIVATE MEMBER DATA: + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +contains + + !----------------------------------------------------------------------- + subroutine truncate_small_values(num_f, filter_f, lb, ub, data_baseline, data) + ! + ! !DESCRIPTION: + ! Truncate relatively small values to 0, within the given filter. + ! + ! "Relatively small" is determined by comparison with some "baseline" version of the + ! data. + ! + ! For example, this can be used after doing a state update. In this case, + ! data_baseline should hold the values before the state update, and data should hold + ! the values after the state update. + ! + ! !ARGUMENTS: + integer , intent(in) :: num_f ! number of points in filter_f + integer , intent(in) :: filter_f(:) ! filter of points in data + integer , intent(in) :: lb ! lower bound of data + integer , intent(in) :: ub ! upper bound of data + real(r8) , intent(in) :: data_baseline(lb:) ! baseline version of data, used to define "relatively close to 0" + real(r8) , intent(inout) :: data(lb:) ! data to operate on + ! + ! !LOCAL VARIABLES: + integer :: fn ! index into filter + integer :: n ! index into data + + character(len=*), parameter :: subname = 'truncate_small_values' + !----------------------------------------------------------------------- + + SHR_ASSERT_ALL((ubound(data_baseline) == (/ub/)), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL((ubound(data) == (/ub/)), errMsg(sourcefile, __LINE__)) + + do fn = 1, num_f + n = filter_f(fn) + if (abs(data(n)) < rel_epsilon * abs(data_baseline(n))) then + data(n) = 0._r8 + end if + end do + + end subroutine truncate_small_values + +end module NumericsMod diff --git a/src/utils/test/CMakeLists.txt b/src/utils/test/CMakeLists.txt index 64e7e89171..2b39a131f7 100644 --- a/src/utils/test/CMakeLists.txt +++ b/src/utils/test/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(clm_time_manager_test) -add_subdirectory(annual_flux_dribbler_test) \ No newline at end of file +add_subdirectory(annual_flux_dribbler_test) +add_subdirectory(numerics_test) diff --git a/src/utils/test/numerics_test/CMakeLists.txt b/src/utils/test/numerics_test/CMakeLists.txt new file mode 100644 index 0000000000..83bf0cc1d0 --- /dev/null +++ b/src/utils/test/numerics_test/CMakeLists.txt @@ -0,0 +1,10 @@ +set (pfunit_sources + test_truncate_small_values.pf) + +set (extra_sources + ) + +create_pFUnit_test(numerics test_numerics_exe + "${pfunit_sources}" "${extra_sources}") + +target_link_libraries(test_numerics_exe clm csm_share esmf_wrf_timemgr) diff --git a/src/utils/test/numerics_test/test_truncate_small_values.pf b/src/utils/test/numerics_test/test_truncate_small_values.pf new file mode 100644 index 0000000000..a989f21c3d --- /dev/null +++ b/src/utils/test/numerics_test/test_truncate_small_values.pf @@ -0,0 +1,120 @@ +module test_truncate_small_values + + ! Tests of NumericsMod: truncate_small_values + + use pfunit_mod + use NumericsMod + use shr_kind_mod , only : r8 => shr_kind_r8 + use unittestSimpleSubgridSetupsMod + use unittestSubgridMod + use unittestFilterBuilderMod, only : filter_from_range + + implicit none + + @TestCase + type, extends(TestCase) :: TestTSV + contains + procedure :: setUp + procedure :: tearDown + end type TestTSV + + real(r8), parameter :: tol = 1.e-13_r8 + +contains + + subroutine setUp(this) + class(TestTSV), intent(inout) :: this + end subroutine setUp + + subroutine tearDown(this) + class(TestTSV), intent(inout) :: this + + call unittest_subgrid_teardown() + end subroutine tearDown + + @Test + subroutine truncates_correct_points(this) + class(TestTSV), intent(inout) :: this + real(r8) :: data_baseline(3) + real(r8) :: data(3) + real(r8) :: data_saved(3) + integer :: num_f + integer, allocatable :: filter_f(:) + + call setup_n_veg_patches(pwtcol = [0.1_r8, 0.8_r8, 0.1_r8], pft_types = [1, 2, 3]) + call filter_from_range(bounds%begp, bounds%endp, num_f, filter_f) + + ! point 2 should be truncated, others should not be truncated + data_baseline = [1._r8, 1._r8, 1._r8] + data = [0.5_r8, 1.e-16_r8, -1._r8] + data_saved = data + + call truncate_small_values( & + num_f = num_f, & + filter_f = filter_f, & + lb = bounds%begp, & + ub = bounds%endp, & + data_baseline = data_baseline, & + data = data) + + @assertEqual(data_saved(1), data(1)) + @assertEqual(data_saved(3), data(3)) + @assertEqual(0._r8, data(2)) + + end subroutine truncates_correct_points + + @Test + subroutine truncates_large_magnitude(this) + ! Make sure we're just relying on relative rather than absolute magnitudes by + ! confirming that it can truncate a value with large magnitude. + class(TestTSV), intent(inout) :: this + real(r8) :: data_baseline(1) + real(r8) :: data(1) + integer :: num_f + integer, allocatable :: filter_f(:) + + call setup_single_veg_patch(pft_type = 1) + call filter_from_range(bounds%begp, bounds%endp, num_f, filter_f) + + data_baseline = [1.e30_r8] + data = [1.e10_r8] + + call truncate_small_values( & + num_f = num_f, & + filter_f = filter_f, & + lb = bounds%begp, & + ub = bounds%endp, & + data_baseline = data_baseline, & + data = data) + + @assertEqual(0._r8, data(1)) + end subroutine truncates_large_magnitude + + @Test + subroutine does_not_truncate_small_magnitude(this) + ! Make sure we're just relying on relative rather than absolute magnitudes by + ! confirming that it does not truncate a value with small magnitude. + class(TestTSV), intent(inout) :: this + real(r8) :: data_baseline(1) + real(r8) :: data(1) + integer :: num_f + integer, allocatable :: filter_f(:) + + call setup_single_veg_patch(pft_type = 1) + call filter_from_range(bounds%begp, bounds%endp, num_f, filter_f) + + data_baseline = [1.e-30_r8] + data = [1.e-31_r8] + + call truncate_small_values( & + num_f = num_f, & + filter_f = filter_f, & + lb = bounds%begp, & + ub = bounds%endp, & + data_baseline = data_baseline, & + data = data) + + @assertEqual(1.e-31_r8, data(1)) + end subroutine does_not_truncate_small_magnitude + +end module test_truncate_small_values