feat(MathUtil): add generic routine for modulo with offset (#1494)

autotest/TestMathUtil.f90

@@ -0,0 +1,44 @@
+module TestMathUtil
+  use KindModule, only: I4B, DP
+  use testdrive, only: check, error_type, new_unittest, test_failed, &
+                       to_string, unittest_type
+  use MathUtilModule, only: mod_offset
+  implicit none
+  private
+  public :: collect_mathutil
+
+contains
+
+  subroutine collect_mathutil(testsuite)
+    type(unittest_type), allocatable, intent(out) :: testsuite(:)
+    testsuite = [ &
+                new_unittest("mod_offset", &
+                             test_mod_offset) &
+                ]
+  end subroutine collect_mathutil
+
+  subroutine test_mod_offset(error)
+    type(error_type), allocatable, intent(out) :: error
+
+    ! with no offset specified, should behave just like mod
+    call check(error, mod_offset(2, 2) == 0)
+    call check(error, mod_offset(2, 3) == 2)
+    call check(error, mod_offset(2.0_DP, 2.0_DP) == 0.0_DP)
+    call check(error, mod_offset(2.0_DP, 3.0_DP) == 2.0_DP)
+
+    ! with offset d specified, if the result x = a mod n falls
+    ! between 0 and n - 1, the new result x = a mod_d n falls
+    ! between d and d + n - 1.
+    call check(error, mod_offset(2, 3, -2) == -1)
+    call check(error, mod_offset(2, 3, -1) == -1)
+    call check(error, mod_offset(2, 3, 0) == 2)
+    call check(error, mod_offset(2, 3, 1) == 2)
+    call check(error, mod_offset(2, 3, 2) == 2)
+    call check(error, mod_offset(2, 3, 3) == 5)
+    call check(error, mod_offset(2, 3, 4) == 5)
+    call check(error, mod_offset(2.0_DP, 3.0_DP, -1.0_DP) == -1.0_DP)
+    call check(error, mod_offset(2.0_DP, 3.0_DP, 2.0_DP) == 2.0_DP)
+    call check(error, mod_offset(2.0_DP, 3.0_DP, 3.0_DP) == 5.0_DP)
+  end subroutine test_mod_offset
+
+end module TestMathUtil

autotest/meson.build

@@ -5,6 +5,7 @@ if test_drive.found() and not fc_id.contains('intel')
     'DevFeature',
     'GeomUtil',
     'InputOutput',
+    'MathUtil',
     'Sim'
   ]
 

autotest/tester.f90

@@ -6,6 +6,7 @@ program tester
   use TestDevFeature, only: collect_dev_feature
   use TestGeomUtil, only: collect_geomutil
   use TestInputOutput, only: collect_inputoutput
+  use TestMathUtil, only: collect_mathutil
   use TestSim, only: collect_sim
   implicit none
   integer :: stat, is
@@ -19,6 +20,7 @@ program tester
                new_testsuite("DevFeature", collect_dev_feature), &
                new_testsuite("GeomUtil", collect_geomutil), &
                new_testsuite("InputOutput", collect_inputoutput), &
+               new_testsuite("MathUtil", collect_mathutil), &
                new_testsuite("Sim", collect_sim) &
                ]
 

make/makefile

@@ -86,6 +86,7 @@ $(OBJDIR)/version.o \
 $(OBJDIR)/Message.o \
 $(OBJDIR)/Sim.o \
 $(OBJDIR)/OpenSpec.o \
+$(OBJDIR)/MathUtil.o \
 $(OBJDIR)/InputOutput.o \
 $(OBJDIR)/TableTerm.o \
 $(OBJDIR)/Table.o \

msvs/mf6core.vfproj

@@ -395,6 +395,7 @@
 		<File RelativePath="..\src\Utilities\List.f90"/>
 		<File RelativePath="..\src\Utilities\ListReader.f90"/>
 		<File RelativePath="..\src\Utilities\LongLineReader.f90"/>
+		<File RelativePath="..\src\Utilities\MathUtil.f90"/>
 		<File RelativePath="..\src\Utilities\Message.f90"/>
 		<File RelativePath="..\src\Utilities\OpenSpec.f90"/>
 		<File RelativePath="..\src\Utilities\PackageBudget.f90"/>

src/Utilities/MathUtil.f90

@@ -0,0 +1,54 @@
+module MathUtilModule
+  use KindModule, only: DP, I4B, LGP
+  use ErrorUtilModule, only: pstop
+  use ConstantsModule, only: MAXCHARLEN, LENHUGELINE, &
+                             DZERO, DPREC, DSAME, &
+                             LINELENGTH, LENHUGELINE, VSUMMARY
+
+  implicit none
+  private
+  public :: mod_offset
+
+  interface mod_offset
+    module procedure :: mod_offset_int, mod_offset_dbl
+  end interface mod_offset
+
+contains
+
+  !> @brief Modulo with offset for integer values.
+  pure function mod_offset_int(a, n, d) result(mo)
+    ! -- dummy
+    integer(I4B), intent(in) :: a !< dividend
+    integer(I4B), intent(in) :: n !< divisor
+    integer(I4B), intent(in), optional :: d !< offset
+    integer(I4B) :: mo
+    ! -- local
+    integer(I4B) :: ld
+
+    if (present(d)) then
+      ld = d
+    else
+      ld = 0
+    end if
+    mo = a - n * floor(real(a - ld) / n)
+  end function mod_offset_int
+
+  !> @brief Modulo with offset for double precision values.
+  pure function mod_offset_dbl(a, n, d) result(mo)
+    ! -- dummy
+    real(DP), intent(in) :: a !< dividend
+    real(DP), intent(in) :: n !< divisor
+    real(DP), intent(in), optional :: d !< offset
+    real(DP) :: mo
+    ! -- local
+    real(DP) :: ld
+
+    if (present(d)) then
+      ld = d
+    else
+      ld = 0
+    end if
+    mo = a - n * floor((a - ld) / n)
+  end function mod_offset_dbl
+
+end module MathUtilModule

src/meson.build

@@ -237,6 +237,7 @@ modflow_sources = files(
     'Utilities' / 'List.f90',
     'Utilities' / 'ListReader.f90',
     'Utilities' / 'LongLineReader.f90',
+    'Utilities' / 'MathUtil.f90',
     'Utilities' / 'Message.f90',
     'Utilities' / 'OpenSpec.f90',
     'Utilities' / 'PackageBudget.f90',