+Add cuberoot function

  Added the new functions cuberoot and intrinsic_functions_unit_tests to the
MOM_intrinsic_functions module, and call intrinsic_functions_unit_tests from
unit_tests to confirm that this new function works as intended.  Separately,
cuberoot was tested by replacing expressions like A**(1./3.) with cuberoot(A) in
MOM_energetic_PBL and verifying that the answers only change at roundoff, but
that it can give bitwise identical results when the argument is scaled by an
integer power of 8 and then unscaled by the corresponding integer power of 2,
but that change will occur in a subsequent commit as it can change answers
depending on an ANSWER_DATE flag.  With this commit, cuberoot is not yet being
used so all answers are bitwise identical, although there are new publicly
visible routines.

src/core/MOM_unit_tests.F90

@@ -4,15 +4,16 @@ module MOM_unit_tests
 ! This file is part of MOM6. See LICENSE.md for the license.
 
 use MOM_error_handler,              only : MOM_error, FATAL, is_root_pe
-
-use MOM_string_functions,           only : string_functions_unit_tests
-use MOM_remapping,                  only : remapping_unit_tests
+use MOM_hor_bnd_diffusion,          only : near_boundary_unit_tests
+use MOM_intrinsic_functions,        only : intrinsic_functions_unit_tests
+use MOM_mixed_layer_restrat,        only : mixedlayer_restrat_unit_tests
 use MOM_neutral_diffusion,          only : neutral_diffusion_unit_tests
 use MOM_random,                     only : random_unit_tests
-use MOM_hor_bnd_diffusion,          only : near_boundary_unit_tests
+use MOM_remapping,                  only : remapping_unit_tests
+use MOM_string_functions,           only : string_functions_unit_tests
 use MOM_CFC_cap,                    only : CFC_cap_unit_tests
 use MOM_EOS,                        only : EOS_unit_tests
-use MOM_mixed_layer_restrat,        only : mixedlayer_restrat_unit_tests
+
 implicit none ; private
 
 public unit_tests
@@ -36,6 +37,8 @@ subroutine unit_tests(verbosity)
        "MOM_unit_tests: EOS_unit_tests FAILED")
     if (remapping_unit_tests(verbose)) call MOM_error(FATAL, &
        "MOM_unit_tests: remapping_unit_tests FAILED")
+    if (intrinsic_functions_unit_tests(verbose)) call MOM_error(FATAL, &
+       "MOM_unit_tests: intrinsic_functions_unit_tests FAILED")
     if (neutral_diffusion_unit_tests(verbose)) call MOM_error(FATAL, &
        "MOM_unit_tests: neutralDiffusionUnitTests FAILED")
     if (random_unit_tests(verbose)) call MOM_error(FATAL, &

src/framework/MOM_intrinsic_functions.F90

@@ -4,9 +4,12 @@ module MOM_intrinsic_functions
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit
+
 implicit none ; private
 
-public :: invcosh
+public :: invcosh, cuberoot
+public :: intrinsic_functions_unit_tests
 
 contains
 
@@ -25,4 +28,137 @@ function invcosh(x)
 
 end function invcosh
 
+!> Returns the cube root of a real argument at roundoff accuracy, in a form that works properly with
+!! rescaling of the argument by integer powers of 8.  If the argument is a NaN, a NaN is returned.
+elemental function cuberoot(x) result(root)
+  real, intent(in) :: x !< The argument of cuberoot in arbitrary units cubed [A3]
+  real :: root !< The real cube root of x in arbitrary units [A]
+
+  real :: asx ! The absolute value of x rescaled by an integer power of 8 to put it into
+              ! the range from 0.125 < asx <= 1.0, in ambiguous units cubed [B3]
+  real :: root_asx ! The cube root of asx [B]
+  real :: num ! The numerator of an expression for the evolving estimate of the cube root of asx
+              ! in arbitrary units that can grow or shrink with each iteration [B C]
+  real :: den ! The denominator of an expression for the evolving estimate of the cube root of asx
+              ! in arbitrary units that can grow or shrink with each iteration [C]
+  real :: num_prev ! The numerator of an expression for the previous iteration of the evolving estimate
+              ! of the cube root of asx in arbitrary units that can grow or shrink with each iteration [B D]
+  real :: den_prev ! The denominator of an expression for the previous iteration of the evolving estimate of
+              ! the cube root of asx in arbitrary units that can grow or shrink with each iteration [D]
+  integer :: ex_3 ! One third of the exponent part of x, used to rescale x to get a.
+  integer :: itt
+
+  if ((x >= 0.0) .eqv. (x <= 0.0)) then
+    ! Return 0 for an input of 0, or NaN for a NaN input.
+    root = x
+  else
+    ex_3 = exponent(x) / 3
+    asx = scale(abs(x), -3*ex_3)
+    if (asx > 1.0) then
+      asx = 0.125 * asx ; ex_3 = ex_3 + 1
+    elseif (asx <= 0.125) then
+      ! This is retained in case ((3.0 / 3.0) == 1.0) is not exactly true, which might lead
+      ! to convergence to slightly different values for root for asx = 0.125 or asx = 1.0.
+      asx = 8.0 * asx ; ex_3 = ex_3 - 1
+    endif
+
+    if (asx == 1.0) then
+      root_asx = 1.0
+    else
+      ! This first estimate is one iteration of Newton's method with a starting guess of 1.  It is
+      ! always an over-estimate of the true solution, but it is a good approximation for asx near 1.
+      num = (2.0 + asx)
+      den = 3.0
+      ! Iteratively determine Root = asx**1/3 using Newton's method, noting that in this case Newton's
+      ! method converges monotonically from above and needs no bounding.  For the range of asx from
+      ! 0.125 to 1.0 with the first guess used above, 6 iterations suffice to converge to roundoff.
+      do itt=1,9
+        ! Newton's method iterates estimates as Root = Root - (Root**3 - asx) / (3.0 * Root**2), or
+        ! equivalently as Root = (2.0*Root**2 + asx) / (3.0 * Root**2).
+        ! Keeping the estimates in a fractional form Root = num / den allows this calculation with
+        ! fewer (or no) real divisions during the iterations before doing a single real division
+        ! at the end, and it is therefore more computationally efficient.
+
+        ! Because successive estimates of the numerator and denominator tend to be the cube of their
+        ! predecessors, the numerator and denominator need to be rescaled by division when they get
+        ! too large or small to avoid overflow or underflow in the convergence test below.  Were
+        ! this being done for 32 bit reals, the values to compare with would be about 1.0e9 and 1.0e-9.
+        if ((den > 1.0e75) .or. (den < 1.0e-75)) then
+          num = num / den ; den = 1.0
+        endif
+
+        num_prev = num ; den_prev = den
+        num = (2.0 * num_prev**3 + asx * den_prev**3)
+        den = 3.0 * (den_prev * num_prev**2)
+
+        !   If successive estimates of root are identical, this is an exact or converged solution.
+        if ((num * den_prev == num_prev * den) .or. (itt == 9)) then
+          root_asx = num / den
+          exit
+        endif
+
+        !   If the estimates are increasing, this also indicates convergence, but for a more subtle
+        ! reason.  Because Newton's method converges monotonically from above (at least for infinite
+        ! precision math), the only reason why this estimate could increase is if the iterations
+        ! have converged to a roundoff-level limit cycle around an irrational or otherwise
+        ! unrepresentable solution, with values only changing in the last bit or two.  If so, we
+        ! should stop iterating and accept the one of the current or previous solutions, both of
+        ! which will be within numerical roundoff of the true solution.
+        if (num * den_prev > num_prev * den) then
+          root_asx = num / den
+          ! Pick the more accurate of the last two iterations.
+          ! Given that both of the two previous iterations are within roundoff of the true
+          ! solution, this next step might be overkill.
+          if ( abs(den_prev**3*root_asx**3 - den_prev**3*asx) > abs(num_prev**3 - den_prev**3*asx) ) then
+            ! The previous iteration was slightly more accurate, so use that for root_asx.
+            root_asx = num_prev / den_prev
+          endif
+          exit
+        endif
+      enddo
+    endif
+
+    root = sign(scale(root_asx, ex_3), x)
+  endif
+
+end function cuberoot
+
+!> Returns true if any unit test of intrinsic_functions fails, or false if they all pass.
+logical function intrinsic_functions_unit_tests(verbose)
+  logical, intent(in) :: verbose !< If true, write results to stdout
+
+  ! Local variables
+  real :: testval  ! A test value for self-consistency testing [nondim]
+  logical :: fail, v
+
+  fail = .false.
+  v = verbose
+  write(stdout,*) '==== MOM_intrinsic_functions: intrinsic_functions_unit_tests ==='
+
+  fail = fail .or. Test_cuberoot(v, 1.2345678901234e9)
+  fail = fail .or. Test_cuberoot(v, -9.8765432109876e-21)
+  fail = fail .or. Test_cuberoot(v, 64.0)
+  fail = fail .or. Test_cuberoot(v, -0.5000000000001)
+  fail = fail .or. Test_cuberoot(v, 0.0)
+
+end function intrinsic_functions_unit_tests
+
+!> True if the cube of cuberoot(val) does not closely match val. False otherwise.
+logical function Test_cuberoot(verbose, val)
+  logical, intent(in) :: verbose !< If true, write results to stdout
+  real, intent(in) :: val  !< The real value to test, in arbitrary units [A]
+  ! Local variables
+  real :: diff ! The difference between val and the cube root of its cube.
+
+  diff = val - cuberoot(val**3)
+  Test_cuberoot = (abs(diff) > 2.0e-15*abs(val))
+
+  if (Test_cuberoot) then
+    write(stdout, '("For val = ",ES22.15,", (val - cuberoot(val**3))) = ",ES9.2," <-- FAIL")') val, diff
+  elseif (verbose) then
+    write(stdout, '("For val = ",ES22.15,", (val - cuberoot(val**3))) = ",ES9.2)') val, diff
+
+  endif
+end function Test_cuberoot
+
 end module MOM_intrinsic_functions