+Add remapping of auxiliary variables in restarts

src/ALE/MOM_ALE.F90

@@ -129,6 +129,8 @@ module MOM_ALE
 public ALE_remap_scalar
 public ALE_remap_tracers
 public ALE_remap_velocities
+public ALE_remap_interface_vals
+public ALE_remap_vertex_vals
 public ALE_PLM_edge_values
 public TS_PLM_edge_values
 public TS_PPM_edge_values
@@ -289,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, &
@@ -549,12 +551,12 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
     endif
   enddo ; enddo
 
-  do j = jsc,jec ; do i=isc,iec
+  do j=jsc,jec ; do i=isc,iec
     if (G%mask2dT(i,j)>0.) then
       if (check_column_integrals(nk, h_src, nk, h_dest)) then
         call MOM_error(FATAL, "ALE_offline_inputs: Kd interpolation columns do not match")
       endif
-      call interpolate_column(nk, h(i,j,:), Kd(i,j,:), nk, h_new(i,j,:), Kd(i,j,:))
+      call interpolate_column(nk, h(i,j,:), Kd(i,j,:), nk, h_new(i,j,:), Kd(i,j,:), .true.)
     endif
   enddo ; enddo
 
@@ -974,6 +976,89 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old, h_new, u, v, OBC, dzInterface,
 
 end subroutine ALE_remap_velocities
 
+!> Interpolate to find an updated array of values at interfaces after remapping.
+subroutine ALE_remap_interface_vals(CS, G, GV, h_old, h_new, int_val)
+  type(ALE_CS),                              intent(in)    :: CS       !< ALE control structure
+  type(ocean_grid_type),                     intent(in)    :: G        !< Ocean 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]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                                             intent(inout) :: int_val  !< The interface values to interpolate [A]
+
+  real :: val_src(GV%ke+1)  ! A column of interface values on the source grid [A]
+  real :: val_tgt(GV%ke+1)  ! A column of interface values on the target grid [A]
+  real :: h_src(GV%ke)      ! A column of source grid layer thicknesses [H ~> m or kg m-2]
+  real :: h_tgt(GV%ke)      ! A column of target grid layer thicknesses [H ~> m or kg m-2]
+  integer :: i, j, k, nz
+
+  nz = GV%ke
+
+  do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j)>0.) then
+    do k=1,nz
+      h_src(k) = h_old(i,j,k)
+      h_tgt(k) = h_new(i,j,k)
+    enddo
+
+    do K=1,nz+1
+      val_src(K) = int_val(i,j,K)
+    enddo
+
+    call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)
+
+    do K=1,nz+1
+      int_val(i,j,K) = val_tgt(K)
+    enddo
+  endif ; enddo ; enddo
+
+end subroutine ALE_remap_interface_vals
+
+!> Interpolate to find an updated array of values at vertices of tracer cells after remapping.
+subroutine ALE_remap_vertex_vals(CS, G, GV, h_old, h_new, vert_val)
+  type(ALE_CS),                              intent(in)    :: CS       !< ALE control structure
+  type(ocean_grid_type),                     intent(in)    :: G        !< Ocean 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]
+  real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
+                                             intent(inout) :: vert_val  !< The interface values to interpolate [A]
+
+  real :: val_src(GV%ke+1)  ! A column of interface values on the source grid [A]
+  real :: val_tgt(GV%ke+1)  ! A column of interface values on the target grid [A]
+  real :: h_src(GV%ke)      ! A column of source grid layer thicknesses [H ~> m or kg m-2]
+  real :: h_tgt(GV%ke)      ! A column of target grid layer thicknesses [H ~> m or kg m-2]
+  real :: I_mask_sum        ! The inverse of the tracer point masks surrounding a corner [nondim]
+  integer :: i, j, k, nz
+
+  nz = GV%ke
+
+  do J=G%JscB,G%JecB ; do I=G%IscB,G%IecB
+    if ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) > 0.0 ) then
+      I_mask_sum = 1.0 / ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)))
+
+    do k=1,nz
+      h_src(k) = ((G%mask2dT(i,j) * h_old(i,j,k) + G%mask2dT(i+1,j+1) * h_old(i+1,j+1,k)) + &
+                  (G%mask2dT(i+1,j) * h_old(i+1,j,k) + G%mask2dT(i,j+1) * h_old(i,j+1,k)) ) * I_mask_sum
+      h_tgt(k) = ((G%mask2dT(i,j) * h_new(i,j,k) + G%mask2dT(i+1,j+1) * h_new(i+1,j+1,k)) + &
+                  (G%mask2dT(i+1,j) * h_new(i+1,j,k) + G%mask2dT(i,j+1) * h_new(i,j+1,k)) ) * I_mask_sum
+    enddo
+
+    do K=1,nz+1
+      val_src(K) = vert_val(I,J,K)
+    enddo
+
+    call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)
+
+    do K=1,nz+1
+      vert_val(I,J,K) = val_tgt(K)
+    enddo
+  endif ; enddo ; enddo
+
+end subroutine ALE_remap_vertex_vals
 
 !> Mask out thicknesses to 0 when their running sum exceeds a specified value.
 subroutine apply_partial_cell_mask(h1, h_mask)

src/ALE/MOM_remapping.F90

@@ -840,78 +840,72 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth
 end subroutine remap_via_sub_cells
 
 !> Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest
-subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest)
+subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, mask_edges)
   integer,                  intent(in)    :: nsrc   !< Number of source cells
   real, dimension(nsrc),    intent(in)    :: h_src  !< Thickness of source cells [H]
   real, dimension(nsrc+1),  intent(in)    :: u_src  !< Values at source cell interfaces [A]
   integer,                  intent(in)    :: ndest  !< Number of destination cells
   real, dimension(ndest),   intent(in)    :: h_dest !< Thickness of destination cells [H]
   real, dimension(ndest+1), intent(inout) :: u_dest !< Interpolated value at destination cell interfaces [A]
+  logical,                  intent(in)    :: mask_edges !< If true, mask the values outside of massless
+                                                    !! layers at the top and bottom of the column.
 
   ! Local variables
-  real :: x_dest ! Relative position of target interface [H]
-  real :: dh ! Source cell thickness [H]
-  real :: u1, u2 ! Values to interpolate between [A]
-  real :: weight_a, weight_b ! Weights for interpolation [nondim]
-  integer :: k_src, k_dest ! Index of cell in src and dest columns
-  logical :: still_vanished ! Used for figuring out what to mask as missing
-
-  ! Initial values for the loop
-  still_vanished = .true.
+  real :: x_dest            ! Relative position of target interface [H]
+  real :: dh                ! Source cell thickness [H]
+  real :: frac_pos(ndest+1) ! Fractional position of the destination interface
+                            ! within the source layer [nondim], 0 <= frac_pos <= 1.
+  integer :: k_src(ndest+1) ! Source grid layer index of destination interface, 1 <= k_src <= ndest.
+  integer :: ks, k_dest     ! Index of cell in src and dest columns
 
   ! The following forces the "do while" loop to do one cycle that will set u1, u2, dh.
-  k_src = 0
+  ks = 0
   dh = 0.
   x_dest = 0.
 
-  do k_dest=1, ndest+1
-    do while (dh<=x_dest .and. k_src<nsrc)
+  ! Find the layer index and fractional position of the interfaces of the target
+  ! grid on the source grid.
+  do k_dest=1,ndest+1
+    do while (dh<=x_dest .and. ks<nsrc)
       ! Move positions pointers forward until the interval 0 .. dh spans x_dest.
       x_dest = x_dest - dh
-      k_src = k_src + 1
-      dh = h_src(k_src) ! Thickness of layer k_src
-
-      ! Values that will be used for the interpolation.
-      u1 = u_src(k_src) ! Value on left of source cell
-      u2 = u_src(k_src+1) ! Value on right of source cell
+      ks = ks + 1
+      dh = h_src(ks) ! Thickness of layer ks
     enddo
+    k_src(k_dest) = ks
 
     if (dh>0.) then
-      weight_b = max(0., min(1., x_dest / dh)) ! Weight of u2
+      frac_pos(k_dest) = max(0., min(1., x_dest / dh)) ! Weight of u2
     else  ! For a vanished source layer we need to do something reasonable...
-      weight_b = 0.5
+      frac_pos(k_dest) = 0.5
     endif
-    weight_a = 1.0 - weight_b   ! Weight of u1
-    ! Linear interpolation between u1 and u2
-    u_dest(k_dest) = weight_a * u1 + weight_b * u2
 
-    ! Mask vanished layers at the surface which would be under an ice-shelf.
-    ! TODO: Need to figure out what to do for an isopycnal coordinate diagnostic that could
-    !       also have vanished layers at the surface.
-    if (k_dest<=ndest) then
+    if (k_dest <= ndest) then
       x_dest = x_dest + h_dest(k_dest) ! Position of interface k_dest+1
-      if (still_vanished .and. h_dest(k_dest)==0.) then
-        ! When the layer k_dest is vanished and all layers above are also vanished, the k_dest
-        ! interface value should be missing.
-        u_dest(k_dest) = 0.0
-      else
-        still_vanished = .false.
-      endif
     endif
+  enddo
 
+  do k_dest=1,ndest+1
+    ! Linear interpolation between surrounding edge values.
+    ks = k_src(k_dest)
+    u_dest(k_dest) = (1.0 - frac_pos(k_dest)) * u_src(ks) + frac_pos(k_dest) * u_src(ks+1)
   enddo
 
-  ! Mask vanished layers on topography
-  still_vanished = .true.
-  do k_dest=ndest, 1, -1
-    if (still_vanished .and. h_dest(k_dest)==0.) then
-      ! When the layer k_dest is vanished and all layers below are also vanished, the k_dest+1
-      ! interface value should be missing.
+  if (mask_edges) then
+    ! Mask vanished layers at the surface which would be under an ice-shelf.
+    ! When the layer k_dest is vanished and all layers above are also vanished,
+    ! the k_dest interface value should be missing.
+    do k_dest=1,ndest
+      if (h_dest(k_dest) > 0.) exit
+      u_dest(k_dest) = 0.0
+    enddo
+
+    ! Mask interfaces below vanished layers at the bottom
+    do k_dest=ndest,1,-1
+      if (h_dest(k_dest) > 0.) exit
       u_dest(k_dest+1) = 0.0
-    else
-      exit
-    endif
-  enddo
+    enddo
+  endif
 
 end subroutine interpolate_column
 
@@ -1717,7 +1711,7 @@ logical function test_interp(verbose, msg, nsrc, h_src, u_src, ndest, h_dest, u_
   real :: error
 
   ! Interpolate from src to dest
-  call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest)
+  call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, .true.)
 
   test_interp = .false.
   do k=1,ndest+1

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 "//&