diff --git a/autotest/TestMathUtil.f90 b/autotest/TestMathUtil.f90
index a1fc90a1233..b28e9050d8e 100644
--- a/autotest/TestMathUtil.f90
+++ b/autotest/TestMathUtil.f90
@@ -1,8 +1,9 @@
module TestMathUtil
use KindModule, only: I4B, DP
+ use ConstantsModule, only: DNODATA, DZERO
use testdrive, only: check, error_type, new_unittest, test_failed, &
to_string, unittest_type
- use MathUtilModule, only: mod_offset
+ use MathUtilModule, only: is_close, mod_offset
implicit none
private
public :: collect_mathutil
@@ -12,6 +13,9 @@ module TestMathUtil
subroutine collect_mathutil(testsuite)
type(unittest_type), allocatable, intent(out) :: testsuite(:)
testsuite = [ &
+ new_unittest("is_close_symmetric", test_is_close_symmetric), &
+ new_unittest("is_close_symmetric_near_0", &
+ test_is_close_symmetric_near_0), &
new_unittest("mod_offset", &
test_mod_offset) &
]
@@ -41,4 +45,123 @@ subroutine test_mod_offset(error)
call check(error, mod_offset(2.0_DP, 3.0_DP, 3.0_DP) == 5.0_DP)
end subroutine test_mod_offset
+ subroutine test_is_close_symmetric(error)
+ type(error_type), allocatable, intent(out) :: error
+ real(DP) :: a, b, rtol
+
+ ! exact match
+ a = 1.0_DP
+ b = 1.0_DP
+ call check(error, is_close(a, b), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", eps=default")
+ if (allocated(error)) return
+
+ ! mismatch with default rtol
+ b = 1.0001_DP
+ call check(error, (.not. (is_close(a, b))), &
+ "exp ne: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", eps=default")
+ if (allocated(error)) return
+
+ ! inexact match with large rtol
+ rtol = 1d-2
+ call check(error, is_close(a, b, rtol=rtol), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", rtol="//to_string(rtol))
+ if (allocated(error)) return
+
+ ! mismatch when we reduce rtol
+ rtol = 0.5d-5
+ call check(error, (.not. is_close(a, b, rtol=rtol)), &
+ "exp ne: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", rtol="//to_string(rtol))
+ if (allocated(error)) return
+
+ ! +/-0
+ call check(error, is_close(0.0_DP, -0.0_DP), &
+ "exp ne: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", eps=default")
+
+ ! DNODATA
+ call check(error, (.not. is_close(0.0_DP, DNODATA)), &
+ "exp ne: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", eps=default")
+ call check(error, is_close(DNODATA, DNODATA), &
+ "exp ne: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", eps=default")
+ call check(error, (.not. is_close(DNODATA, DNODATA / 10)), &
+ "exp ne: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", eps=default")
+ call check(error, (.not. is_close(DNODATA, DNODATA * 10)), &
+ "exp ne: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", eps=default")
+
+ end subroutine test_is_close_symmetric
+
+ subroutine test_is_close_symmetric_near_0(error)
+ type(error_type), allocatable, intent(out) :: error
+ real(DP) :: a, b, rtol, atol
+
+ a = 0.0_DP
+ b = 0.0_DP
+ call check(error, is_close(a, b), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", rtol=default")
+ if (allocated(error)) return
+
+ a = DZERO
+ b = DZERO
+ call check(error, is_close(a, b), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", rtol=default")
+ if (allocated(error)) return
+
+ b = 1d-4
+ call check(error, (.not. is_close(a, b)), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", rtol=default")
+ if (allocated(error)) return
+
+ rtol = 0.999_DP
+ call check(error, &
+ ! expect failure, see above
+ (.not. is_close(a, b, rtol=rtol)), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", rtol="//to_string(rtol))
+ if (allocated(error)) return
+
+ ! absolute comparison is appropriate when a and/or b are near or equal to 0
+ b = 1d-4
+ atol = 1d-3
+ call check(error, is_close(a, b, atol=atol), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", atol="//to_string(atol))
+ if (allocated(error)) return
+
+ ! make sure the absolute tolerance is applied
+ b = 1d-4
+ atol = 1d-5
+ call check(error, (.not. is_close(a, b, atol=atol)), &
+ "exp eq: a="//to_string(a)// &
+ ", b="//to_string(b)// &
+ ", atol="//to_string(atol))
+ if (allocated(error)) return
+
+ end subroutine test_is_close_symmetric_near_0
+
end module TestMathUtil
diff --git a/make/makefile b/make/makefile
index cd2e35d8612..e0c00d9cbf2 100644
--- a/make/makefile
+++ b/make/makefile
@@ -85,7 +85,7 @@ $(OBJDIR)/compilerversion.o \
$(OBJDIR)/version.o \
$(OBJDIR)/Sim.o \
$(OBJDIR)/OpenSpec.o \
-$(OBJDIR)/genericutils.o \
+$(OBJDIR)/MathUtil.o \
$(OBJDIR)/InputOutput.o \
$(OBJDIR)/TableTerm.o \
$(OBJDIR)/Table.o \
@@ -324,7 +324,6 @@ $(OBJDIR)/BaseGeometry.o \
$(OBJDIR)/mf6.o \
$(OBJDIR)/StringList.o \
$(OBJDIR)/MemorySetHandler.o \
-$(OBJDIR)/MathUtil.o \
$(OBJDIR)/ilut.o \
$(OBJDIR)/sparsekit.o \
$(OBJDIR)/rcm.o \
diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj
index 6c78c5313af..7d547f60d79 100644
--- a/msvs/mf6core.vfproj
+++ b/msvs/mf6core.vfproj
@@ -384,7 +384,6 @@
-
diff --git a/msvs/mf6lib.vfproj b/msvs/mf6lib.vfproj
index 7186195f1e3..dee2f209199 100644
--- a/msvs/mf6lib.vfproj
+++ b/msvs/mf6lib.vfproj
@@ -143,7 +143,6 @@
-
diff --git a/pymake/makefile b/pymake/makefile
index 484e7798cf0..8e7065a0d6d 100644
--- a/pymake/makefile
+++ b/pymake/makefile
@@ -114,7 +114,6 @@ $(OBJDIR)/GwfNpfOptions.o \
$(OBJDIR)/GwfBuyInputData.o \
$(OBJDIR)/ims8misc.o \
$(OBJDIR)/LinearSolverBase.o \
-$(OBJDIR)/genericutils.o \
$(OBJDIR)/ArrayHandlers.o \
$(OBJDIR)/IndexMap.o \
$(OBJDIR)/version.o \
diff --git a/src/Model/Connection/GridSorting.f90 b/src/Model/Connection/GridSorting.f90
index 3360b62809d..5d15fa2382c 100644
--- a/src/Model/Connection/GridSorting.f90
+++ b/src/Model/Connection/GridSorting.f90
@@ -2,7 +2,7 @@ module GridSorting
use KindModule, only: I4B, DP, LGP
use ConstantsModule, only: DHALF
use CellWithNbrsModule, only: GlobalCellType
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use BaseDisModule, only: dis_transform_xy
implicit none
private
@@ -75,11 +75,11 @@ function lessThan(n, m) result(isLess)
dis_bot_m(gcm%index))
! compare
- if (.not. is_same(zn, zm, 10 * epsilon(zn))) then
+ if (.not. is_close(zn, zm, 10 * epsilon(zn))) then
isLess = zn > zm
- else if (.not. is_same(yn, ym, 10 * epsilon(yn))) then
+ else if (.not. is_close(yn, ym, 10 * epsilon(yn))) then
isLess = yn > ym
- else if (.not. is_same(xn, xm, 10 * epsilon(xn))) then
+ else if (.not. is_close(xn, xm, 10 * epsilon(xn))) then
isLess = xn < xm
else
isLess = .false.
diff --git a/src/Model/GroundWaterFlow/gwf3csub8.f90 b/src/Model/GroundWaterFlow/gwf3csub8.f90
index 73029f582db..db87b8ff411 100644
--- a/src/Model/GroundWaterFlow/gwf3csub8.f90
+++ b/src/Model/GroundWaterFlow/gwf3csub8.f90
@@ -18,7 +18,7 @@ module GwfCsubModule
TABLEFT, TABCENTER, TABRIGHT, &
TABSTRING, TABUCSTRING, TABINTEGER, TABREAL
use MemoryHelperModule, only: create_mem_path
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use MessageModule, only: write_message
use SmoothingModule, only: sQuadraticSaturation, &
sQuadraticSaturationDerivative, &
@@ -2875,7 +2875,7 @@ subroutine csub_fc(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
rhs(node) = rhs(node) + rhsterm
!
! -- calculate interbed water compressibility terms
- if (this%brg /= DZERO .and. idelay == 0) then
+ if (.not. is_close(this%brg, DZERO) .and. idelay == 0) then
call this%csub_nodelay_wcomp_fc(ib, node, tled, area, &
hnew(node), hold(node), &
hcof, rhsterm)
diff --git a/src/Model/GroundWaterFlow/gwf3lak8.f90 b/src/Model/GroundWaterFlow/gwf3lak8.f90
index 52ffc9e2916..3834f640263 100644
--- a/src/Model/GroundWaterFlow/gwf3lak8.f90
+++ b/src/Model/GroundWaterFlow/gwf3lak8.f90
@@ -28,7 +28,7 @@ module LakModule
use BaseDisModule, only: DisBaseType
use SimModule, only: count_errors, store_error, store_error_unit, &
deprecation_warning
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use BlockParserModule, only: BlockParserType
use BaseDisModule, only: DisBaseType
use SimVariablesModule, only: errmsg, warnmsg
@@ -833,7 +833,7 @@ subroutine lak_read_lake_connections(this)
warnmsg, this%parser%GetUnit())
case default
read (keyword, *) rval
- if (is_same(rval, DNODATA)) then
+ if (is_close(rval, DNODATA)) then
is_lake_bed = .FALSE.
else
is_lake_bed = .TRUE.
@@ -1832,7 +1832,7 @@ subroutine lak_read_initial_attr(this)
end if
length = this%connlength(j)
end if
- if (is_same(this%bedleak(j), DNODATA)) then
+ if (is_close(this%bedleak(j), DNODATA)) then
clb(j) = DNODATA
else if (this%bedleak(j) > DZERO) then
clb(j) = DONE / this%bedleak(j)
@@ -1844,7 +1844,7 @@ subroutine lak_read_initial_attr(this)
else
caq(j) = DZERO
end if
- if (is_same(this%bedleak(j), DNODATA)) then
+ if (is_close(this%bedleak(j), DNODATA)) then
this%satcond(j) = area / caq(j)
else if (clb(j) * caq(j) > DZERO) then
this%satcond(j) = area / (clb(j) + caq(j))
@@ -1879,7 +1879,7 @@ subroutine lak_read_initial_attr(this)
nn = this%cellid(j)
area = this%warea(j)
c1 = DZERO
- if (is_same(clb(j), DNODATA)) then
+ if (is_close(clb(j), DNODATA)) then
cbedleak = ' NONE '
cbedcond = ' NONE '
else if (clb(j) > DZERO) then
@@ -2604,7 +2604,7 @@ subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
call this%lak_calculate_sarea(ilak, stage, sa)
ev = sa * this%evaporation(ilak)
if (ev > avail) then
- if (is_same(avail, DPREC)) then
+ if (is_close(avail, DPREC)) then
ev = DZERO
else
ev = -avail
diff --git a/src/Solution/LinearMethods/ims8base.f90 b/src/Solution/LinearMethods/ims8base.f90
index 72bc98eb32b..4c82d9f87af 100644
--- a/src/Solution/LinearMethods/ims8base.f90
+++ b/src/Solution/LinearMethods/ims8base.f90
@@ -9,7 +9,7 @@ MODULE IMSLinearBaseModule
use KindModule, only: DP, I4B
use ConstantsModule, only: LINELENGTH, IZERO, &
DZERO, DPREC, DEM6, DEM3, DHALF, DONE
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use BlockParserModule, only: BlockParserType
use IMSReorderingModule, only: ims_odrv
use ConvergenceSummaryModule
@@ -211,7 +211,7 @@ SUBROUTINE ims_base_cg(ICNVG, ITMAX, INNERIT, &
IF (ICNVG .NE. 0) EXIT INNER
!
! -- CHECK THAT CURRENT AND PREVIOUS rho ARE DIFFERENT
- lsame = is_same(rho, rho0)
+ lsame = is_close(rho, rho0)
IF (lsame) THEN
EXIT INNER
END IF
@@ -513,15 +513,15 @@ SUBROUTINE ims_base_bcgs(ICNVG, ITMAX, INNERIT, &
!
! -- CHECK THAT CURRENT AND PREVIOUS rho, alpha, AND omega ARE
! DIFFERENT
- lsame = is_same(rho, rho0)
+ lsame = is_close(rho, rho0)
IF (lsame) THEN
EXIT INNER
END IF
- lsame = is_same(alpha, alpha0)
+ lsame = is_close(alpha, alpha0)
IF (lsame) THEN
EXIT INNER
END IF
- lsame = is_same(omega, omega0)
+ lsame = is_close(omega, omega0)
IF (lsame) THEN
EXIT INNER
END IF
diff --git a/src/Solution/NumericalSolution.f90 b/src/Solution/NumericalSolution.f90
index 456b58d54bd..41bc43e576b 100644
--- a/src/Solution/NumericalSolution.f90
+++ b/src/Solution/NumericalSolution.f90
@@ -13,8 +13,8 @@ module NumericalSolutionModule
LENMEMPATH
use MemoryHelperModule, only: create_mem_path
use TableModule, only: TableType, table_cr
- use GenericUtilitiesModule, only: is_same
Use MessageModule, only: write_message
+ use MathUtilModule, only: is_close
use VersionModule, only: IDEVELOPMODE
use BaseModelModule, only: BaseModelType
use BaseExchangeModule, only: BaseExchangeType
@@ -2476,7 +2476,7 @@ subroutine sln_ls(this, kiter, kstp, kper, in_iter, iptc, ptcf)
end if
end if
else
- lsame = is_same(l2norm, this%l2norm0)
+ lsame = is_close(l2norm, this%l2norm0)
if (lsame) then
iptc = 0
end if
diff --git a/src/Utilities/InputOutput.f90 b/src/Utilities/InputOutput.f90
index 7253b996f3e..913ae044d5b 100644
--- a/src/Utilities/InputOutput.f90
+++ b/src/Utilities/InputOutput.f90
@@ -11,7 +11,6 @@ module InputOutputModule
TABLEFT, TABCENTER, TABRIGHT, &
TABSTRING, TABUCSTRING, TABINTEGER, TABREAL, &
DZERO
- use GenericUtilitiesModule, only: is_same
use MessageModule, only: write_message
private
public :: GetUnit, &
diff --git a/src/Utilities/MathUtil.f90 b/src/Utilities/MathUtil.f90
index a33d3332d39..ad6af25ddd0 100644
--- a/src/Utilities/MathUtil.f90
+++ b/src/Utilities/MathUtil.f90
@@ -7,7 +7,7 @@ module MathUtilModule
implicit none
private
- public :: mod_offset
+ public :: mod_offset, is_close
interface mod_offset
module procedure :: mod_offset_int, mod_offset_dbl
@@ -15,6 +15,67 @@ module MathUtilModule
contains
+ !> @brief Check if a real value is approximately equal to another.
+ !!
+ !! By default the determination is symmetric in a and b, as in
+ !! Python's math.isclose, with relative difference scaled by a
+ !! factor of the larger absolute value of a and b. The formula
+ !! is: abs(a - b) <= max(rtol * max(abs(a), abs(b)), atol).
+ !!
+ !! If symmetric is set to false the test is asymmetric in a and
+ !! b, with b taken to be the reference value, and the alternate
+ !! formula (abs(a - b) <= (atol + rtol * abs(b))) is used. This
+ !! is the approach taken by numpy.allclose.
+ !!
+ !! Defaults for rtol and atol are DSAME and DZERO, respectively.
+ !! If a or b are near 0 (especially if either is 0), an absolute
+ !! tolerance suitable for the particular case should be provided.
+ !! For a justification of a zero absolute tolerance default see:
+ !! https://peps.python.org/pep-0485/#absolute-tolerance-default
+ !<
+ pure logical function is_close(a, b, rtol, atol, symmetric)
+ ! dummy
+ real(DP), intent(in) :: a !< first real
+ real(DP), intent(in) :: b !< second real (reference value if asymmetric)
+ real(DP), intent(in), optional :: rtol !< relative tolerance (default=DSAME)
+ real(DP), intent(in), optional :: atol !< absolute tolerance (default=DZERO)
+ logical(LGP), intent(in), optional :: symmetric !< toggle (a)symmetric comparison
+ ! local
+ real(DP) :: lrtol, latol
+ logical(LGP) :: lsymmetric
+
+ ! check for exact equality
+ if (a == b) then
+ is_close = .true.
+ return
+ end if
+
+ ! process optional arguments
+ if (.not. present(rtol)) then
+ lrtol = DSAME
+ else
+ lrtol = rtol
+ end if
+ if (.not. present(atol)) then
+ latol = DZERO
+ else
+ latol = atol
+ end if
+ if (.not. present(symmetric)) then
+ lsymmetric = .true.
+ else
+ lsymmetric = symmetric
+ end if
+
+ if (lsymmetric) then
+ ! "weak" symmetric test, https://peps.python.org/pep-0485/#which-symmetric-test
+ is_close = abs(a - b) <= max(lrtol * max(abs(a), abs(b)), latol)
+ else
+ ! asymmetric, https://numpy.org/doc/stable/reference/generated/numpy.isclose.html
+ is_close = (abs(a - b) <= (latol + lrtol * abs(b)))
+ end if
+ end function is_close
+
!> @brief Modulo with offset for integer values.
pure function mod_offset_int(a, n, d) result(mo)
! -- dummy
diff --git a/src/Utilities/TimeSeries/TimeArraySeries.f90 b/src/Utilities/TimeSeries/TimeArraySeries.f90
index 24e1125c31c..7c8c4607bbc 100644
--- a/src/Utilities/TimeSeries/TimeArraySeries.f90
+++ b/src/Utilities/TimeSeries/TimeArraySeries.f90
@@ -4,7 +4,7 @@ module TimeArraySeriesModule
use BlockParserModule, only: BlockParserType
use ConstantsModule, only: LINELENGTH, UNDEFINED, STEPWISE, LINEAR, &
LENTIMESERIESNAME, LENMODELNAME, DZERO, DONE
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use InputOutputModule, only: GetUnit, openfile
use KindModule, only: DP, I4B
use ListModule, only: ListType, ListNodeType
@@ -474,7 +474,7 @@ subroutine get_values_at_time(this, nvals, values, time)
ierr = 1
end if
else
- if (is_same(taEarlier%taTime, time)) then
+ if (is_close(taEarlier%taTime, time)) then
values = taEarlier%taArray
else
! -- Only earlier time is available, and it is not time of interest;
@@ -488,7 +488,7 @@ subroutine get_values_at_time(this, nvals, values, time)
end if
else
if (associated(taLater)) then
- if (is_same(taLater%taTime, time)) then
+ if (is_close(taLater%taTime, time)) then
values = taLater%taArray
else
! -- only later time is available, and it is not time of interest
@@ -717,7 +717,7 @@ subroutine get_latest_preceding_node(this, time, tslNode)
if (associated(currNode%nextNode)) then
obj => currNode%nextNode%GetItem()
ta => CastAsTimeArrayType(obj)
- if (ta%taTime < time .or. is_same(ta%taTime, time)) then
+ if (ta%taTime < time .or. is_close(ta%taTime, time)) then
currNode => currNode%nextNode
else
exit
diff --git a/src/Utilities/TimeSeries/TimeSeries.f90 b/src/Utilities/TimeSeries/TimeSeries.f90
index 50904b0b9cd..b4a5d9d48f2 100644
--- a/src/Utilities/TimeSeries/TimeSeries.f90
+++ b/src/Utilities/TimeSeries/TimeSeries.f90
@@ -5,7 +5,7 @@ module TimeSeriesModule
use ConstantsModule, only: LINELENGTH, UNDEFINED, STEPWISE, LINEAR, &
LINEAREND, LENTIMESERIESNAME, LENHUGELINE, &
DZERO, DONE, DNODATA
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use InputOutputModule, only: GetUnit, openfile, ParseLine, upcase
use ListModule, only: ListType, ListNodeType
use SimVariablesModule, only: errmsg
@@ -319,7 +319,7 @@ subroutine get_surrounding_records(this, time, tsrecEarlier, tsrecLater)
if (associated(currNode%nextNode)) then
obj => currNode%nextNode%GetItem()
tsr => CastAsTimeSeriesRecordType(obj)
- if (tsr%tsrTime < time .and. .not. is_same(tsr%tsrTime, time)) then
+ if (tsr%tsrTime < time .and. .not. is_close(tsr%tsrTime, time)) then
currNode => currNode%nextNode
else
exit
@@ -356,7 +356,7 @@ subroutine get_surrounding_records(this, time, tsrecEarlier, tsrecLater)
obj => tsNode1%GetItem()
tsrec1 => CastAsTimeSeriesRecordType(obj)
time1 = tsrec1%tsrTime
- do while (time1 < time .and. .not. is_same(time1, time))
+ do while (time1 < time .and. .not. is_close(time1, time))
if (associated(tsNode1%nextNode)) then
tsNode1 => tsNode1%nextNode
obj => tsNode1%GetItem()
@@ -373,8 +373,8 @@ subroutine get_surrounding_records(this, time, tsrecEarlier, tsrecLater)
!
end if
!
- if (time0 < time .or. is_same(time0, time)) tsrecEarlier => tsrec0
- if (time1 > time .or. is_same(time1, time)) tsrecLater => tsrec1
+ if (time0 < time .or. is_close(time0, time)) tsrecEarlier => tsrec0
+ if (time1 > time .or. is_close(time1, time)) tsrecLater => tsrec1
!
! -- Return
return
@@ -418,7 +418,7 @@ subroutine get_surrounding_nodes(this, time, nodeEarlier, nodeLater)
if (associated(currNode%nextNode)) then
obj => currNode%nextNode%GetItem()
tsr => CastAsTimeSeriesRecordType(obj)
- if (tsr%tsrTime < time .and. .not. is_same(tsr%tsrTime, time)) then
+ if (tsr%tsrTime < time .and. .not. is_close(tsr%tsrTime, time)) then
currNode => currNode%nextNode
else
exit
@@ -454,7 +454,7 @@ subroutine get_surrounding_nodes(this, time, nodeEarlier, nodeLater)
obj => tsNode1%GetItem()
tsrec1 => CastAsTimeSeriesRecordType(obj)
time1 = tsrec1%tsrTime
- do while (time1 < time .and. .not. is_same(time1, time))
+ do while (time1 < time .and. .not. is_close(time1, time))
if (associated(tsNode1%nextNode)) then
tsNode1 => tsNode1%nextNode
obj => tsNode1%GetItem()
@@ -467,11 +467,11 @@ subroutine get_surrounding_nodes(this, time, nodeEarlier, nodeLater)
!
end if
!
- if (time0 < time .or. is_same(time0, time)) then
+ if (time0 < time .or. is_close(time0, time)) then
tsrecEarlier => tsrec0
nodeEarlier => tsNode0
end if
- if (time1 > time .or. is_same(time1, time)) then
+ if (time1 > time .or. is_close(time1, time)) then
tsrecLater => tsrec1
nodeLater => tsNode1
end if
@@ -552,7 +552,7 @@ function get_value_at_time(this, time, extendToEndOfSimulation)
ierr = 1
end if
else
- if (extendToEndOfSimulation .or. is_same(tsrEarlier%tsrTime, time)) then
+ if (extendToEndOfSimulation .or. is_close(tsrEarlier%tsrTime, time)) then
get_value_at_time = tsrEarlier%tsrValue
else
! -- Only earlier time is available, and it is not time of interest;
@@ -566,7 +566,7 @@ function get_value_at_time(this, time, extendToEndOfSimulation)
end if
else
if (associated(tsrLater)) then
- if (is_same(tsrLater%tsrTime, time)) then
+ if (is_close(tsrLater%tsrTime, time)) then
get_value_at_time = tsrLater%tsrValue
else
! -- only later time is available, and it is not time of interest
@@ -624,7 +624,7 @@ function get_integrated_value(this, time0, time1, extendToEndOfSimulation)
currObj => currNode%GetItem()
currRecord => CastAsTimeSeriesRecordType(currObj)
currTime = currRecord%tsrTime
- if (is_same(currTime, time1)) then
+ if (is_close(currTime, time1)) then
! Current node time = time1 so should be ldone
ldone = .true.
elseif (currTime < time1) then
@@ -657,12 +657,12 @@ function get_integrated_value(this, time0, time1, extendToEndOfSimulation)
if (lprocess) then
! -- determine lower and upper limits of time span of interest
! within current interval
- if (currTime > time0 .or. is_same(currTime, time0)) then
+ if (currTime > time0 .or. is_close(currTime, time0)) then
t0 = currTime
else
t0 = time0
end if
- if (nextTime < time1 .or. is_same(nextTime, time1)) then
+ if (nextTime < time1 .or. is_close(nextTime, time1)) then
t1 = nextTime
else
t1 = time1
@@ -697,7 +697,7 @@ function get_integrated_value(this, time0, time1, extendToEndOfSimulation)
! -- Are we done yet?
if (t1 > time1) then
ldone = .true.
- elseif (is_same(t1, time1)) then
+ elseif (is_close(t1, time1)) then
ldone = .true.
else
! -- We are not done yet
@@ -796,7 +796,7 @@ subroutine get_latest_preceding_node(this, time, tslNode)
if (associated(currNode%nextNode)) then
obj => currNode%nextNode%GetItem()
tsr => CastAsTimeSeriesRecordType(obj)
- if (tsr%tsrTime < time .or. is_same(tsr%tsrTime, time)) then
+ if (tsr%tsrTime < time .or. is_close(tsr%tsrTime, time)) then
currNode => currNode%nextNode
else
exit
@@ -829,7 +829,7 @@ subroutine get_latest_preceding_node(this, time, tslNode)
end do
end if
!
- if (time0 < time .or. is_same(time0, time)) tslNode => tsNode0
+ if (time0 < time .or. is_close(time0, time)) tslNode => tsNode0
!
! -- Return
return
@@ -946,7 +946,7 @@ function GetTimeSeriesRecord(this, time, epsi) result(res)
do
tsr => this%GetNextTimeSeriesRecord()
if (associated(tsr)) then
- if (is_same(tsr%tsrTime, time)) then
+ if (is_close(tsr%tsrTime, time)) then
res => tsr
exit
end if
diff --git a/src/Utilities/genericutils.f90 b/src/Utilities/genericutils.f90
deleted file mode 100644
index e6921da9dd6..00000000000
--- a/src/Utilities/genericutils.f90
+++ /dev/null
@@ -1,66 +0,0 @@
-!> @brief This module contains generic utilties
-!!
-!! This module contains generic utilities that have
-!! limited dependencies.
-!!
-!<
-module GenericUtilitiesModule
- use KindModule, only: DP, I4B, LGP
- use ConstantsModule, only: MAXCHARLEN, LENHUGELINE, &
- DZERO, DPREC, DSAME, &
- LINELENGTH, LENHUGELINE, VSUMMARY
- use SimVariablesModule, only: istdout, isim_level
- !
- implicit none
-
- private
-
- public :: is_same
-
-contains
-
- !> @brief Function to determine if two reals are the same
- !!
- !! Function to evaluate if the difference between a and b are less than eps
- !! (i.e. a and b are the same).
- !!
- !<
- function is_same(a, b, eps) result(lvalue)
- ! -- return variable
- logical(LGP) :: lvalue !< boolean indicating if a and b are the same
- ! -- dummy variables
- real(DP), intent(in) :: a !< first number to evaluate
- real(DP), intent(in) :: b !< second number to evaluate
- real(DP), intent(in), optional :: eps !< optional maximum difference between a abd b (default=DSAME)
- ! -- local variables
- real(DP) :: epsloc
- real(DP) :: denom
- real(DP) :: rdiff
- !
- ! -- evaluate optioanl arguments
- if (present(eps)) then
- epsloc = eps
- else
- epsloc = DSAME
- end if
- lvalue = .FALSE.
- if (a == b) then
- lvalue = .TRUE.
- else
- if (abs(b) > abs(a)) then
- denom = b
- else
- denom = a
- if (abs(denom) == DZERO) then
- denom = DPREC
- end if
- end if
- rdiff = abs((a - b) / denom)
- if (rdiff <= epsloc) then
- lvalue = .TRUE.
- end if
- end if
-
- end function is_same
-
-end module GenericUtilitiesModule
diff --git a/src/meson.build b/src/meson.build
index 800accc6553..e77d9891412 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -227,7 +227,6 @@ modflow_sources = files(
'Utilities' / 'defmacro.F90',
'Utilities' / 'DevFeature.f90',
'Utilities' / 'ErrorUtil.f90',
- 'Utilities' / 'genericutils.f90',
'Utilities' / 'GeomUtil.f90',
'Utilities' / 'HashTable.f90',
'Utilities' / 'HeadFileReader.f90',
diff --git a/utils/mf5to6/make/makefile b/utils/mf5to6/make/makefile
index 00f87156629..fd923f00ba0 100644
--- a/utils/mf5to6/make/makefile
+++ b/utils/mf5to6/make/makefile
@@ -37,14 +37,14 @@ $(OBJDIR)/compilerversion.o \
$(OBJDIR)/version.o \
$(OBJDIR)/OpenSpec.o \
$(OBJDIR)/GlobalVariables.o \
+$(OBJDIR)/ErrorUtil.o \
$(OBJDIR)/SimPHMF.o \
-$(OBJDIR)/genericutils.o \
+$(OBJDIR)/MathUtil.o \
$(OBJDIR)/InputOutput.o \
$(OBJDIR)/TableTerm.o \
$(OBJDIR)/Table.o \
$(OBJDIR)/MemoryHelper.o \
$(OBJDIR)/CharString.o \
-$(OBJDIR)/ErrorUtil.o \
$(OBJDIR)/Memory.o \
$(OBJDIR)/List.o \
$(OBJDIR)/MemoryList.o \
diff --git a/utils/mf5to6/msvs/mf5to6.vfproj b/utils/mf5to6/msvs/mf5to6.vfproj
index f206d7d30aa..10ff7233751 100644
--- a/utils/mf5to6/msvs/mf5to6.vfproj
+++ b/utils/mf5to6/msvs/mf5to6.vfproj
@@ -103,7 +103,6 @@
-
diff --git a/utils/mf5to6/pymake/extrafiles.txt b/utils/mf5to6/pymake/extrafiles.txt
index 6b5a9b976c4..f8a3888f945 100644
--- a/utils/mf5to6/pymake/extrafiles.txt
+++ b/utils/mf5to6/pymake/extrafiles.txt
@@ -9,7 +9,6 @@
../../../src/Utilities/Constants.f90
../../../src/Utilities/SimVariables.f90
../../../src/Utilities/compilerversion.F90
-../../../src/Utilities/genericutils.f90
../../../src/Utilities/Message.f90
../../../src/Utilities/ErrorUtil.f90
../../../src/Utilities/GeomUtil.f90
@@ -17,6 +16,7 @@
../../../src/Utilities/kind.f90
../../../src/Utilities/List.f90
../../../src/Utilities/LongLineReader.f90
+../../../src/Utilities/MathUtil.f90
../../../src/Utilities/OpenSpec.f90
../../../src/Utilities/defmacro.F90
../../../src/Utilities/version.f90
diff --git a/utils/mf5to6/src/MultiLayerObsModule.f90 b/utils/mf5to6/src/MultiLayerObsModule.f90
index f8c77064f46..e2c87430a8e 100644
--- a/utils/mf5to6/src/MultiLayerObsModule.f90
+++ b/utils/mf5to6/src/MultiLayerObsModule.f90
@@ -1,8 +1,8 @@
module MultiLayerObs
-
- use ConstantsModule, only: DONE, MAXCHARLEN
- use ConstantsPHMFModule, only: LENOBSNAMENEW
- use GenericUtilitiesModule, only: is_same
+
+ use ConstantsModule, only: DONE, MAXCHARLEN
+ use ConstantsPHMFModule, only: LENOBSNAMENEW
+ use MathUtilModule, only: is_close
use ListModule, only: ListType
use SimModule, only: store_error, ustop
@@ -26,7 +26,7 @@ module MultiLayerObs
procedure, public :: CheckWeightSum
end type
- contains
+contains
! Non-type-bound procedures
@@ -43,7 +43,7 @@ function CastAsLayerObsType(obj) result(res)
end select
return
end function CastAsLayerObsType
-
+
subroutine ConstructLayerObs(newLayerObs, layobname, layer, weight)
! dummy
type(LayerObsType), pointer :: newLayerObs
@@ -51,7 +51,7 @@ subroutine ConstructLayerObs(newLayerObs, layobname, layer, weight)
integer, intent(in) :: layer
double precision, intent(in) :: weight
!
- allocate(newLayerObs)
+ allocate (newLayerObs)
newLayerObs%lobsname = layobname
newLayerObs%layer = layer
newLayerObs%weight = weight
@@ -71,7 +71,7 @@ subroutine AddLayerObsToList(list, layerobs)
!
return
end subroutine AddLayerObsToList
-
+
function GetLayerObsFromList(list, indx) result(res)
! dummy
type(ListType), intent(inout) :: list
@@ -88,18 +88,18 @@ function GetLayerObsFromList(list, indx) result(res)
!
return
end function GetLayerObsFromList
-
+
subroutine ConstructMLObs(newMLObs, obsname)
! dummy
type(MLObsType), pointer, intent(inout) :: newMLObs
character(len=LENOBSNAMENEW), intent(in) :: obsname
!
- allocate(newMLObs)
+ allocate (newMLObs)
newMLObs%mlobsname = obsname
!
return
end subroutine ConstructMLObs
-
+
subroutine AddMLObsToList(list, mlobs)
! dummy
type(ListType), intent(inout) :: list
@@ -112,7 +112,7 @@ subroutine AddMLObsToList(list, mlobs)
!
return
end subroutine AddMLObsToList
-
+
function GetMLObsFromList(list, indx) result(res)
! dummy
type(ListType), intent(inout) :: list
@@ -129,37 +129,37 @@ function GetMLObsFromList(list, indx) result(res)
!
return
end function GetMLObsFromList
-
+
! Type-bound procedures of MLObsType
-
+
subroutine CheckWeightSum(this)
! dummy
class(MLObsType) :: this
! local
double precision :: weightsum
integer :: i, nlayers
- type(LayerObsType), pointer :: layobs => null()
+ type(LayerObsType), pointer :: layobs => null()
character(len=MAXCHARLEN) :: ermsg
! formats
- 10 format('Weights of layer observations do not sum to 1.0 for', &
- ' multilayer observation: ',a)
+10 format('Weights of layer observations do not sum to 1.0 for', &
+ ' multilayer observation: ', a)
!
if (this%summ) return
!
weightsum = 0.0d0
nlayers = this%LayerObsList%Count()
- do i=1,nlayers
+ do i = 1, nlayers
layobs => GetLayerObsFromList(this%LayerObsList, i)
weightsum = weightsum + layobs%weight
- enddo
+ end do
!
- if (.not. is_same(weightsum, DONE)) then
- write(ermsg,10)trim(this%mlobsname)
+ if (.not. is_close(weightsum, DONE)) then
+ write (ermsg, 10) trim(this%mlobsname)
call store_error(ermsg)
call ustop()
- endif
+ end if
!
return
end subroutine
-
+
end module MultiLayerObs
diff --git a/utils/mf5to6/src/Preproc/ObsBlock.f90 b/utils/mf5to6/src/Preproc/ObsBlock.f90
index 61f3e5bad78..c4e51a3562e 100644
--- a/utils/mf5to6/src/Preproc/ObsBlock.f90
+++ b/utils/mf5to6/src/Preproc/ObsBlock.f90
@@ -3,7 +3,7 @@ module ObsBlockModule
use BlockParserModule, only: BlockParserType
use ConstantsModule, only: DONE, DZERO, &
LINELENGTH, MAXCHARLEN, LENOBSNAME
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use ConstantsPHMFModule, only: CONTINUOUS, SINGLE, LENOBSNAMENEW
use DnmDis3dModule, only: Dis3dType
use GlobalVariablesPHMFModule, only: verbose
@@ -184,7 +184,7 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser)
iadjrow = 0
jadjcol = 0
!
- if (.not. is_same(xoff, DZERO)) then
+ if (.not. is_close(xoff, DZERO)) then
if (xoff > DZERO) then
if (jcol < ncol) then
if (dis3d%idomain(jcol+1, irow, layer) == 1) then
@@ -202,7 +202,7 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser)
endif
endif
!
- if (.not. is_same(yoff, DZERO)) then
+ if (.not. is_close(yoff, DZERO)) then
if (yoff > DZERO) then
if (irow > 1) then
if (dis3d%idomain(jcol, irow-1, layer) == 1) then
diff --git a/utils/mf5to6/src/Preproc/ObservePHMF.f90 b/utils/mf5to6/src/Preproc/ObservePHMF.f90
index d35fd18079d..c50c7ad5f8d 100644
--- a/utils/mf5to6/src/Preproc/ObservePHMF.f90
+++ b/utils/mf5to6/src/Preproc/ObservePHMF.f90
@@ -15,7 +15,7 @@ module ObserveModule
use ConstantsModule, only: DONE, DZERO, LENOBSNAME, &
LENOBSTYPE, MAXCHARLEN
use ConstantsPHMFModule, only: LENOBSNAMENEW, HUGEDBL, HDRYDEFAULT
- use GenericUtilitiesModule, only: is_same
+ use MathUtilModule, only: is_close
use ListModule, only: ListType
use SimModule, only: store_warning, store_error, &
store_error_unit, ustop
@@ -206,7 +206,7 @@ subroutine CalcSimVal(this, itime)
sumweights = DZERO
k = 0
do i=1,nsrc
- if (is_same(this%srcvals(itime, i), this%hdry)) then
+ if (is_close(this%srcvals(itime, i), this%hdry)) then
k = k + 1
weights(i) = DZERO
else
diff --git a/utils/zonebudget/make/makefile b/utils/zonebudget/make/makefile
index db96d3c3471..0a848b02b10 100644
--- a/utils/zonebudget/make/makefile
+++ b/utils/zonebudget/make/makefile
@@ -25,7 +25,7 @@ $(OBJDIR)/compilerversion.o \
$(OBJDIR)/version.o \
$(OBJDIR)/Sim.o \
$(OBJDIR)/OpenSpec.o \
-$(OBJDIR)/genericutils.o \
+$(OBJDIR)/MathUtil.o \
$(OBJDIR)/InputOutput.o \
$(OBJDIR)/LongLineReader.o \
$(OBJDIR)/DevFeature.o \
diff --git a/utils/zonebudget/msvs/zonebudget.vfproj b/utils/zonebudget/msvs/zonebudget.vfproj
index f24872d4846..f6e8613efd2 100644
--- a/utils/zonebudget/msvs/zonebudget.vfproj
+++ b/utils/zonebudget/msvs/zonebudget.vfproj
@@ -44,7 +44,6 @@
-
diff --git a/utils/zonebudget/pymake/extrafiles.txt b/utils/zonebudget/pymake/extrafiles.txt
index 846b560f0d5..54d11147ff3 100644
--- a/utils/zonebudget/pymake/extrafiles.txt
+++ b/utils/zonebudget/pymake/extrafiles.txt
@@ -4,9 +4,9 @@
../../../src/Utilities/Budget.f90
../../../src/Utilities/Constants.f90
../../../src/Utilities/compilerversion.F90
-../../../src/Utilities/genericutils.f90
../../../src/Utilities/ErrorUtil.f90
../../../src/Utilities/GeomUtil.f90
+../../../src/Utilities/MathUtil.f90
../../../src/Utilities/InputOutput.f90
../../../src/Utilities/kind.f90
../../../src/Utilities/LongLineReader.f90