Skip to content

Commit

Permalink
Merge pull request ESCOMP#13 from NCAR/truncate_h2osfc
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
billsacks committed Dec 28, 2017
2 parents dea2996 + 7372be2 commit dd8f549
Show file tree
Hide file tree
Showing 7 changed files with 246 additions and 36 deletions.
57 changes: 23 additions & 34 deletions src/biogeophys/SoilHydrologyMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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( &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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()

Expand All @@ -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

Expand Down
2 changes: 1 addition & 1 deletion src/unit_test_shr/unittestFilterBuilderMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/utils/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ list(APPEND clm_sources
clm_nlUtilsMod.F90
clm_time_manager.F90
fileutils.F90
NumericsMod.F90
)

sourcelist_to_parent(clm_sources)
89 changes: 89 additions & 0 deletions src/utils/NumericsMod.F90
Original file line number Diff line number Diff line change
@@ -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
3 changes: 2 additions & 1 deletion src/utils/test/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
10 changes: 10 additions & 0 deletions src/utils/test/numerics_test/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
120 changes: 120 additions & 0 deletions src/utils/test/numerics_test/test_truncate_small_values.pf
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit dd8f549

Please sign in to comment.