diff --git a/ChangeLog_branch b/ChangeLog_branch index 8f65aed1f1..ace21cfdb9 100644 --- a/ChangeLog_branch +++ b/ChangeLog_branch @@ -1,4 +1,216 @@ =============================================================== +Tag name: ctsm_n10_clm4_5_16_r249 +Originator(s): sacks +Date: Sep 14, 2017 +One-line Summary: Truncate small h2osfc values to zero + +Purpose of changes +------------------ + +This replaces the earlier code that only did this truncation in one +particular circumstance. + +This seems mainly important to truncate small negative numbers to 0, but +it also seems like a good idea to truncate small positive numbers that +should have been 0. + +With some intermediate commits on this branch, I did some careful tests +to ensure that these diffs introduce no more than roundoff-level +changes. + + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- + +CLM testing: + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + yellowstone - not run + + unit-tests (components/clm/src): + + yellowstone - pass + + tools-tests (components/clm/test/tools): + + yellowstone - not run + + PTCLM testing (components/clm/tools/shared/PTCLM/test): + + yellowstone - not run + + regular tests (aux_clm): + + yellowstone_intel - ok + yellowstone_pgi - ok + yellowstone_gnu - ok + cheyenne_intel - ok + hobart_nag - ok + hobart_pgi - ok + hobart_intel - ok + + ok means tests pass, change answers as expected + +CLM tag used for the baseline comparisons: ctsm_n09_clm4_5_16_r249 + + +Answer changes +-------------- + +Changes answers relative to baseline: YES + + If a tag changes answers relative to baseline comparison the + following should be filled in (otherwise remove this section): + + Summarize any changes to answers, i.e., + - what code configurations: all clm45/clm50 + - what platforms/compilers: all + - nature of change (roundoff; larger than roundoff/same climate; new climate): + roundoff + + If bitwise differences were observed, how did you show they were no worse + than roundoff? + + I introduced some temporary code to ensure that, where there are + differences between the old and new code, these differences are + only roundoff-level (e.g., a roundoff-level value being rounded to + 0 now where it wasn't before). + + If this tag changes climate describe the run(s) done to evaluate the new + climate (put details of the simulations in the experiment database) + - casename: N/A + + URL for LMWG diagnostics output used to validate new climate: N/A + + +Detailed list of changes +------------------------ + +List any svn externals directories updated (cime, rtm, mosart, cism, etc.): none + +List all files eliminated: none + +List all files added and what they do: + +========= New utility for truncating small values to 0, and associated + unit tests +A components/clm/src/utils/NumericsMod.F90 +A components/clm/src/utils/test/numerics_test/CMakeLists.txt +A components/clm/src/utils/test/numerics_test/test_truncate_small_values.pf + +List all existing files that have been modified, and describe the changes: + +========= Main changes +M components/clm/src/biogeophys/SoilHydrologyMod.F90 + +========= New unit tests +M components/clm/src/utils/CMakeLists.txt +M components/clm/src/utils/test/CMakeLists.txt + +========= Unrelated fix to subroutine name +M components/clm/src/unit_test_shr/unittestFilterBuilderMod.F90 + +========= Shorten a test that was taking > 2 hours +M components/clm/cime_config/testdefs/testlist_clm.xml + +=============================================================== +=============================================================== +Tag name: ctsm_n09_clm4_5_16_r249 +Originator(s): sacks +Date: Sep 7, 2017 +One-line Summary: Use full qflx_surf in BGC code + +Purpose of changes +------------------ + +We had been using a separate variable to avoid changing answers. But +Dave Lawrence gave his approval to change this. + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- + +CLM testing: + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + yellowstone - not run + + unit-tests (components/clm/src): + + yellowstone - not run + + tools-tests (components/clm/test/tools): + + yellowstone - not run + + PTCLM testing (components/clm/tools/shared/PTCLM/test): + + yellowstone - not run + + regular tests (aux_clm): + + yellowstone_intel - ok + yellowstone_pgi - ok + yellowstone_gnu - ok + cheyenne_intel - ok + hobart_nag - ok + hobart_pgi - ok + hobart_intel - ok + + ok means tests pass, answers change as expected + +CLM tag used for the baseline comparisons: ctsm_n08_clm4_5_16_r249 + + +Answer changes +-------------- + +Changes answers relative to baseline: YES + + If a tag changes answers relative to baseline comparison the + following should be filled in (otherwise remove this section): + + Summarize any changes to answers, i.e., + - what code configurations: All Bgc configurations + - what platforms/compilers: all + - nature of change (roundoff; larger than roundoff/same climate; new climate): + not investigated carefully; larger than roundoff, but expected to + be same climate + + If bitwise differences were observed, how did you show they were no worse + than roundoff? N/A + + If this tag changes climate describe the run(s) done to evaluate the new + climate (put details of the simulations in the experiment database) + - casename: N/A + + URL for LMWG diagnostics output used to validate new climate: N/A + + +Detailed list of changes +------------------------ + +List any svn externals directories updated (cime, rtm, mosart, cism, etc.): none + +List all files eliminated: none + +List all files added and what they do: none + +List all existing files that have been modified, and describe the changes: + +========= Main changes, as described above +M components/clm/src/biogeochem/ch4Mod.F90 +M components/clm/src/biogeophys/SoilHydrologyMod.F90 +M components/clm/src/biogeophys/WaterfluxType.F90 +M components/clm/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 + +=============================================================== +=============================================================== Tag name: ctsm_n08_clm4_5_16_r249 Originator(s): sacks Date: Aug 24, 2017 diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 53e942b2d0..a6f2480a8f 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -485,8 +485,8 @@ cheyenne yellowstone - - yellowstone + + yellowstone diff --git a/src/biogeochem/ch4Mod.F90 b/src/biogeochem/ch4Mod.F90 index 787aebbc24..575269deb9 100644 --- a/src/biogeochem/ch4Mod.F90 +++ b/src/biogeochem/ch4Mod.F90 @@ -1725,7 +1725,7 @@ subroutine ch4 (bounds, num_soilc, filter_soilc, num_lakec, filter_lakec, & frac_h2osfc => waterstate_inst%frac_h2osfc_col , & ! Input: [real(r8) (:) ] fraction of ground covered by surface water (0 to 1) snow_depth => waterstate_inst%snow_depth_col , & ! Input: [real(r8) (:) ] snow height (m) tws => waterstate_inst%tws_grc , & ! Input: [real(r8) (:) ] total water storage (kg m-2) - qflx_surf => waterflux_inst%qflx_surf_for_bgc_col , & ! Input: [real(r8) (:) ] surface runoff for input into BGC code (mm H2O /s) + qflx_surf => waterflux_inst%qflx_surf_col , & ! Input: [real(r8) (:) ] total surface runoff (mm H2O /s) conc_o2_sat => ch4_inst%conc_o2_sat_col , & ! Input: [real(r8) (:,:) ] O2 conc in each soil layer (mol/m3) (nlevsoi) zwt0 => ch4_inst%zwt0_col , & ! Input: [real(r8) (:) ] decay factor for finundated (m) diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90 index 19b6bfdebe..30cc8a60a7 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 @@ -231,12 +232,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) @@ -362,7 +363,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( & @@ -400,36 +400,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 @@ -569,16 +565,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 @@ -589,7 +582,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: type(drainPond_type) :: drainPond_inst @@ -608,7 +600,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() @@ -617,7 +608,6 @@ 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 ! define maximum drainage @@ -654,7 +644,6 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & ! 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 @@ -710,7 +699,6 @@ subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, & snl => col%snl , & ! Input: [integer (:) ] minus number of snow layers qflx_surf => waterflux_inst%qflx_surf_col , & ! Output: [real(r8) (:) ] total surface runoff (mm H2O /s) - qflx_surf_for_bgc => waterflux_inst%qflx_surf_for_bgc_col , & ! Input: [real(r8) (:) ] surface runoff for input into BGC code (mm H2O /s) qflx_infl_excess_surf => waterflux_inst%qflx_infl_excess_surf_col, & ! Input: [real(r8) (:) ] surface runoff due to infiltration excess (mm H2O /s) qflx_h2osfc_surf => waterflux_inst%qflx_h2osfc_surf_col, & ! Input: [real(r8) (:) ] surface water runoff (mm H2O /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) @@ -735,16 +723,6 @@ subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, & ! Depending on whether h2osfcflag is 0 or 1, one of qflx_infl_excess or ! qflx_h2osfc_surf will always be 0. But it's safe to just add them both. qflx_surf(c) = qflx_sat_excess_surf(c) + qflx_infl_excess_surf(c) + qflx_h2osfc_surf(c) - - ! TODO(wjs, 2017-07-11) I'm distinguishing between qflx_surf and qflx_surf_for_bgc - ! simply to maintain answers as they were before. But I have a feeling that the - ! BGC code should really be using the total qflx_surf in its calculations. Once - ! Dave Lawrence or someone else signs off on this change, we should change the BGC - ! code to use qflx_surf and remove this qflx_surf_for_bgc variable. - ! Alternatively, if we deem the current implementation correct, we should - ! consider renaming this something better than qflx_surf_for_bgc, or simply - ! making the BGC code depend on qflx_sat_excess_surf, if that's what's intended. - qflx_surf_for_bgc(c) = qflx_sat_excess_surf(c) + qflx_infl_excess_surf(c) end do ! ------------------------------------------------------------------------ diff --git a/src/biogeophys/WaterfluxType.F90 b/src/biogeophys/WaterfluxType.F90 index ff55279596..8d6e76eaa1 100644 --- a/src/biogeophys/WaterfluxType.F90 +++ b/src/biogeophys/WaterfluxType.F90 @@ -76,7 +76,6 @@ module WaterfluxType real(r8), pointer :: qflx_infl_excess_surf_col(:) ! col surface runoff due to infiltration excess (mm H2O /s) real(r8), pointer :: qflx_h2osfc_surf_col (:) ! col surface water runoff (mm H2O /s) real(r8), pointer :: qflx_surf_col (:) ! col total surface runoff (mm H2O /s) - real(r8), pointer :: qflx_surf_for_bgc_col (:) ! col total surface runoff for input into BGC code (mm H2O /s) real(r8), pointer :: qflx_drain_col (:) ! col sub-surface runoff (mm H2O /s) real(r8), pointer :: qflx_rain_plus_snomelt_col(:) ! col rain plus snow melt falling on the soil (mm/s) real(r8), pointer :: qflx_top_soil_col (:) ! col net water input into soil from top (mm/s) @@ -223,7 +222,6 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_rootsoi_col (begc:endc,1:nlevsoi)) ; this%qflx_rootsoi_col (:,:) = nan allocate(this%qflx_infl_col (begc:endc)) ; this%qflx_infl_col (:) = nan allocate(this%qflx_surf_col (begc:endc)) ; this%qflx_surf_col (:) = nan - allocate(this%qflx_surf_for_bgc_col (begc:endc)) ; this%qflx_surf_for_bgc_col (:) = nan allocate(this%qflx_drain_col (begc:endc)) ; this%qflx_drain_col (:) = nan allocate(this%qflx_rain_plus_snomelt_col(begc:endc)) ; this%qflx_rain_plus_snomelt_col(:) = nan allocate(this%qflx_top_soil_col (begc:endc)) ; this%qflx_top_soil_col (:) = nan diff --git a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 index 2cc77aef17..72db6363bb 100644 --- a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 @@ -109,7 +109,7 @@ subroutine SoilBiogeochemNLeaching(bounds, num_soilc, filter_soilc, & h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) (new) (-nlevsno+1:nlevgrnd) qflx_drain => waterflux_inst%qflx_drain_col , & ! Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s) - qflx_surf => waterflux_inst%qflx_surf_for_bgc_col , & ! Input: [real(r8) (:) ] surface runoff for input into BGC code (mm H2O /s) + qflx_surf => waterflux_inst%qflx_surf_col , & ! Input: [real(r8) (:) ] total surface runoff (mm H2O /s) sminn_vr => soilbiogeochem_nitrogenstate_inst%sminn_vr_col , & ! Input: [real(r8) (:,:) ] (gN/m3) soil mineral N smin_no3_vr => soilbiogeochem_nitrogenstate_inst%smin_no3_vr_col , & ! Input: [real(r8) (:,:) ] 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