Merge remote-tracking branch 'gfdl/dev/gfdl' into esmg_work

config_src/infra/FMS1/MOM_domain_infra.F90

@@ -19,7 +19,8 @@ module MOM_domain_infra
 use mpp_domains_mod, only : mpp_compute_block_extent
 use mpp_domains_mod, only : mpp_broadcast_domain, mpp_redistribute, mpp_global_field
 use mpp_domains_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM
-use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE
+use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN
+use mpp_domains_mod, only : FOLD_NORTH_EDGE, FOLD_SOUTH_EDGE, FOLD_EAST_EDGE, FOLD_WEST_EDGE
 use mpp_domains_mod, only : To_East => WUPDATE, To_West => EUPDATE, Omit_Corners => EDGEUPDATE
 use mpp_domains_mod, only : To_North => SUPDATE, To_South => NUPDATE
 use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE => NORTH, EAST_FACE => EAST
@@ -1553,6 +1554,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain
     call get_layout_extents(MD_in, exnj, exni)
 
     MOM_dom%X_FLAGS = MD_in%Y_FLAGS ; MOM_dom%Y_FLAGS = MD_in%X_FLAGS
+    ! Correct the position of a tripolar grid, assuming that flags are not additive.
+    if (qturns == 1) then
+      if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE
+      if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE
+      if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE
+      if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE
+    elseif (qturns == 3) then
+      if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE
+      if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE
+      if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE
+      if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE
+    endif
+
     MOM_dom%layout(:) = MD_in%layout(2:1:-1)
     MOM_dom%io_layout(:) = io_layout_in(2:1:-1)
   else
@@ -1561,11 +1575,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain
     call get_layout_extents(MD_in, exni, exnj)
 
     MOM_dom%X_FLAGS = MD_in%X_FLAGS ; MOM_dom%Y_FLAGS = MD_in%Y_FLAGS
+    ! Correct the position of a tripolar grid, assuming that flags are not additive.
+    if (qturns == 2) then
+      if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE
+      if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE
+      if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE
+      if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE
+    endif
+
     MOM_dom%layout(:) = MD_in%layout(:)
     MOM_dom%io_layout(:) = io_layout_in(:)
   endif
 
-  ! Ensure that the points per processor are the same on the source and densitation grids.
+  ! Ensure that the points per processor are the same on the source and destination grids.
   select case (qturns)
     case (1) ; call invert(exni)
     case (2) ; call invert(exni) ; call invert(exnj)

config_src/infra/FMS2/MOM_domain_infra.F90

@@ -19,7 +19,8 @@ module MOM_domain_infra
 use mpp_domains_mod, only : mpp_compute_block_extent
 use mpp_domains_mod, only : mpp_broadcast_domain, mpp_redistribute, mpp_global_field
 use mpp_domains_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM
-use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE
+use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN
+use mpp_domains_mod, only : FOLD_NORTH_EDGE, FOLD_SOUTH_EDGE, FOLD_EAST_EDGE, FOLD_WEST_EDGE
 use mpp_domains_mod, only : To_East => WUPDATE, To_West => EUPDATE, Omit_Corners => EDGEUPDATE
 use mpp_domains_mod, only : To_North => SUPDATE, To_South => NUPDATE
 use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE => NORTH, EAST_FACE => EAST
@@ -1555,6 +1556,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain
     call get_layout_extents(MD_in, exnj, exni)
 
     MOM_dom%X_FLAGS = MD_in%Y_FLAGS ; MOM_dom%Y_FLAGS = MD_in%X_FLAGS
+    ! Correct the position of a tripolar grid, assuming that flags are not additive.
+    if (modulo(qturns, 4) == 1) then
+      if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE
+      if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE
+      if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE
+      if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE
+    elseif (modulo(qturns, 4) == 3) then
+      if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE
+      if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE
+      if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE
+      if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE
+    endif
+
     MOM_dom%layout(:) = MD_in%layout(2:1:-1)
     MOM_dom%io_layout(:) = io_layout_in(2:1:-1)
   else
@@ -1563,11 +1577,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain
     call get_layout_extents(MD_in, exni, exnj)
 
     MOM_dom%X_FLAGS = MD_in%X_FLAGS ; MOM_dom%Y_FLAGS = MD_in%Y_FLAGS
+    ! Correct the position of a tripolar grid, assuming that flags are not additive.
+    if (modulo(qturns, 4) == 2) then
+      if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE
+      if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE
+      if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE
+      if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE
+    endif
+
     MOM_dom%layout(:) = MD_in%layout(:)
     MOM_dom%io_layout(:) = io_layout_in(:)
   endif
 
-  ! Ensure that the points per processor are the same on the source and densitation grids.
+  ! Ensure that the points per processor are the same on the source and destination grids.
   select case (qturns)
     case (1) ; call invert(exni)
     case (2) ; call invert(exni) ; call invert(exnj)

src/core/MOM.F90

@@ -2904,6 +2904,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   if (CS%rotate_index) then
     G_in%ke = GV%ke
 
+    ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
+    if (CS%debug .or. G_in%symmetric) then
+      call clone_MOM_domain(G_in%Domain, G_in%Domain_aux, symmetric=.false.)
+    else ; G_in%Domain_aux => G_in%Domain ; endif
+
     allocate(u_in(G_in%IsdB:G_in%IedB, G_in%jsd:G_in%jed, nz), source=0.0)
     allocate(v_in(G_in%isd:G_in%ied, G_in%JsdB:G_in%JedB, nz), source=0.0)
     allocate(h_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz), source=GV%Angstrom_H)

src/core/MOM_dynamics_split_RK2.F90

@@ -851,7 +851,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
   call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, &
                             MEKE, Varmix, G, GV, US, CS%hor_visc, &
                             OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
-                            ADp=CS%ADp)
+                            ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v)
   call cpu_clock_end(id_clock_horvisc)
   if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
 
@@ -1518,7 +1518,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. &
       .not. query_initialized(CS%diffv, "diffv", restart_CS)) then
     call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, G, GV, US, CS%hor_visc, &
-                              OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp)
+                              OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
+                              hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v)
     call set_initialized(CS%diffu, "diffu", restart_CS)
     call set_initialized(CS%diffv, "diffv", restart_CS)
   endif

src/core/MOM_forcing_type.F90

@@ -6,6 +6,7 @@ module MOM_forcing_type
 use MOM_array_transform, only : rotate_array, rotate_vector, rotate_array_pair
 use MOM_coupler_types, only : coupler_2d_bc_type, coupler_type_destructor
 use MOM_coupler_types, only : coupler_type_increment_data, coupler_type_initialized
+use MOM_coupler_types, only : coupler_type_copy_data
 use MOM_cpu_clock,     only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
 use MOM_debugging,     only : hchksum, uvchksum
 use MOM_diag_mediator, only : post_data, register_diag_field, register_scalar_field
@@ -3702,9 +3703,11 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns)
   if (associated(fluxes_in%ustar_tidal)) &
     call rotate_array(fluxes_in%ustar_tidal, turns, fluxes%ustar_tidal)
 
-  ! TODO: tracer flux rotation
-  if (coupler_type_initialized(fluxes%tr_fluxes)) &
-    call MOM_error(FATAL, "Rotation of tracer BC fluxes not yet implemented.")
+  ! NOTE: Tracer fields are handled by FMS, so are left unrotated.  Any
+  ! reads/writes to tr_fields must be appropriately rotated.
+  if (coupler_type_initialized(fluxes%tr_fluxes)) then
+    call coupler_type_copy_data(fluxes_in%tr_fluxes, fluxes%tr_fluxes)
+  endif
 
   ! Scalars and flags
   fluxes%accumulate_p_surf = fluxes_in%accumulate_p_surf
@@ -3757,10 +3760,18 @@ subroutine rotate_mech_forcing(forces_in, turns, forces)
   endif
 
   if (do_press) then
-    ! NOTE: p_surf_SSH either points to p_surf or p_surf_full
     call rotate_array(forces_in%p_surf, turns, forces%p_surf)
     call rotate_array(forces_in%p_surf_full, turns, forces%p_surf_full)
     call rotate_array(forces_in%net_mass_src, turns, forces%net_mass_src)
+
+    ! p_surf_SSH points to either p_surf or p_surf_full
+    if (associated(forces_in%p_surf_SSH, forces_in%p_surf)) then
+      forces%p_surf_SSH => forces%p_surf
+    else if (associated(forces_in%p_surf_SSH, forces_in%p_surf_full)) then
+      forces%p_surf_SSH => forces%p_surf_full
+    else
+      forces%p_surf_SSH => null()
+    endif
   endif
 
   if (do_iceberg) then

src/core/MOM_porous_barriers.F90

@@ -168,9 +168,10 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt)
 
   if (CS%debug) then
     call uvchksum("Interface height used by porous barrier for layer weights", &
-                  eta_u, eta_v, G%HI, haloshift=0)
+                  eta_u, eta_v, G%HI, haloshift=0, scalar_pair=.true.)
     call uvchksum("Porous barrier layer-averaged weights: por_face_area[UV]", &
-                  pbv%por_face_areaU, pbv%por_face_areaV, G%HI, haloshift=0)
+                  pbv%por_face_areaU, pbv%por_face_areaV, G%HI, haloshift=0, &
+                  scalar_pair=.true.)
   endif
 
   if (CS%id_por_face_areaU > 0) call post_data(CS%id_por_face_areaU, pbv%por_face_areaU, CS%diag)
@@ -256,9 +257,10 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt)
 
   if (CS%debug) then
     call uvchksum("Interface height used by porous barrier for interface weights", &
-                  eta_u, eta_v, G%HI, haloshift=0)
+                  eta_u, eta_v, G%HI, haloshift=0, scalar_pair=.true.)
     call uvchksum("Porous barrier weights at the layer-interface: por_layer_width[UV]", &
-                  pbv%por_layer_widthU, pbv%por_layer_widthV, G%HI, haloshift=0)
+                  pbv%por_layer_widthU, pbv%por_layer_widthV, G%HI, &
+                  haloshift=0, scalar_pair=.true.)
   endif
 
   if (CS%id_por_layer_widthU > 0) call post_data(CS%id_por_layer_widthU, pbv%por_layer_widthU, CS%diag)

src/core/MOM_variables.F90

@@ -6,6 +6,7 @@ module MOM_variables
 use MOM_array_transform, only : rotate_array, rotate_vector
 use MOM_coupler_types, only : coupler_1d_bc_type, coupler_2d_bc_type
 use MOM_coupler_types, only : coupler_type_spawn, coupler_type_destructor, coupler_type_initialized
+use MOM_coupler_types, only : coupler_type_copy_data
 use MOM_debugging,     only : hchksum
 use MOM_domains,       only : MOM_domain_type, get_domain_extent, group_pass_type
 use MOM_EOS,           only : EOS_type
@@ -499,9 +500,11 @@ subroutine rotate_surface_state(sfc_state_in, sfc_state, G, turns)
   sfc_state%T_is_conT = sfc_state_in%T_is_conT
   sfc_state%S_is_absS = sfc_state_in%S_is_absS
 
-  ! TODO: tracer field rotation
-  if (coupler_type_initialized(sfc_state_in%tr_fields)) &
-    call MOM_error(FATAL, "Rotation of surface state tracers is not yet implemented.")
+  ! NOTE: Tracer fields are handled by FMS, so are left unrotated.  Any
+  ! reads/writes to tr_fields must be appropriately rotated.
+  if (coupler_type_initialized(sfc_state_in%tr_fields)) then
+    call coupler_type_copy_data(sfc_state_in%tr_fields, sfc_state%tr_fields)
+  endif
 end subroutine rotate_surface_state
 
 !> Allocates the arrays contained within a BT_cont_type and initializes them to 0.

src/framework/MOM_io_file.F90

@@ -34,6 +34,7 @@ module MOM_io_file
 use MOM_netcdf, only : write_netcdf_attribute
 use MOM_netcdf, only : get_netcdf_size
 use MOM_netcdf, only : get_netcdf_fields
+use MOM_netcdf, only : get_netcdf_filename
 use MOM_netcdf, only : read_netcdf_field
 
 use MOM_error_handler, only : MOM_error, FATAL
@@ -1326,7 +1327,13 @@ subroutine open_file_nc(handle, filename, action, MOM_domain, threading, fileset
 
   if (present(MOM_domain)) then
     handle%domain_decomposed = .true.
-    call hor_index_init(MOM_domain, handle%HI)
+
+    ! Input files use unrotated indexing.
+    if (associated(MOM_domain%domain_in)) then
+      call hor_index_init(MOM_domain%domain_in, handle%HI)
+    else
+      call hor_index_init(MOM_domain, handle%HI)
+    endif
   endif
 
   call handle%axes%init()
@@ -1751,8 +1758,9 @@ subroutine get_field_nc(handle, label, values, rescale)
   ! NOTE: Data on face and vertex points is not yet supported.  This is a
   ! temporary check to detect such cases, but may be removed in the future.
   if (.not. (compute_domain .or. data_domain)) &
-    call MOM_error(FATAL, 'get_field_nc: Only compute and data domains ' // &
-        'are currently supported.')
+    call MOM_error(FATAL, 'get_field_nc trying to read '//trim(label)//' from '//&
+                   trim(get_netcdf_filename(handle%handle_nc))//&
+                   ': Only compute and data domains are currently supported.')
 
   field_nc = handle%fields%get(label)
 

src/framework/MOM_netcdf.F90

@@ -39,6 +39,7 @@ module MOM_netcdf
 public :: write_netcdf_attribute
 public :: get_netcdf_size
 public :: get_netcdf_fields
+public :: get_netcdf_filename
 public :: read_netcdf_field
 
 
@@ -722,6 +723,14 @@ subroutine get_netcdf_fields(handle, axes, fields)
   fields(:) = vars(:nfields)
 end subroutine get_netcdf_fields
 
+!> Return the name of a file from a netCDF handle
+function get_netcdf_filename(handle)
+  type(netcdf_file_type), intent(in) :: handle !< A netCDF file handle
+  character(len=:), allocatable :: get_netcdf_filename !< The name of the file that this handle refers to.
+
+  get_netcdf_filename = handle%filename
+
+end function
 
 !> Read the values of a field from a netCDF file
 subroutine read_netcdf_field(handle, field, values, bounds)

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -111,7 +111,7 @@ module MOM_hor_visc
                              !! limit the grid Reynolds number [L2 T-1 ~> m2 s-1]
   real    :: min_grid_Ah     !< Minimun horizontal biharmonic viscosity used to
                              !! limit grid Reynolds number [L4 T-1 ~> m4 s-1]
-
+  logical :: use_cont_thick  !< If true, thickness at velocity points adopts h[uv] in BT_cont from continuity solver.
   type(ZB2020_CS) :: ZB2020  !< Zanna-Bolton 2020 control structure.
   logical :: use_ZB2020      !< If true, use Zanna-Bolton 2020 parameterization.
 
@@ -239,7 +239,7 @@ module MOM_hor_visc
 !!   v[is-2:ie+2,js-2:je+2]
 !!   h[is-1:ie+1,js-1:je+1]
 subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
-                                CS, OBC, BT, TD, ADp)
+                                CS, OBC, BT, TD, ADp, hu_cont, hv_cont)
   type(ocean_grid_type),         intent(in)  :: G      !< The ocean's grid structure.
   type(verticalGrid_type),       intent(in)  :: GV     !< The ocean's vertical grid structure.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
@@ -263,6 +263,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   type(barotropic_CS), intent(in), optional  :: BT     !< Barotropic control structure
   type(thickness_diffuse_CS), intent(in), optional :: TD  !< Thickness diffusion control structure
   type(accel_diag_ptrs), intent(in), optional :: ADp   !< Acceleration diagnostics
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
+                        optional, intent(in) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
+                        optional, intent(in) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2].
 
   ! Local variables
   real, dimension(SZIB_(G),SZJ_(G)) :: &
@@ -391,6 +395,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   logical :: apply_OBC = .false.
   logical :: use_MEKE_Ku
   logical :: use_MEKE_Au
+  logical :: use_cont_huv
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
   integer :: i, j, k, n
   real :: inv_PI3, inv_PI2, inv_PI6 ! Powers of the inverse of pi [nondim]
@@ -445,6 +450,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   use_MEKE_Ku = allocated(MEKE%Ku)
   use_MEKE_Au = allocated(MEKE%Au)
 
+  use_cont_huv = CS%use_cont_thick .and. present(hu_cont) .and. present(hv_cont)
+
   rescale_Kh = .false.
   if (VarMix%use_variable_mixing) then
     rescale_Kh = VarMix%Resoln_scaled_Kh
@@ -554,12 +561,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   !$OMP   CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, &
   !$OMP   is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
   !$OMP   apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
-  !$OMP   use_MEKE_Ku, use_MEKE_Au, &
+  !$OMP   use_MEKE_Ku, use_MEKE_Au, use_cont_huv, &
   !$OMP   backscat_subround, GME_effic_h, GME_effic_q, &
   !$OMP   h_neglect, h_neglect3, inv_PI3, inv_PI6, &
   !$OMP   diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, &
   !$OMP   div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, &
-  !$OMP   KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt &
+  !$OMP   KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, hu_cont, hv_cont &
   !$OMP ) &
   !$OMP private( &
   !$OMP   i, j, k, n, &
@@ -658,6 +665,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       enddo ; enddo
     endif
 
+    ! The following should obviously be combined with the previous block if adopted.
+    if (use_cont_huv) then
+      do j=js-2,je+2 ; do I=Isq-1,Ieq+1
+        h_u(I,j) = hu_cont(I,j,k)
+      enddo ; enddo
+      do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2
+        h_v(i,J) = hv_cont(i,J,k)
+      enddo ; enddo
+    endif
+
     ! Adjust contributions to shearing strain and interpolated values of
     ! thicknesses on open boundaries.
     if (apply_OBC) then ; do n=1,OBC%number_of_segments
@@ -1969,6 +1986,9 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
   if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
 
   call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)
+  call get_param(param_file, mdl, "USE_CONT_THICKNESS", CS%use_cont_thick, &
+                 "If true, use thickness at velocity points from continuity solver. This option"//&
+                 "currently only works with split mode.", default=.false.)
   call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, &
                  "If true, use a Laplacian horizontal viscosity.", &
                  default=.false.)