+Add REMAP_AUXILIARY_VARS to remap accelerations

  Added REMAP_AUXILIARY_VARS to control whether to remap the accelerations that
are used in the predictor step of the split RK2 time stepping scheme.  Also
added the new routines remap_dyn_split_rk2_aux_vars, remap_OBC_fields and
remap_vertvisc_aux_vars to do the remapping, and code to call these routines
when REMAP_AUXILIARY_VARS is true. By default, REMAP_AUXILIARY_VARS is false,
and all answers are bitwise identical, but the entire MOM6-examples regression
suite has been run with this set to true, and they do appear to give physically
plausible answers in all cases, partially addressing the issue noted at
github.com//issues/203.  New entries are added to the
MOM_parameter_doc files, and there are three new publicly visible routines, but
by default answers do not change.

src/ALE/MOM_ALE.F90

@@ -291,10 +291,10 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   call set_regrid_params(CS%regridCS, depth_of_time_filter_shallow=filter_shallow_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 "//&
+                 "If true, the regridding integrates upwards from the bottom for "//&
                  "interface positions, much as the main model does. If false "//&
-                 "regridding integrates downward, consistant with the remapping "//&
-                 "code.", default=.true., do_not_log=.true.)
+                 "regridding integrates downward, consistent with the remapping code.", &
+                 default=.true., do_not_log=.true.)
   call set_regrid_params(CS%regridCS, integrate_downward_for_e=.not.local_logical)
 
   call get_param(param_file, mdl, "REMAP_VEL_MASK_BBL_THICK", CS%BBL_h_vel_mask, &

src/core/MOM.F90

@@ -74,7 +74,7 @@ module MOM
 use MOM_dynamics_unsplit,      only : MOM_dyn_unsplit_CS
 use MOM_dynamics_split_RK2,    only : step_MOM_dyn_split_RK2, register_restarts_dyn_split_RK2
 use MOM_dynamics_split_RK2,    only : initialize_dyn_split_RK2, end_dyn_split_RK2
-use MOM_dynamics_split_RK2,    only : MOM_dyn_split_RK2_CS
+use MOM_dynamics_split_RK2,    only : MOM_dyn_split_RK2_CS, remap_dyn_split_rk2_aux_vars
 use MOM_dynamics_unsplit_RK2,  only : step_MOM_dyn_unsplit_RK2, register_restarts_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : MOM_dyn_unsplit_RK2_CS
@@ -103,14 +103,13 @@ module MOM
 use MOM_mixed_layer_restrat,   only : mixedlayer_restrat_register_restarts
 use MOM_obsolete_diagnostics,  only : register_obsolete_diagnostics
 use MOM_open_boundary,         only : ocean_OBC_type, OBC_registry_type
-use MOM_open_boundary,         only : register_temp_salt_segments
-use MOM_open_boundary,         only : open_boundary_register_restarts
-use MOM_open_boundary,         only : update_segment_tracer_reservoirs
+use MOM_open_boundary,         only : register_temp_salt_segments, update_segment_tracer_reservoirs
+use MOM_open_boundary,         only : open_boundary_register_restarts, remap_OBC_fields
 use MOM_open_boundary,         only : rotate_OBC_config, rotate_OBC_init
 use MOM_porous_barriers,       only : porous_widths_layer, porous_widths_interface, porous_barriers_init
 use MOM_porous_barriers,       only : porous_barrier_CS
-use MOM_set_visc,              only : set_viscous_BBL, set_viscous_ML
-use MOM_set_visc,              only : set_visc_register_restarts, set_visc_CS
+use MOM_set_visc,              only : set_viscous_BBL, set_viscous_ML, set_visc_CS
+use MOM_set_visc,              only : set_visc_register_restarts, remap_vertvisc_aux_vars
 use MOM_set_visc,              only : set_visc_init, set_visc_end
 use MOM_shared_initialization, only : write_ocean_geometry_file
 use MOM_sponge,                only : init_sponge_diags, sponge_CS
@@ -250,6 +249,10 @@ module MOM
   logical :: use_ALE_algorithm  !< If true, use the ALE algorithm rather than layered
                     !! isopycnal/stacked shallow water mode. This logical is set by calling the
                     !! function useRegridding() from the MOM_regridding module.
+  logical :: remap_aux_vars     !< If true, apply ALE remapping to all of the auxiliary 3-D
+                    !! variables that are needed to reproduce across restarts,
+                    !! similarly to what is done with the primary state variables.
+
   type(MOM_stoch_eos_CS) :: stoch_eos_CS !< structure containing random pattern for stoch EOS
   logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction
                     !! updates occur first in directionally split parts of the calculation.
@@ -1547,6 +1550,16 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
       call ALE_remap_tracers(CS%ALE_CSp, G, GV, h, h_new, CS%tracer_Reg, showCallTree, dtdia, PCM_cell)
       call ALE_remap_velocities(CS%ALE_CSp, G, GV, h, h_new, u, v, CS%OBC, dzRegrid, showCallTree, dtdia)
 
+      if (CS%remap_aux_vars) then
+        if (CS%split) &
+          call remap_dyn_split_RK2_aux_vars(G, GV, CS%dyn_split_RK2_CSp, h, h_new, CS%ALE_CSp, CS%OBC, dzRegrid)
+
+        if (associated(CS%OBC)) &
+          call remap_OBC_fields(G, GV, h, h_new, CS%OBC, PCM_cell=PCM_cell)
+
+        call remap_vertvisc_aux_vars(G, GV, CS%visc, h, h_new, CS%ALE_CSp, CS%OBC)
+      endif
+
       ! Replace the old grid with new one.  All remapping must be done by this point in the code.
       !$OMP parallel do default(shared)
       do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
@@ -2072,6 +2085,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   call get_param(param_file, "MOM", "USE_REGRIDDING", CS%use_ALE_algorithm, &
                  "If True, use the ALE algorithm (regridding/remapping). "//&
                  "If False, use the layered isopycnal algorithm.", default=.false. )
+  call get_param(param_file, "MOM", "REMAP_AUXILIARY_VARS", CS%remap_aux_vars, &
+                 "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&
+                 "variables that are needed to reproduce across restarts, similarly to "//&
+                 "what is already being done with the primary state variables.  "//&
+                 "The default should be changed to true.", default=.false., &
+                 do_not_log=.not.CS%use_ALE_algorithm)
   call get_param(param_file, "MOM", "BULKMIXEDLAYER", bulkmixedlayer, &
                  "If true, use a Kraus-Turner-like bulk mixed layer "//&
                  "with transitional buffer layers.  Layers 1 through "//&

src/core/MOM_dynamics_split_RK2.F90

@@ -36,7 +36,7 @@ module MOM_dynamics_split_RK2
 use MOM_time_manager,      only : time_type, time_type_to_real, operator(+)
 use MOM_time_manager,      only : operator(-), operator(>), operator(*), operator(/)
 
-use MOM_ALE,                   only : ALE_CS
+use MOM_ALE,                   only : ALE_CS, ALE_remap_velocities
 use MOM_barotropic,            only : barotropic_init, btstep, btcalc, bt_mass_source
 use MOM_barotropic,            only : register_barotropic_restarts, set_dtbt, barotropic_CS
 use MOM_barotropic,            only : barotropic_end
@@ -160,6 +160,9 @@ module MOM_dynamics_split_RK2
                                   !! predictor step.  This is used to accomodate various generations
                                   !! of restart files.
   logical :: use_tides            !< If true, tidal forcing is enabled.
+  logical :: remap_aux            !< If true, apply ALE remapping to all of the auxiliary 3-D
+                                  !! variables that are needed to reproduce across restarts,
+                                  !! similarly to what is done with the primary state variables.
 
   real    :: be      !< A nondimensional number from 0.5 to 1 that controls
                      !! the backward weighting of the time stepping scheme [nondim]
@@ -256,6 +259,7 @@ module MOM_dynamics_split_RK2
 public step_MOM_dyn_split_RK2
 public register_restarts_dyn_split_RK2
 public initialize_dyn_split_RK2
+public remap_dyn_split_RK2_aux_vars
 public end_dyn_split_RK2
 
 !>@{ CPU time clock IDs
@@ -1160,6 +1164,32 @@ subroutine register_restarts_dyn_split_RK2(HI, GV, US, param_file, CS, restart_C
 
 end subroutine register_restarts_dyn_split_RK2
 
+!> This subroutine does remapping for the auxiliary restart variables that are used
+!! with the split RK2 time stepping scheme.
+subroutine remap_dyn_split_RK2_aux_vars(G, GV, CS, h_old, h_new, ALE_CSp, OBC, dzRegrid)
+  type(ocean_grid_type),            intent(inout) :: G        !< ocean grid structure
+  type(verticalGrid_type),          intent(in)    :: GV       !< ocean vertical grid structure
+  type(MOM_dyn_split_RK2_CS),       pointer       :: CS       !< module control structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                                    intent(in)    :: h_old    !< Thickness of source grid  [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                                    intent(in)    :: h_new    !< Thickness of destination grid [H ~> m or kg m-2]
+  type(ALE_CS),                     pointer       :: ALE_CSp  !< ALE control structure to use when remapping
+  type(ocean_OBC_type),             pointer       :: OBC      !< OBC control structure to use when remapping
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                          optional, intent(in)    :: dzRegrid !< Change in interface position [H ~> m or kg m-2]
+
+  if (.not.CS%remap_aux) return
+
+  if (CS%store_CAu) then
+    call ALE_remap_velocities(ALE_CSp, G, GV, h_old, h_new, CS%u_av, CS%v_av, OBC, dzRegrid)
+    call ALE_remap_velocities(ALE_CSp, G, GV, h_old, h_new, CS%CAu_pred, CS%CAv_pred, OBC, dzRegrid)
+  endif
+
+  call ALE_remap_velocities(ALE_CSp, G, GV, h_old, h_new, CS%diffu, CS%diffv, OBC, dzRegrid)
+
+end subroutine remap_dyn_split_RK2_aux_vars
+
 !> This subroutine initializes all of the variables that are used by this
 !! dynamic core, including diagnostics and the cpu clocks.
 subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, &
@@ -1276,6 +1306,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
                  "If true, calculate the Coriolis accelerations at the end of each "//&
                  "timestep for use in the predictor step of the next split RK2 timestep.", &
                  default=.true.)
+  call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", CS%remap_aux, &
+                 "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&
+                 "variables that are needed to reproduce across restarts, similarly to "//&
+                 "what is already being done with the primary state variables.  "//&
+                 "The default should be changed to true.", default=.false., do_not_log=.true.)
+  if (CS%remap_aux .and. .not.CS%store_CAu) call MOM_error(FATAL, &
+      "REMAP_AUXILIARY_VARS requires that STORE_CORIOLIS_ACCEL = True.")
   call get_param(param_file, mdl, "DEBUG", CS%debug, &
                  "If true, write out verbose debugging data.", &
                  default=.false., debuggingParam=.true.)

src/core/MOM_open_boundary.F90

@@ -65,6 +65,7 @@ module MOM_open_boundary
 public open_boundary_register_restarts
 public update_segment_tracer_reservoirs
 public update_OBC_ramp
+public remap_OBC_fields
 public rotate_OBC_config
 public rotate_OBC_init
 public initialize_segment_data
@@ -5408,6 +5409,184 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
 
 end subroutine update_segment_tracer_reservoirs
 
+!> Vertically remap the OBC tracer reservoirs and radiation rates that are filtered in time.
+subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell)
+  type(ocean_grid_type),                     intent(in) :: G     !< The ocean's grid structure
+  type(verticalGrid_type),                   intent(in) :: GV    !<  Ocean vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid [H ~> m or kg m-2]
+  type(ocean_OBC_type),                      pointer    :: OBC   !< Open boundary structure
+  logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                                   optional, intent(in) :: PCM_cell !< Use PCM remapping in cells where true
+
+  ! Local variables
+  type(OBC_segment_type), pointer :: segment => NULL() ! A pointer to the various segments, used just for shorthand.
+
+  real :: tr_column(GV%ke)  ! A column of updated tracer concentrations [CU ~> Conc]
+  real :: r_norm_col(GV%ke) ! A column of updated radiation rates, in grid points per timestep [nondim]
+  real :: rxy_col(GV%ke) ! A column of updated radiation rates for oblique OBCs [L2 T-2 ~> m2 s-2]
+  real :: h1(GV%ke)     ! A column of source grid layer thicknesses [H ~> m or kg m-2]
+  real :: h2(GV%ke)     ! A column of target grid layer thicknesses [H ~> m or kg m-2]
+  real :: I_scale       ! The inverse of the scaling factor for the tracers.
+                        ! For salinity the units would be [ppt S-1 ~> 1].
+  real :: h_neglect     ! Tiny thickness used in remapping [H ~> m or kg m-2]
+  logical :: PCM(GV%ke) ! If true, do PCM remapping from a cell.
+  integer :: i, j, k, m, n, ntr, nz
+
+  if (.not.associated(OBC)) return
+
+  nz = GV%ke
+  ntr = OBC%ntr
+  h_neglect = GV%H_subroundoff
+
+  if (.not.present(PCM_cell)) PCM(:) = .false.
+
+  if (associated(OBC)) then ; if (OBC%OBC_pe) then ; do n=1,OBC%number_of_segments
+    segment => OBC%segment(n)
+    if (.not.associated(segment%tr_Reg)) cycle
+
+    if (segment%is_E_or_W) then
+      I = segment%HI%IsdB
+      do j=segment%HI%jsd,segment%HI%jed
+
+        ! Store a column of the start and final grids
+        if (segment%direction == OBC_DIRECTION_W) then
+          if (G%mask2dT(i+1,j) == 0.0) cycle
+          h1(:) = h_old(i+1,j,:)
+          h2(:) = h_new(i+1,j,:)
+          if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i+1,j,:) ; endif
+        else
+          if (G%mask2dT(i,j) == 0.0) cycle
+          h1(:) = h_old(i,j,:)
+          h2(:) = h_new(i,j,:)
+          if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i,j,:) ; endif
+        endif
+
+        ! Vertically remap the reservoir tracer concentrations
+        do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then
+          I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale
+
+          if (present(PCM_cell)) then
+            call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column, &
+                                  h_neglect, h_neglect, PCM_cell=PCM)
+          else
+            call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column, &
+                                  h_neglect, h_neglect)
+          endif
+
+          ! Possibly underflow any very tiny tracer concentrations to 0?
+
+          ! Update tracer concentrations
+          segment%tr_Reg%Tr(m)%tres(I,j,:) = tr_column(:)
+          if (allocated(OBC%tres_x)) then ; do k=1,nz
+            OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k)
+          enddo ; endif
+
+        endif ; enddo
+
+        if (segment%radiation .and. (OBC%gamma_uv < 1.0)) then
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_rad(I,j,:), nz, h2, r_norm_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+
+          do k=1,nz
+            segment%rx_norm_rad(I,j,k) = r_norm_col(k)
+            OBC%rx_normal(I,j,k) = segment%rx_norm_rad(I,j,k)
+          enddo
+        endif
+
+        if (segment%oblique .and. (OBC%gamma_uv < 1.0)) then
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_obl(I,j,:), nz, h2, rxy_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+          segment%rx_norm_obl(I,j,:) = rxy_col(:)
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_obl(I,j,:), nz, h2, rxy_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+          segment%ry_norm_obl(I,j,:) = rxy_col(:)
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%cff_normal(I,j,:), nz, h2, rxy_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+          segment%cff_normal(I,j,:) = rxy_col(:)
+
+          do k=1,nz
+            OBC%rx_oblique_u(I,j,k) = segment%rx_norm_obl(I,j,k)
+            OBC%ry_oblique_u(I,j,k) = segment%ry_norm_obl(I,j,k)
+            OBC%cff_normal_u(I,j,k) = segment%cff_normal(I,j,k)
+          enddo
+        endif
+
+      enddo
+    elseif (segment%is_N_or_S) then
+      J = segment%HI%JsdB
+      do i=segment%HI%isd,segment%HI%ied
+
+        ! Store a column of the start and final grids
+        if (segment%direction == OBC_DIRECTION_S) then
+          if (G%mask2dT(i,j+1) == 0.0) cycle
+          h1(:) = h_old(i,j+1,:)
+          h2(:) = h_new(i,j+1,:)
+          if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i,j+1,:) ; endif
+        else
+          if (G%mask2dT(i,j) == 0.0) cycle
+          h1(:) = h_old(i,j,:)
+          h2(:) = h_new(i,j,:)
+          if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i,j,:) ; endif
+        endif
+
+        ! Vertically remap the reservoir tracer concentrations
+        do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then
+          I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale
+
+          if (present(PCM_cell)) then
+            call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column, &
+                                  h_neglect, h_neglect, PCM_cell=PCM)
+          else
+            call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column, &
+                                  h_neglect, h_neglect)
+          endif
+
+          ! Possibly underflow any very tiny tracer concentrations to 0?
+
+          ! Update tracer concentrations
+          segment%tr_Reg%Tr(m)%tres(i,J,:) = tr_column(:)
+          if (allocated(OBC%tres_y)) then ; do k=1,nz
+            OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k)
+          enddo ; endif
+
+        endif ; enddo
+
+        if (segment%radiation .and. (OBC%gamma_uv < 1.0)) then
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_rad(i,J,:), nz, h2, r_norm_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+
+          do k=1,nz
+            segment%ry_norm_rad(i,J,k) = r_norm_col(k)
+            OBC%ry_normal(i,J,k) = segment%ry_norm_rad(i,J,k)
+          enddo
+        endif
+
+        if (segment%oblique .and. (OBC%gamma_uv < 1.0)) then
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_obl(i,J,:), nz, h2, rxy_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+          segment%rx_norm_obl(i,J,:) = rxy_col(:)
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_obl(i,J,:), nz, h2, rxy_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+          segment%ry_norm_obl(i,J,:) = rxy_col(:)
+          call remapping_core_h(OBC%remap_CS, nz, h1, segment%cff_normal(i,J,:), nz, h2, rxy_col, &
+                                h_neglect, h_neglect, PCM_cell=PCM)
+          segment%cff_normal(i,J,:) = rxy_col(:)
+
+          do k=1,nz
+            OBC%rx_oblique_v(i,J,k) = segment%rx_norm_obl(i,J,k)
+            OBC%ry_oblique_v(i,J,k) = segment%ry_norm_obl(i,J,k)
+            OBC%cff_normal_v(i,J,k) = segment%cff_normal(i,J,k)
+          enddo
+        endif
+
+      enddo
+    endif
+  enddo ; endif ; endif
+
+end subroutine remap_OBC_fields
+
+
 !> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
 !!
 !! If the bottom most interface is below the topography then the bottom-most