Merge pull request #13 from NCAR/truncate_h2osfc

Truncate small h2osfc values to zero

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.

ChangeLog_branch

@@ -1,5 +1,124 @@
 ===============================================================
-Tag name: ctsm_n08_clm4_5_16_r249
+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

cime_config/testdefs/testlist_clm.xml

@@ -485,8 +485,8 @@
         <machine compiler="intel" testtype="aux_clm" testmods="clm/default">cheyenne</machine>
         <machine compiler="intel" testtype="aux_clm" testmods="clm/default">yellowstone</machine>
       </test>
-      <test name="SMS_Lm25">
-        <machine compiler="gnu" testtype="aux_clm" testmods="clm/cropMonthOutput">yellowstone</machine>
+      <test name="SMS_Lm13">
+        <machine compiler="gnu" testtype="aux_clm" testmods="clm/cropMonthOutput" comment="include a relatively long crop test at relatively high resolution">yellowstone</machine>
       </test>
     </grid>
     <grid name="f19_g17_gl4">

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
 

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

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)

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

src/utils/test/CMakeLists.txt

@@ -1,2 +1,3 @@
 add_subdirectory(clm_time_manager_test)
-add_subdirectory(annual_flux_dribbler_test)
+add_subdirectory(annual_flux_dribbler_test)
+add_subdirectory(numerics_test)

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)