Merge branch 'dev/gfdl' into convert_temp_salt_for_TEOS10_bugfix

docs/forcing.rst

@@ -1,5 +1,33 @@
 Forcing
 =======
+Data Override
+-------
+When running MOM6 with the Flexible Modelling System (FMS) coupler, forcing can be specified by a `data_table` file. This is particularly useful when running MOM6 with a data atmosphere, as paths to the relevent atmospheric forcing products (eg. JRA55-do or ERA5) can be provided here. Each item in the data table must be separated by a new line, and contains the following information:
+
+| ``gridname``: The component of the model this data applies to. eg. `atm` `ocn` `lnd` `ice`.
+| ``fieldname_code``: The field name according to the model component. eg. `salt`
+| ``fieldname_file``: The name of the field within the source file. 
+| ``file_name``: Path to the source file.
+| ``interpol_method``: Interpolation method eg. `bilinear`
+| ``factor``: A scalar by which to multiply the field ahead of passing it onto the model. This is a quick way to do unit conversions for example. 
+
+| 
+The data table is commonly formatted by specifying each of the fields in the order listed above, with a new line for each entry.
+
+Example Format:
+    "ATM", "t_bot",  "t2m", "./INPUT/2t_ERA5.nc", "bilinear", 1.0
+
+A `yaml` format is also possible if you prefer. This is outlined in the `FMS data override <https://github.com/NOAA-GFDL/FMS/tree/main/data_override>`_ github page, along with other details. 
+
+Speficying a constant value:
+    Rather than overriding with data from a file, one can also set a field to constant. To do this, pass empty strings to `fieldname_file` and `file_name`. The `factor` now corresponds to the override value. For example, the following sets the temperature at the bottom of the atmosphere to 290 Kelvin.
+
+
+    "ATM", "t_bot",  "", "", "bilinear", 290.0
+
+Which units do I need?
+    For configurations using SIS2 and MOM, a list of available surface flux variables along with the expected units can be found in the `flux_exchange <https://github.com/NOAA-GFDL/FMScoupler/blob/f4782c2c033df086eeac29fbbefb4a0bdac1649f/full/flux_exchange.F90>`_ file. 
+
 
 .. toctree::
     :maxdepth: 2

src/diagnostics/MOM_spatial_means.F90

@@ -211,11 +211,13 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale)
   ! Local variables
   ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
   ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
-  real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: tmpForSumming  ! An unscaled cell integral [a m3]
-  real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: weight  ! The volume of each cell, used as a weight [m3]
+  real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: tmpForSumming  ! An unscaled cell integral [a m3] or [a kg]
+  real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: weight  ! The volume or mass of each cell, depending on
+                                                    ! whether the model is Boussinesq, used as a weight [m3] or [kg]
   type(EFP_type), dimension(2*SZK_(GV)) :: laysums
-  real, dimension(SZK_(GV)) :: global_temp_scalar    ! The global integral of the tracer in each layer [a m3]
-  real, dimension(SZK_(GV)) :: global_weight_scalar  ! The global integral of the volume of each layer [m3]
+  real, dimension(SZK_(GV)) :: global_temp_scalar   ! The global integral of the tracer in each layer [a m3] or [a kg]
+  real, dimension(SZK_(GV)) :: global_weight_scalar ! The global integral of the volume or mass of each
+                                                    ! layer [m3] or [kg]
   real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
   real :: scalefac  ! A scaling factor for the variable [a A-1 ~> 1]
   integer :: i, j, k, is, ie, js, je, nz
@@ -226,7 +228,7 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale)
   tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0.
 
   do k=1,nz ; do j=js,je ; do i=is,ie
-    weight(i,j,k)  =  (GV%H_to_m * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
+    weight(i,j,k)  =  (GV%H_to_MKS * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
     tmpForSumming(i,j,k) =  scalefac * var(i,j,k) * weight(i,j,k)
   enddo ; enddo ; enddo
 
@@ -262,9 +264,9 @@ function global_volume_mean(var, h, G, GV, scale, tmp_scale)
   ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
   real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
   real :: scalefac   ! A scaling factor for the variable [a A-1 ~> 1]
-  real :: weight_here ! The volume of a grid cell [m3]
-  real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! The volume integral of the variable in a column [a m3]
-  real, dimension(SZI_(G),SZJ_(G)) :: sum_weight  ! The volume of each column of water [m3]
+  real :: weight_here ! The volume or mass of a grid cell [m3] or [kg]
+  real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! The volume integral of the variable in a column [a m3] or [a kg]
+  real, dimension(SZI_(G),SZJ_(G)) :: sum_weight  ! The volume or mass of each column of water [m3] or [kg]
   integer :: i, j, k, is, ie, js, je, nz
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
@@ -273,7 +275,7 @@ function global_volume_mean(var, h, G, GV, scale, tmp_scale)
   tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0.
 
   do k=1,nz ; do j=js,je ; do i=is,ie
-    weight_here  =  (GV%H_to_m * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
+    weight_here  =  (GV%H_to_MKS * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
     tmpForSumming(i,j) = tmpForSumming(i,j) + scalefac * var(i,j,k) * weight_here
     sum_weight(i,j) = sum_weight(i,j) + weight_here
   enddo ; enddo ; enddo

src/framework/MOM_diag_remap.F90

@@ -658,8 +658,15 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
   logical, dimension(:), intent(inout) :: averaged_mask  !< Mask for horizontally averaged field [nondim]
 
   ! Local variables
-  real, dimension(G%isc:G%iec, G%jsc:G%jec, size(field,3)) :: volume, stuff
-  real, dimension(size(field, 3)) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages
+  real :: volume(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area [m2], volume [m3] or mass [kg] of each cell.
+  real :: stuff(G%isc:G%iec, G%jsc:G%jec, size(field,3))  ! The area, volume or mass-weighted integral of the
+                                             ! field being averaged in each cell, in [m2 A], [m3 A] or [kg A],
+                                             ! depending on the weighting for the averages and whether the
+                                             ! model makes the Boussinesq approximation.
+  real, dimension(size(field, 3)) :: vol_sum   ! The global sum of the areas [m2], volumes [m3] or mass [kg]
+                                               ! in the cells that used in the weighted averages.
+  real, dimension(size(field, 3)) :: stuff_sum ! The global sum of the weighted field in all cells, in
+                                               ! [A m2], [A m3] or [A kg]
   type(EFP_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer
   real :: height  ! An average thickness attributed to an velocity point [H ~> m or kg m-2]
   integer :: i, j, k, nz
@@ -688,7 +695,7 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
             I1 = i - G%isdB + 1
             height = 0.5 * (h(i,j,k) + h(i+1,j,k))
             volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) &
-                * (GV%H_to_m * height) * G%mask2dCu(I,j)
+                * (GV%H_to_MKS * height) * G%mask2dCu(I,j)
             stuff(I,j,k) = volume(I,j,k) * field(I1,j,k)
           enddo ; enddo
         endif
@@ -717,7 +724,7 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
             J1 = J - G%jsdB + 1
             height = 0.5 * (h(i,j,k) + h(i,j+1,k))
             volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) &
-                * (GV%H_to_m * height) * G%mask2dCv(i,J)
+                * (GV%H_to_MKS * height) * G%mask2dCv(i,J)
             stuff(i,J,k) = volume(i,J,k) * field(i,J1,k)
           enddo ; enddo
         endif
@@ -748,7 +755,7 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
         else ! Intensive
           do j=G%jsc, G%jec ; do i=G%isc, G%iec
             volume(i,j,k) = (G%US%L_to_m**2 * G%areaT(i,j)) &
-                * (GV%H_to_m * h(i,j,k)) * G%mask2dT(i,j)
+                * (GV%H_to_MKS * h(i,j,k)) * G%mask2dT(i,j)
             stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
           enddo ; enddo
         endif

src/framework/MOM_file_parser.F90

@@ -56,6 +56,8 @@ module MOM_file_parser
 !> Specify the active parameter block
 type, private :: parameter_block ; private
   character(len=240) :: name = ''   !< The active parameter block name
+  logical :: log_access = .true.
+    !< Log the entry and exit of the block (but not its contents)
 end type parameter_block
 
 !> A structure that can be parsed to read and document run-time parameters.
@@ -2082,17 +2084,29 @@ subroutine clearParameterBlock(CS)
 end subroutine clearParameterBlock
 
 !> Tags blockName onto the end of the active parameter block name
-subroutine openParameterBlock(CS,blockName,desc)
+subroutine openParameterBlock(CS, blockName, desc, do_not_log)
   type(param_file_type),      intent(in) :: CS      !< The control structure for the file_parser module,
                                          !! it is also a structure to parse for run-time parameters
   character(len=*),           intent(in) :: blockName !< The name of a parameter block being added
   character(len=*), optional, intent(in) :: desc    !< A description of the parameter block being added
+  logical, optional, intent(in) :: do_not_log
+    !< Log block entry if true.  This only prevents logging of entry to the block, and not the contents.
 
   type(parameter_block), pointer :: block => NULL()
+  logical :: do_log
+
+  do_log = .true.
+  if (present(do_not_log)) do_log = .not. do_not_log
+
   if (associated(CS%blockName)) then
     block => CS%blockName
     block%name = pushBlockLevel(block%name,blockName)
-    call doc_openBlock(CS%doc,block%name,desc)
+    if (do_log) then
+      call doc_openBlock(CS%doc, block%name, desc)
+      block%log_access = .true.
+    else
+      block%log_access = .false.
+    endif
   else
     if (is_root_pe()) call MOM_error(FATAL, &
       'openParameterBlock: A push was attempted before allocation.')
@@ -2111,7 +2125,7 @@ subroutine closeParameterBlock(CS)
     if (is_root_pe().and.len_trim(block%name)==0) call MOM_error(FATAL, &
       'closeParameterBlock: A pop was attempted on an empty stack. ("'//&
       trim(block%name)//'")')
-    call doc_closeBlock(CS%doc,block%name)
+    if (block%log_access) call doc_closeBlock(CS%doc, block%name)
   else
     if (is_root_pe()) call MOM_error(FATAL, &
       'closeParameterBlock: A pop was attempted before allocation.')

src/framework/MOM_intrinsic_functions.F90

@@ -5,12 +5,26 @@ 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
+use iso_fortran_env, only : int64, real64
 
 implicit none ; private
 
 public :: invcosh, cuberoot
 public :: intrinsic_functions_unit_tests
 
+! Floating point model, if bit layout from high to low is (sign, exp, frac)
+
+integer, parameter :: bias = maxexponent(1.) - 1
+  !< The double precision exponent offset
+integer, parameter :: signbit = storage_size(1.) - 1
+  !< Position of sign bit
+integer, parameter :: explen = 1 + ceiling(log(real(bias))/log(2.))
+  !< Bit size of exponent
+integer, parameter :: expbit = signbit - explen
+  !< Position of lowest exponent bit
+integer, parameter :: fraclen = expbit
+  !< Length of fractional part
+
 contains
 
 !> Evaluate the inverse cosh, either using a math library or an
@@ -28,6 +42,7 @@ 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)
@@ -37,24 +52,28 @@ elemental function cuberoot(x) result(root)
   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 :: ra_3 ! root_asx 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 :: 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 :: np_3 ! num_prev cubed  [B3 D3]
   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.
+  real :: dp_3 ! den_prev cubed  [C3]
+  real :: r0  ! Initial value of the iterative solver. [B C]
+  real :: r0_3 ! r0 cubed [B3 C3]
   integer :: itt
 
+  integer(kind=int64) :: e_x, s_x
+
   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 = ceiling(exponent(x) / 3.)
-    ! Here asx is in the range of 0.125 <= asx < 1.0
-    asx = scale(abs(x), -3*ex_3)
+    call rescale_cbrt(x, asx, e_x, s_x)
 
     !   Iteratively determine root_asx = asx**1/3 using Halley's method and then Newton's method,
     ! noting that Halley's method onverges monotonically and needs no bounding.  Halley's method is
@@ -65,35 +84,115 @@ elemental function cuberoot(x) result(root)
 
     ! This first estimate gives the same magnitude of errors for 0.125 and 1.0 after two iterations.
     ! The first iteration is applied explicitly.
-    num = 0.707106 * (0.707106**3 + 2.0 * asx)
-    den = 2.0 * (0.707106**3) + asx
+    r0 = 0.707106
+    r0_3 = r0 * r0 * r0
+    num = r0 * (r0_3 + 2.0 * asx)
+    den = 2.0 * r0_3 + asx
 
     do itt=1,2
       ! Halley's method iterates estimates as Root = Root * (Root**3 + 2.*asx) / (2.*Root**3 + asx).
       num_prev = num ; den_prev = den
-      num = num_prev * (num_prev**3 + 2.0 * asx * (den_prev**3))
-      den = den_prev * (2.0 * num_prev**3 + asx * (den_prev**3))
+
+      ! Pre-compute these as integer powers, to avoid `pow()`-like intrinsics.
+      np_3 = num_prev * num_prev * num_prev
+      dp_3 = den_prev * den_prev * den_prev
+
+      num = num_prev * (np_3 + 2.0 * asx * dp_3)
+      den = den_prev * (2.0 * np_3 + asx * dp_3)
       ! Equivalent to:  root_asx = root_asx * (root_asx**3 + 2.*asx) / (2.*root_asx**3 + asx)
     enddo
     ! At this point the error in root_asx is better than 1 part in 3e14.
     root_asx = num / den
 
     ! One final iteration with Newton's method polishes up the root and gives a solution
     ! that is within the last bit of the true solution.
-    root_asx = root_asx - (root_asx**3 - asx) / (3.0 * (root_asx**2))
+    ra_3 = root_asx * root_asx * root_asx
+    root_asx = root_asx - (ra_3 - asx) / (3.0 * (root_asx * root_asx))
 
-    root = sign(scale(root_asx, ex_3), x)
+    root = descale(root_asx, e_x, s_x)
   endif
-
 end function cuberoot
 
+
+!> Rescale `a` to the range [0.125, 1) and compute its cube-root exponent.
+pure subroutine rescale_cbrt(a, x, e_r, s_a)
+  real, intent(in) :: a
+    !< The real parameter to be rescaled for cube root
+  real, intent(out) :: x
+    !< The rescaled value of a
+  integer(kind=int64), intent(out) :: e_r
+    !< Cube root of the exponent of the rescaling of `a`
+  integer(kind=int64), intent(out) :: s_a
+    !< The sign bit of a
+
+  integer(kind=int64) :: xb
+    ! Floating point value of a, bit-packed as an integer
+  integer(kind=int64) :: e_a
+    ! Unscaled exponent of a
+  integer(kind=int64) :: e_x
+    ! Exponent of x
+  integer(kind=int64) :: e_div, e_mod
+    ! Quotient and remainder of e in e = 3*(e/3) + modulo(e,3).
+
+  ! Pack bits of a into xb and extract its exponent and sign.
+  xb = transfer(a, 1_int64)
+  s_a = ibits(xb, signbit, 1)
+  e_a = ibits(xb, expbit, explen) - bias
+
+  ! Compute terms of exponent decomposition e = 3*(e/3) + modulo(e,3).
+  ! (Fortran division is round-to-zero, so we must emulate floor division.)
+  e_mod = modulo(e_a, 3_int64)
+  e_div = (e_a - e_mod)/3
+
+  ! Our scaling decomposes e_a into e = {3*(e/3) + 3} + {modulo(e,3) - 3}.
+
+  ! The first term is a perfect cube, whose cube root is computed below.
+  e_r = e_div + 1
+
+  ! The second term ensures that x is shifted to [0.125, 1).
+  e_x = e_mod - 3
+
+  ! Insert the new 11-bit exponent into xb and write to x and extend the
+  ! bitcount to 12, so that the sign bit is zero and x is always positive.
+  call mvbits(e_x + bias, 0, explen + 1, xb, fraclen)
+  x = transfer(xb, 1.)
+end subroutine rescale_cbrt
+
+
+!> Undo the rescaling of a real number back to its original base.
+pure function descale(x, e_a, s_a) result(a)
+  real, intent(in) :: x
+    !< The rescaled value which is to be restored.
+  integer(kind=int64), intent(in) :: e_a
+    !< Exponent of the unscaled value
+  integer(kind=int64), intent(in) :: s_a
+    !< Sign bit of the unscaled value
+  real :: a
+    !< Restored value with the corrected exponent and sign
+
+  integer(kind=int64) :: xb
+    ! Bit-packed real number into integer form
+  integer(kind=int64) :: e_x
+    ! Biased exponent of x
+
+  ! Apply the corrected exponent and sign to x.
+  xb = transfer(x, 1_8)
+  e_x = ibits(xb, expbit, explen)
+  call mvbits(e_a + e_x, 0, explen, xb, expbit)
+  call mvbits(s_a, 0, 1, xb, signbit)
+  a = transfer(xb, 1.)
+end function descale
+
+
 !> Returns true if any unit test of intrinsic_functions fails, or false if they all pass.
-logical function intrinsic_functions_unit_tests(verbose)
+function intrinsic_functions_unit_tests(verbose) result(fail)
   logical, intent(in) :: verbose !< If true, write results to stdout
+  logical :: fail !< True if any of the unit tests fail
 
   ! Local variables
   real :: testval  ! A test value for self-consistency testing [nondim]
-  logical :: fail, v
+  logical :: v
+  integer :: n
 
   fail = .false.
   v = verbose
@@ -107,7 +206,15 @@ logical function intrinsic_functions_unit_tests(verbose)
   fail = fail .or. Test_cuberoot(v, 1.0)
   fail = fail .or. Test_cuberoot(v, 0.125)
   fail = fail .or. Test_cuberoot(v, 0.965)
-
+  fail = fail .or. Test_cuberoot(v, 1.0 - epsilon(1.0))
+  fail = fail .or. Test_cuberoot(v, 1.0 - 0.5*epsilon(1.0))
+
+  testval = 1.0e-99
+  v = .false.
+  do n=-160,160
+    fail = fail .or. Test_cuberoot(v, testval)
+    testval = (-2.908 * (1.414213562373 + 1.2345678901234e-5*n)) * testval
+  enddo
 end function intrinsic_functions_unit_tests
 
 !> True if the cube of cuberoot(val) does not closely match val. False otherwise.
@@ -117,7 +224,7 @@ logical function Test_cuberoot(verbose, val)
   ! Local variables
   real :: diff ! The difference between val and the cube root of its cube.
 
-  diff = val - cuberoot(val**3)
+  diff = val - cuberoot(val)**3
   Test_cuberoot = (abs(diff) > 2.0e-15*abs(val))
 
   if (Test_cuberoot) then