fix(UzfCellGroup): SQUARE_GWET option records ET as inflow instead of…

… outflow (#1951) * fix(UzfCellGroup): SQUARE_GWET option records ET as inflow instead of outflow * Add new test and add a release note * minor cleanup to keep in line with standard conventions * A more straight forward fix * supplement the autotest a bit more * transfer nlin function to UzfEtUtil.f90 * try unit testing with the SQUARE_GWET fix * fprettify and remove unused variables * took out 1 to many variables in previous commit * comma not necessary (debugging on CI) * further augmenting square_gwet unit tests * Update doc/ReleaseNotes/develop.tex Co-authored-by: langevin-usgs <[email protected]> --------- Co-authored-by: langevin-usgs <[email protected]>
MODFLOW-USGS · Nov 19, 2024 · 68f7f87 · 68f7f87
1 parent 8d839c8
commit 68f7f87
Show file tree

Hide file tree

Showing 5 changed files with 730 additions and 38 deletions.
diff --git a/autotest/TestUzfEtUtil.f90 b/autotest/TestUzfEtUtil.f90
@@ -1,11 +1,12 @@
 module TestUzfEtUtil
 
   use KindModule, only: I4B, DP
-  use ConstantsModule, only: DZERO, DONE, DTWO, DEM1, DHALF, D1P1, DSIX
+  use ConstantsModule, only: DZERO, DONE, DTWO, DEM1, DEM3, DHALF, D1P1, &
+                             DFIVETHIRDS, DSIX
   use MathUtilModule, only: is_close
   use testdrive, only: check, error_type, new_unittest, test_failed, &
                        to_string, unittest_type
-  use UzfETUtilModule, only: etfunc_lin, calc_lin_scaling_fac
+  use UzfETUtilModule, only: etfunc_lin, etfunc_nlin, calc_lin_scaling_fac
 
   implicit none
   private
@@ -18,7 +19,8 @@ subroutine collect_uzfetutil(testsuite)
     testsuite = [ &
                 new_unittest("etfunc_lin", test_etfunc_lin), &
                 new_unittest("calc_lin_scaling_fac", &
-                             test_calc_lin_scaling_fac) &
+                             test_calc_lin_scaling_fac), &
+                new_unittest("etfunc_nlin", test_etfunc_nlin) &
                 ]
   end subroutine collect_uzfetutil
 
@@ -107,4 +109,64 @@ subroutine test_etfunc_lin(error)
     if (allocated(error)) return
   end subroutine test_etfunc_lin
 
+  subroutine test_etfunc_nlin(error)
+    type(error_type), allocatable, intent(out) :: error
+    real(DP) :: rate1, rate2 !< calculated pET rates
+    real(DP) :: range !< value used for smoothing bottom of square_gwet interval
+    real(DP) :: atol !< user specified tolerance for comparing calculated values
+    ! local
+    real(DP) :: deriv_et !< derivative of gw ET for Newton addition to equations in _fn()
+    real(DP) :: extdp !< extinction depth
+    real(DP) :: hgwf !< groundwater head
+    real(DP) :: pET !< potential evapotranspiration
+    real(DP) :: trhs !< total uzf rhs contribution to GWF model
+    real(DP) :: thcof !< total uzf hcof contribution to GWF model
+    real(DP) :: celtop !< elevation of the top of the cell
+
+    ! water table exactly in the middle of the extinction depth, should return pET
+    atol = 1d-12
+    deriv_et = DZERO
+    extdp = DONE
+    pET = DEM1
+    trhs = DZERO
+    thcof = DZERO
+    celtop = DTWO
+    hgwf = 1.5_DP
+    deriv_et = DZERO
+    trhs = DZERO
+    thcof = DZERO
+
+    rate1 = etfunc_nlin(celtop, extdp, pET, deriv_et, trhs, thcof, hgwf)
+    call check(error, is_close(rate1, pET))
+    if (allocated(error)) return
+
+    ! similar to previous check, but testing a different groundwater depth
+    hgwf = DFIVETHIRDS
+    deriv_et = DZERO
+    trhs = DZERO
+    thcof = DZERO
+    rate2 = etfunc_nlin(celtop, extdp, pET, deriv_et, trhs, thcof, hgwf)
+    call check(error, is_close(rate1, rate2))
+    if (allocated(error)) return
+
+    ! even with the SQUARE_GWET function, there is smoothing close to the extinction depth
+    ! the following tests that the smoothing near the EXTDP scales pET by half
+    range = DEM3 * extdp ! an intermediate calc
+    hgwf = celtop - extdp + (range * DHALF)
+    rate1 = etfunc_nlin(celtop, extdp, pET, deriv_et, trhs, thcof, hgwf)
+    call check(error, is_close(rate1, pET * DHALF, atol=atol))
+    ! however, if water table is above (or in this case, at) the calculated "range",
+    ! there should be no scaling
+    hgwf = celtop - extdp + (range * DONE)
+    rate2 = etfunc_nlin(celtop, extdp, pET, deriv_et, trhs, thcof, hgwf)
+    call check(error, is_close(rate2, rate1 * DTWO, atol=atol))
+
+    ! when water table is at the extinction depth, scaling of gwet should result
+    ! in no gwet
+    hgwf = celtop - extdp
+    rate1 = etfunc_nlin(celtop, extdp, pET, deriv_et, trhs, thcof, hgwf)
+    call check(error, is_close(rate1, DZERO, atol=atol))
+
+  end subroutine test_etfunc_nlin
+
 end module TestUzfEtUtil