+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,103 @@ 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 < a <= 1.0, in ambiguous units cubed [B3]
+  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 :: np2 ! The previous value of num squared [B2 C2]
+  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
+      asx = 8.0 * asx
+    endif
+
+    num = (2.0 + asx)
+    den = 3.0
+    if (asx < 1.0) then
+      ! Iteratively determine R = 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 a first guess of 1.0, 6 iterations suffice to converge to roundoff.
+      do itt=1,6
+        ! Newton's method iterates estimates as Root = Root - (Root**3 - 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 efficienty.
+
+        ! 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.  Were this being done for 32 bit reals,
+        ! the values to compare with would be about 1.0e12 and 1.0e-12.
+        if ((den > 1.0e100) .or. (den < 1.0e-100)) then
+          num = num / den ; den = 1.0
+        endif
+
+        ! Test whether this iteration would change anything.
+        if (abs(num**3 - asx * den**3) < 5.e-16 * num**3) exit
+
+        np2 = num**2
+        num = (2.0 * num**3 + asx * den**3)
+        den = 3.0 * (den * np2)
+      enddo
+    endif
+    root = sign(scale(num / den, 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