Merge in latest dev/gfdl updates (#40)

* (*)Fixed dimensional inconsistency in P3M_functions

  Corrected dimensionally inconsistent expressions in P3M_functions.F90,
notably in P3M_limiter and monotonize_cubic and a complete rewrite and
simplification of is_cubic_monotonic.  Also added comments documenting the
units of all real variables in this module, and changed the code to use logical
variables in place of integer "booleans", including in the return value from
is_cubic_monotonic.  These changes will change (fix) the answers when remapping
variables with small numerical values, but no answers change in the
MOM6-examples test cases.

* +Added REMAPPING_2018 runtime option

  Added a new runtime option, REMAPPING_2018, which if set to false triggers the
use of new, more accurate expressions in various parts of the ALE remapping
code.  By default, the older expressions are used, and all answers are bitwise
identical, but there are new optional arguments to various routines related to
remapping to trigger the use of new mathematically equivalent expressions.  By
default all answers are bitwise identical, but there are new and reordered
entries in the MOM6_parameter_doc files.

* Corrected the formatting of a doxygen comment

* Added conversion factors to forcing diagnostics

  Added conversion factors to 4 mass-flux diagnostics and comments to 4 others
on why no conversion factors are needed.  All answers are bitwise identical.

* Added correct scaling factors to chksum calls

  Added scale arguments to 5 chksum calls and grouped another two chksum calls
while also adding the right scaling argument. All answers are bitwise identical.

* +Unscales area before taking global sum

  Undoes the dimensional scaling of the cell areas before taking their global
sum, so that the reproducing sum does not overflow when there is dimensional
rescaling.  All answers are bitwise identical when there is no rescaling, but
this eliminates a source of inadvertent overflows or underflows in the global
sums, and there is a new optional argument to compute_global_grid_integrals.

* (*)Correct dimensionally inconsistent advective CFL

  Corrects the dimensionally inconsistent expressions for the CFL number in
the tracer advection code, in which a negligible thickness had been added to
the cell volume to avoid division by zero.  This change does not alter the
solutions in the MOM6-examples test cases, but now it permits dimensional
rescaling of lengths over a much larger range, and it could change answers if
the minimum layer thicknesses are small enough.

* Unscale sea level before averaging

  Unscale interface heights before taking a global average via a reproducing sum
in non-Boussinesq mode global diagnostics to permit dimensional consistency
testing over a larger range.  All answers are bitwise identical.

* +Added an optional tmp_scale arg to global_i_mean

  Added an optional tmp_scale argument to global_i_mean and global_j_mean to
specify an internal rescaling of variables being averaged before the reproducing
sum.  All answers are bitwise identical, but there are new optional arguments
to two public interfaces.

* Expand consistency testing with i-mean sponges

  Use tmp_scale when taking the i-mean interface heights for i-mean sponges, to
give a greatly expanded range of dimensional consistency testing.  All answers
are bitwise identical.

Author: wrongkindofdoctor

src/ALE/MOM_ALE.F90

@@ -63,19 +63,23 @@ module MOM_ALE
 
 !> ALE control structure
 type, public :: ALE_CS ; private
-  logical :: remap_uv_using_old_alg              !< If true, uses the old "remapping via a delta z"
-                                                 !! method. If False, uses the new method that
-                                                 !! remaps between grids described by h.
+  logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z"
+                                    !! method. If False, uses the new method that
+                                    !! remaps between grids described by h.
 
   real :: regrid_time_scale !< The time-scale used in blending between the current (old) grid
                             !! and the target (new) grid [T ~> s]
 
   type(regridding_CS) :: regridCS !< Regridding parameters and work arrays
   type(remapping_CS)  :: remapCS  !< Remapping parameters and work arrays
 
-  integer :: nk              !< Used only for queries, not directly by this module
+  integer :: nk             !< Used only for queries, not directly by this module
 
-  logical :: remap_after_initialization !<   Indicates whether to regrid/remap after initializing the state.
+  logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state.
+
+  logical :: answers_2018   !< If true, use the order of arithmetic and expressions for remapping
+                            !! that recover the answers from the end of 2018.  Otherwise, use more
+                            !! robust and accurate forms of mathematically equivalent expressions.
 
   logical :: show_call_tree !< For debugging
 
@@ -145,6 +149,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   character(len=40)               :: mdl = "MOM_ALE" ! This module's name.
   character(len=80)               :: string ! Temporary strings
   real                            :: filter_shallow_depth, filter_deep_depth
+  logical       :: default_2018_answers
   logical                         :: check_reconstruction
   logical                         :: check_remapping
   logical                         :: force_bounds_in_subcell
@@ -192,11 +197,19 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
                  "If true, values at the interfaces of boundary cells are "//&
                  "extrapolated instead of piecewise constant", default=.false.)
+  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
+                 "This sets the default value for the various _2018_ANSWERS parameters.", &
+                 default=.true.)
+  call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", CS%answers_2018, &
+                 "If true, use the order of arithmetic and expressions that recover the "//&
+                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
+                 "forms of the same expressions.", default=default_2018_answers)
   call initialize_remapping( CS%remapCS, string, &
                              boundary_extrapolation=remap_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
-                             force_bounds_in_subcell=force_bounds_in_subcell)
+                             force_bounds_in_subcell=force_bounds_in_subcell, &
+                             answers_2018=CS%answers_2018)
 
   call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", CS%remap_after_initialization, &
                  "If true, applies regridding and remapping immediately after "//&
@@ -220,7 +233,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
                  "REGRID_FILTER_SHALLOW_DEPTH the filter weights adopt a cubic profile.", &
                  units="m", default=0., scale=GV%m_to_H)
   call set_regrid_params(CS%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &
-                                      depth_of_time_filter_deep=filter_deep_depth)
+                         depth_of_time_filter_deep=filter_deep_depth)
   call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, &
                  "If true, the regridding ntegrates upwards from the bottom for "//&
                  "interface positions, much as the main model does. If false "//&
@@ -1089,13 +1102,13 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
 
   ! Local variables
   integer :: i, j, k
-  real    :: hTmp(GV%ke)
-  real    :: tmp(GV%ke)
+  real    :: hTmp(GV%ke) ! A 1-d copy of h [H ~> m or kg m-2]
+  real    :: tmp(GV%ke)  ! A 1-d copy of a column of temperature [degC] or salinity [ppt]
   real, dimension(CS%nk,2) :: &
-      ppol_E            !Edge value of polynomial
+      ppol_E            ! Edge value of polynomial in [degC] or [ppt]
   real, dimension(CS%nk,3) :: &
-      ppol_coefs !Coefficients of polynomial
-  real :: h_neglect, h_neglect_edge
+      ppol_coefs        ! Coefficients of polynomial, all in [degC] or [ppt]
+  real :: h_neglect, h_neglect_edge ! Tiny thicknesses [H ~> m or kg m-2]
 
   !### Try replacing both of these with GV%H_subroundoff
   if (GV%Boussinesq) then
@@ -1116,7 +1129,8 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
     ppol_E(:,:) = 0.0
     ppol_coefs(:,:) = 0.0
     !### Try to replace the following value of h_neglect with GV%H_subroundoff.
-    call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=h_neglect_edge )
+    call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=h_neglect_edge, &
+                                  answers_2018=CS%answers_2018 )
     call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
     if (bdry_extrap) &
       call PPM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
@@ -1131,7 +1145,8 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
     ppol_coefs(:,:) = 0.0
     tmp(:) = tv%T(i,j,:)
     !### Try to replace the following value of h_neglect with GV%H_subroundoff.
-    call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=1.0e-10*GV%m_to_H )
+    call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=1.0e-10*GV%m_to_H, &
+                                  answers_2018=CS%answers_2018 )
     call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
     if (bdry_extrap) &
       call PPM_boundary_extrapolation(GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )

src/ALE/MOM_remapping.F90

@@ -33,6 +33,8 @@ module MOM_remapping
   logical :: check_remapping = .false.
   !> If true, the intermediate values used in remapping are forced to be bounded.
   logical :: force_bounds_in_subcell = .false.
+  !> If true use older, less acccurate expressions.
+  logical :: answers_2018 = .true.
 end type
 
 ! The following routines are visible to the outside world
@@ -84,13 +86,14 @@ module MOM_remapping
 
 !> Set parameters within remapping object
 subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation,  &
-               check_reconstruction, check_remapping, force_bounds_in_subcell)
+               check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
   type(remapping_CS),         intent(inout) :: CS !< Remapping control structure
   character(len=*), optional, intent(in)    :: remapping_scheme !< Remapping scheme to use
   logical, optional,          intent(in)    :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
   logical, optional,          intent(in)    :: check_reconstruction !< Indicate to check reconstructions
   logical, optional,          intent(in)    :: check_remapping !< Indicate to check results of remapping
   logical, optional,          intent(in)    :: force_bounds_in_subcell !< Force subcells values to be bounded
+  logical, optional,          intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
 
   if (present(remapping_scheme)) then
     call setReconstructionType( remapping_scheme, CS )
@@ -107,6 +110,9 @@ subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation,  &
   if (present(force_bounds_in_subcell)) then
     CS%force_bounds_in_subcell = force_bounds_in_subcell
   endif
+  if (present(answers_2018)) then
+    CS%answers_2018 = answers_2018
+  endif
 end subroutine remapping_set_param
 
 subroutine extract_member_remapping_CS(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, &
@@ -392,31 +398,31 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, &
       endif
       iMethod = INTEGRATION_PLM
     case ( REMAPPING_PPM_H4 )
-      call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge )
+      call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
       call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect )
       if ( CS%boundary_extrapolation ) then
         call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect )
       endif
       iMethod = INTEGRATION_PPM
     case ( REMAPPING_PPM_IH4 )
-      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge )
+      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
       call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect )
       if ( CS%boundary_extrapolation ) then
         call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect )
       endif
       iMethod = INTEGRATION_PPM
     case ( REMAPPING_PQM_IH4IH3 )
-      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge )
-      call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect )
+      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
+      call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect, answers_2018=CS%answers_2018 )
       call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect )
       if ( CS%boundary_extrapolation ) then
         call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, &
                                             ppoly_r_coefs, h_neglect )
       endif
       iMethod = INTEGRATION_PQM
     case ( REMAPPING_PQM_IH6IH5 )
-      call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge )
-      call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S, h_neglect )
+      call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
+      call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S, h_neglect, answers_2018=CS%answers_2018 )
       call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect )
       if ( CS%boundary_extrapolation ) then
         call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, &
@@ -1537,19 +1543,20 @@ end subroutine dzFromH1H2
 
 !> Constructor for remapping control structure
 subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, &
-                check_reconstruction, check_remapping, force_bounds_in_subcell)
+                check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
   ! Arguments
   type(remapping_CS), intent(inout) :: CS !< Remapping control structure
   character(len=*),   intent(in)    :: remapping_scheme !< Remapping scheme to use
   logical, optional,  intent(in)    :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
   logical, optional,  intent(in)    :: check_reconstruction !< Indicate to check reconstructions
   logical, optional,  intent(in)    :: check_remapping !< Indicate to check results of remapping
   logical, optional,  intent(in)    :: force_bounds_in_subcell !< Force subcells values to be bounded
+  logical, optional,  intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
 
-  ! Note that remapping_scheme is mandatory fir initialize_remapping()
+  ! Note that remapping_scheme is mandatory for initialize_remapping()
   call remapping_set_param(CS, remapping_scheme=remapping_scheme, boundary_extrapolation=boundary_extrapolation,  &
                check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
-               force_bounds_in_subcell=force_bounds_in_subcell)
+               force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018)
 
 end subroutine initialize_remapping
 
@@ -1615,13 +1622,15 @@ logical function remapping_unit_tests(verbose)
   data h2 /6*0.5/  ! 6 uniform layers with total depth of 3
   type(remapping_CS) :: CS !< Remapping control structure
   real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S, ppoly0_coefs
+  logical :: answers_2018 !  If true use older, less acccurate expressions.
   integer :: i
   real :: err, h_neglect, h_neglect_edge
   logical :: thisTest, v
 
   v = verbose
   h_neglect = hNeglect_dflt
   h_neglect_edge = 1.0e-10
+  answers_2018 = .true.
 
   write(*,*) '==== MOM_remapping: remapping_unit_tests ================='
   remapping_unit_tests = .false. ! Normally return false
@@ -1643,7 +1652,7 @@ logical function remapping_unit_tests(verbose)
   remapping_unit_tests = remapping_unit_tests .or. thisTest
 
   thisTest = .false.
-  call initialize_remapping(CS, 'PPM_H4')
+  call initialize_remapping(CS, 'PPM_H4', answers_2018=answers_2018)
   if (verbose) write(*,*) 'h0 (test data)'
   if (verbose) call dumpGrid(n0,h0,x0,u0)
 
@@ -1667,7 +1676,7 @@ logical function remapping_unit_tests(verbose)
   ppoly0_S(:,:) = 0.0
   ppoly0_coefs(:,:) = 0.0
 
-  call edge_values_explicit_h4( n0, h0, u0, ppoly0_E, h_neglect=1e-10 )
+  call edge_values_explicit_h4( n0, h0, u0, ppoly0_E, h_neglect=1e-10, answers_2018=answers_2018 )
   call PPM_reconstruction( n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect )
   call PPM_boundary_extrapolation( n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect )
   u1(:) = 0.
@@ -1798,7 +1807,7 @@ logical function remapping_unit_tests(verbose)
     test_answer(v, 3, ppoly0_coefs(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1')
 
   call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_E, &
-                                h_neglect=1e-10 )
+                                h_neglect=1e-10, answers_2018=answers_2018 )
   ! The next two tests currently fail due to roundoff.
   thisTest = test_answer(v, 5, ppoly0_E(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges')
   thisTest = test_answer(v, 5, ppoly0_E(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges')
@@ -1814,7 +1823,7 @@ logical function remapping_unit_tests(verbose)
     test_answer(v, 5, ppoly0_coefs(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2')
 
   call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_E, &
-                                h_neglect=1e-10 )
+                                h_neglect=1e-10, answers_2018=answers_2018 )
   ! The next two tests currently fail due to roundoff.
   thisTest = test_answer(v, 5, ppoly0_E(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges')
   thisTest = test_answer(v, 5, ppoly0_E(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges')