Fix remapping of last layer for mismatched ocean depths

- The original remap_via_subcells() assumed that the first and last
  sub-layers would be vanished and thus needed to only assign the
  edge values of the source reconstructions to the sub-layers.
  However, this is only a valid assumption (and correct) if the
  total column thickness of source/target are the same. That is
  generally true doing the main loop. However, when initializing
  from WOA, the ocean model depth and source-data depth often differ,
  and this assumption that we can ignore the top/bottom reconstructions
  breaks.
- The original fix was developed with @claireyung and took the form
  claireyung@924b7ac
  In this patch, I have unrolled the loop inside remap_src_to_sub_grid() to
  avoid adding an additional `if` test on the loop index.
- The fix is implemented in remap_src_to_sub_grid() and the original
  method is reached by setting a logical inside the remapping_CS. Default
  is to use the OM4 alorithm (original method with bug).
- New runtime parameters are added in to recover original algorithm selectively:
  - OBC_REMAPPING_USE_OM4_SUBCELLS,
  - Z_INIT_REMAPPING_USE_OM4_SUBCELLS,
  - EBT_REMAPPING_USE_OM4_SUBCELLS,
  - SPONGE_REMAPPING_USE_OM4_SUBCELLS,
  - SPONGE_REMAPPING_USE_OM4_SUBCELLS,
  - DIAG_REMAPPING_USE_OM4_SUBCELLS,
  - NDIFF_REMAPPING_USE_OM4_SUBCELLS,
  - HBD_REMAPPING_USE_OM4_SUBCELLS,
  - INTWAVE_REMAPPING_USE_OM4_SUBCELLS, and
  - REMAPPING_USE_OM4_SUBCELLS
  all of which default to True.
- No answer changes.

src/ALE/MOM_ALE.F90

@@ -178,6 +178,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   logical           :: local_logical
   logical           :: remap_boundary_extrap
   logical           :: init_boundary_extrap
+  logical           :: om4_remap_via_sub_cells
   type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding
                                                          ! for sharing parameters.
 
@@ -234,6 +235,11 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
                  default=99991231)
+  call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 "This selects the remapping algorithm used in OM4 that does not use "//&
+                 "the full reconstruction for the top- and lower-most sub-layers, but instead "//&
+                 "assumes they are always vanished (untrue) and so just uses their edge values. "//&
+                 "We recommend setting this option to false.", default=.true.)
   call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", CS%answer_date, &
                  "The vintage of the expressions and order of arithmetic to use for remapping.  "//&
                  "Values below 20190101 result in the use of older, less accurate expressions "//&
@@ -247,12 +253,14 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
+                             om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
                              answer_date=CS%answer_date)
   call initialize_remapping( CS%vel_remapCS, vel_string, &
                              boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
+                             om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
                              answer_date=CS%answer_date)
 
   call get_param(param_file, mdl, "PARTIAL_CELL_VELOCITY_REMAP", CS%partial_cell_vel_remap, &
@@ -326,6 +334,21 @@ subroutine ALE_set_extrap_boundaries( param_file, CS)
   call remapping_set_param(CS%remapCS, boundary_extrapolation=remap_boundary_extrap)
 end subroutine ALE_set_extrap_boundaries
 
+!> Sets the remapping algorithm to that of OM4
+!!
+!! The remapping aglorithm used in OM4 made poor assumptions about the reconstructions
+!! in the top/bottom layers, namely that they were always vanished and could be
+!! represented solely by their upper/lower edge value respectively.
+!! Passing .false. here uses the full reconstruction of those top and bottom layers
+!! and properly sample those layers.
+subroutine ALE_set_OM4_remap_algorithm( CS, om4_remap_via_sub_cells )
+  type(ALE_CS), pointer :: CS !< Module control structure
+  logical, intent(in)   :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm
+
+  call remapping_set_param(CS%remapCS, om4_remap_via_sub_cells =om4_remap_via_sub_cells )
+
+end subroutine ALE_set_OM4_remap_algorithm
+
 !> Initialize diagnostics for the ALE module.
 subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
   type(time_type),target,     intent(in)  :: Time  !< Time structure

src/core/MOM_open_boundary.F90

@@ -441,6 +441,7 @@ subroutine open_boundary_config(G, US, param_file, OBC)
   real               :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m]
   integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags.
   logical :: check_reconstruction, check_remapping, force_bounds_in_subcell
+  logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
   character(len=64)  :: remappingScheme
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
@@ -695,10 +696,15 @@ subroutine open_boundary_config(G, US, param_file, OBC)
                  "that were in use at the end of 2018.  Higher values result in the use of more "//&
                  "robust and accurate forms of mathematically equivalent expressions.", &
                  default=default_answer_date)
+    call get_param(param_file, mdl, "OBC_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//&
+                 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
+                 "We recommend setting this option to false.", default=.true.)
 
     allocate(OBC%remap_CS)
     call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., &
                check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
+               om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
                force_bounds_in_subcell=force_bounds_in_subcell, answer_date=OBC%remap_answer_date)
 
   endif ! OBC%number_of_segments > 0

src/diagnostics/MOM_diagnostics.F90

@@ -1580,6 +1580,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   logical :: better_speed_est ! If true, use a more robust estimate of the first
                               ! mode wave speed as the starting point for iterations.
   logical :: split            ! True if using the barotropic-baroclinic split algorithm
+  logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "MOM_diagnostics" ! This module's name.
@@ -1617,6 +1618,10 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
                  "If true, use a more robust estimate of the first mode wave speed as the "//&
                  "starting point for iterations.", default=.true.)
+  call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//&
+                 "See REMAPPING_USE_OM4_SUBCELLS for details. "//&
+                 "We recommend setting this option to false.", default=.true.)
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
                  default=99991231)
@@ -1858,7 +1863,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       (CS%id_cg_ebt>0) .or. (CS%id_Rd_ebt>0) .or. (CS%id_p_ebt>0)) then
     call wave_speed_init(CS%wave_speed, remap_answer_date=remap_answer_date, &
                          better_speed_est=better_speed_est, min_speed=wave_speed_min, &
-                         wave_speed_tol=wave_speed_tol)
+                         wave_speed_tol=wave_speed_tol, om4_remap_via_sub_cells=om4_remap_via_sub_cells)
   endif
 
   CS%id_mass_wt = register_diag_field('ocean_model', 'mass_wt', diag%axesT1, Time, &

src/diagnostics/MOM_wave_speed.F90

@@ -1611,7 +1611,8 @@ end subroutine tridiag_det
 
 !> Initialize control structure for MOM_wave_speed
 subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, &
-                           remap_answer_date, better_speed_est, min_speed, wave_speed_tol, c1_thresh)
+                           remap_answer_date, better_speed_est, om4_remap_via_sub_cells, &
+                           min_speed, wave_speed_tol, c1_thresh)
   type(wave_speed_CS), intent(inout) :: CS  !< Wave speed control struct
   logical, optional, intent(in) :: use_ebt_mode  !< If true, use the equivalent
                                      !! barotropic mode instead of the first baroclinic mode.
@@ -1630,6 +1631,8 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de
                                       !! forms of the same remapping expressions.
   logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
                                      !! mode speed as the starting point for iterations.
+  logical, optional, intent(in) :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells
+                                     !! for calculating the EBT structure
   real,    optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
                                      !! below which 0 is returned [L T-1 ~> m s-1].
   real,    optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
@@ -1656,6 +1659,7 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de
 
   ! The remap_answers_2018 argument here is irrelevant, because remapping is hard-coded to use PLM.
   call initialize_remapping(CS%remapping_CS, 'PLM', boundary_extrapolation=.false., &
+                            om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
                             answer_date=CS%remap_answer_date)
 
 end subroutine wave_speed_init

src/framework/MOM_diag_mediator.F90

@@ -3251,6 +3251,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
                                   ! answers from 2018, while higher values use more robust
                                   ! forms of the same remapping expressions.
   integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags.
+  logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for diagnostics
   character(len=8)   :: this_pe
   character(len=240) :: doc_file, doc_file_dflt, doc_path
   character(len=240), allocatable :: diag_coords(:)
@@ -3282,6 +3283,10 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
                  default=99991231)
+  call get_param(param_file, mdl, "DIAG_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 "If true, use the OM4 remapping-via-subcells algorithm for diagnostics. "//&
+                 "See REMAPPING_USE_OM4_SUBCELLS for details. "//&
+                 "We recommend setting this option to false.", default=.true.)
   call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
                  "The vintage of the expressions and order of arithmetic to use for remapping.  "//&
                  "Values below 20190101 result in the use of older, less accurate expressions "//&
@@ -3311,7 +3316,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
     allocate(diag_cs%diag_remap_cs(diag_cs%num_diag_coords))
     ! Initialize each diagnostic vertical coordinate
     do i=1, diag_cs%num_diag_coords
-      call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), answer_date=remap_answer_date, GV=GV)
+      call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), om4_remap_via_sub_cells, remap_answer_date, GV)
     enddo
     deallocate(diag_coords)
   endif

src/framework/MOM_diag_remap.F90

@@ -92,6 +92,7 @@ module MOM_diag_remap
                                       !! vertical extents in [Z ~> m] for remapping extensive variables
   integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces
   integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers
+  logical :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells
   integer :: answer_date      !< The vintage of the order of arithmetic and expressions
                               !! to use for remapping.  Values below 20190101 recover
                               !! the answers from 2018, while higher values use more
@@ -102,10 +103,11 @@ module MOM_diag_remap
 contains
 
 !> Initialize a diagnostic remapping type with the given vertical coordinate.
-subroutine diag_remap_init(remap_cs, coord_tuple, answer_date, GV)
+subroutine diag_remap_init(remap_cs, coord_tuple, om4_remap_via_sub_cells, answer_date, GV)
   type(diag_remap_ctrl), intent(inout) :: remap_cs    !< Diag remapping control structure
   character(len=*),      intent(in)    :: coord_tuple !< A string in form of
                                                       !! MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME
+  logical,               intent(in)    :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells
   integer,               intent(in)    :: answer_date !< The vintage of the order of arithmetic and expressions
                                                       !! to use for remapping.  Values below 20190101 recover
                                                       !! the answers from 2018, while higher values use more
@@ -127,6 +129,7 @@ subroutine diag_remap_init(remap_cs, coord_tuple, answer_date, GV)
   remap_cs%configured = .false.
   remap_cs%initialized = .false.
   remap_cs%used = .false.
+  remap_cs%om4_remap_via_sub_cells = om4_remap_via_sub_cells
   remap_cs%answer_date = answer_date
   remap_cs%nz = 0
 
@@ -309,6 +312,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe
   if (.not. remap_cs%initialized) then
     ! Initialize remapping and regridding on the first call
     call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false., &
+                              om4_remap_via_sub_cells=remap_cs%om4_remap_via_sub_cells, &
                               answer_date=remap_cs%answer_date)
     remap_cs%initialized = .true.
   endif

src/initialization/MOM_state_initialization.F90

@@ -2504,6 +2504,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
   character(len=64) :: remappingScheme
   real :: tempAvg  ! Spatially averaged temperatures on a layer [C ~> degC]
   real :: saltAvg  ! Spatially averaged salinities on a layer [S ~> ppt]
+  logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm (only used if useALEremapping)
   logical :: do_conv_adj, ignore
   integer :: nPoints
   integer :: id_clock_routine, id_clock_ALE
@@ -2583,6 +2584,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
                  "that were in use at the end of 2018.  Higher values result in the use of more "//&
                  "robust and accurate forms of mathematically equivalent expressions.", &
                  default=default_answer_date, do_not_log=just_read.or.(.not.GV%Boussinesq))
+    call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//&
+                 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
+                 "We recommend setting this option to false.", default=.true.)
     if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
   endif
   call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, &
@@ -2753,7 +2758,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
     ! Build the target grid (and set the model thickness to it)
 
     call ALE_initRegridding( GV, US, G%max_depth, PF, mdl, regridCS ) ! sets regridCS
-    call initialize_remapping( remapCS, remappingScheme, boundary_extrapolation=.false., answer_date=remap_answer_date )
+    call initialize_remapping( remapCS, remappingScheme, boundary_extrapolation=.false., &
+                               om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date )
 
     ! Now remap from source grid to target grid, first setting reconstruction parameters
     if (remap_general) then

src/initialization/MOM_tracer_initialization_from_Z.F90

@@ -96,6 +96,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_
                                   ! remapping cell reconstructions [Z ~> m]
   real :: dz_neglect_edge         ! A negligibly small vertical layer extent used in
                                   ! remapping edge value calculations [Z ~> m]
+  logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
   integer :: nPoints    ! The number of valid input data points in a column
   integer :: id_clock_routine, id_clock_ALE
   integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags.
@@ -137,6 +138,10 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_
                  "that were in use at the end of 2018.  Higher values result in the use of more "//&
                  "robust and accurate forms of mathematically equivalent expressions.", &
                  default=default_answer_date, do_not_log=.not.GV%Boussinesq)
+    call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//&
+                 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
+                 "We recommend setting this option to false.", default=.true.)
     if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
   endif
   call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, &
@@ -174,7 +179,8 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_
     allocate( dzSrc(isd:ied,jsd:jed,kd) )
     allocate( hSrc(isd:ied,jsd:jed,kd) )
     ! Set parameters for reconstructions
-    call initialize_remapping( remapCS, remapScheme, boundary_extrapolation=.false., answer_date=remap_answer_date )
+    call initialize_remapping( remapCS, remapScheme, boundary_extrapolation=.false., &
+                               om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date )
     ! Next we initialize the regridding package so that it knows about the target grid
 
     do j = js, je ; do i = is, ie

src/ocean_data_assim/MOM_oda_driver.F90

@@ -183,6 +183,7 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS)
   character(len=80) :: remap_scheme
   character(len=80) :: bias_correction_file, inc_file
   integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags.
+  logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
 
   if (associated(CS)) call MOM_error(FATAL, 'Calling oda_init with associated control structure')
   allocate(CS)
@@ -320,8 +321,10 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS)
   call get_param(PF, 'oda_driver', "REGRIDDING_COORDINATE_MODE", coord_mode, &
        "Coordinate mode for vertical regridding.", &
        default="ZSTAR", fail_if_missing=.false.)
+  call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 do_not_log=.true., default=.true.)
   call initialize_regridding(CS%regridCS, CS%GV, CS%US, dG%max_depth,PF,'oda_driver',coord_mode,'','')
-  call initialize_remapping(CS%remapCS,remap_scheme)
+  call initialize_remapping(CS%remapCS, remap_scheme, om4_remap_via_sub_cells=om4_remap_via_sub_cells)
   call set_regrid_params(CS%regridCS, min_thickness=0.)
   isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed