diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90
index d01e7f9639..9c59a4768a 100644
--- a/src/core/MOM_checksum_packages.F90
+++ b/src/core/MOM_checksum_packages.F90
@@ -52,16 +52,23 @@ module MOM_checksum_packages
 ! =============================================================================
 
 subroutine MOM_state_chksum_5arg(mesg, u, v, h, uh, vh, G, GV, haloshift, symmetric)
-  character(len=*),                          intent(in) :: mesg
-  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_(G)), intent(in) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uh
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vh
-  integer, optional,                         intent(in) :: haloshift
-  logical, optional,                         intent(in) :: symmetric
+  character(len=*),                          &
+                           intent(in) :: mesg !< A message that appears on the chksum lines.
+  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_(G)), &
+                           intent(in) :: u    !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                           intent(in) :: v    !< The meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                           intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2).
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                           intent(in) :: uh   !< Volume flux through zonal faces = u*h*dy, m3 s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                           intent(in) :: vh   !< Volume flux through meridional
+                                              !! faces = v*h*dx, in m3 s-1.
+  integer, optional,       intent(in) :: haloshift
+  logical, optional,       intent(in) :: symmetric
 !   This subroutine writes out chksums for the model's basic state variables.
 ! Arguments: mesg - A message that appears on the chksum lines.
 !  (in)      u - Zonal velocity, in m s-1.
@@ -89,14 +96,17 @@ end subroutine MOM_state_chksum_5arg
 ! =============================================================================
 
 subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, haloshift, symmetric)
-  character(len=*),                          intent(in) :: mesg
-  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_(G)), intent(in) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  integer, optional,                         intent(in) :: haloshift
-  logical, optional,                         intent(in) :: symmetric
+  character(len=*),        intent(in) :: mesg !< A message that appears on the chksum lines.
+  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_(G)), &
+                           intent(in) :: u    !< Zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                           intent(in) :: v    !< Meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                           intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2).
+  integer, optional,       intent(in) :: haloshift
+  logical, optional,       intent(in) :: symmetric
 !   This subroutine writes out chksums for the model's basic state variables.
 ! Arguments: mesg - A message that appears on the chksum lines.
 !  (in)      u - Zonal velocity, in m s-1.
@@ -122,9 +132,10 @@ end subroutine MOM_state_chksum_3arg
 ! =============================================================================
 
 subroutine MOM_thermo_chksum(mesg, tv, G, haloshift)
-  character(len=*),         intent(in) :: mesg
-  type(thermo_var_ptrs),    intent(in) :: tv   !< A structure pointing to various thermodynamic variables
-  type(ocean_grid_type),    intent(in) :: G    !< The ocean's grid structure
+  character(len=*),         intent(in) :: mesg !< A message that appears on the chksum lines.
+  type(thermo_var_ptrs),    intent(in) :: tv   !< A structure pointing to various
+                                               !! thermodynamic variables.
+  type(ocean_grid_type),    intent(in) :: G    !< The ocean's grid structure.
   integer, optional,        intent(in) :: haloshift
 !   This subroutine writes out chksums for the model's thermodynamic state
 ! variables.
@@ -147,19 +158,39 @@ end subroutine MOM_thermo_chksum
 
 subroutine MOM_accel_chksum(mesg, CAu, CAv, PFu, PFv, diffu, diffv, G, GV, pbce, &
                             u_accel_bt, v_accel_bt, symmetric)
-  character(len=*),                          intent(in) :: mesg
-  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_(G)), intent(in) :: CAu
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: CAv
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: PFu
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: PFv
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: diffu
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: diffv
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  optional, intent(in) :: pbce
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: u_accel_bt
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), optional, intent(in) :: v_accel_bt
-  logical, optional,                         intent(in) :: symmetric
+  character(len=*),         intent(in) :: mesg !< A message that appears on the chksum lines.
+  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_(G)), &
+                            intent(in) :: CAu  !< Zonal acceleration due to Coriolis
+                                               !! and momentum advection terms, in m s-2.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                            intent(in) :: CAv  !< Meridional acceleration due to Coriolis
+                                               !! and momentum advection terms, in m s-2.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                            intent(in) :: PFu  !< Zonal acceleration due to pressure gradients
+                                               !! (equal to -dM/dx) in m s-2.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                            intent(in) :: PFv  !< Meridional acceleration due to pressure gradients
+                                               !! (equal to -dM/dy) in m s-2.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                            intent(in) :: diffu !< Zonal acceleration due to convergence of the
+                                                !! along-isopycnal stress tensor, in m s-2.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                            intent(in) :: diffv !< Meridional acceleration due to convergence of
+                                                !! the along-isopycnal stress tensor, in m s-2.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                  optional, intent(in) :: pbce !< The baroclinic pressure anomaly in each layer
+                                               !! due to free surface height anomalies, in
+                                               !! m s-2.pbce points to a space with nz layers
+                                               !! or NULL.                                                                            !! NULL.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                  optional, intent(in) :: u_accel_bt !< The zonal acceleration from terms in the
+                                                     !! barotropic solver,in m s-2.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                  optional, intent(in) :: v_accel_bt !< The meridional acceleration from terms in
+                                                     !! the barotropic solver,in m s-2.
+  logical, optional,        intent(in) :: symmetric
 
 !   This subroutine writes out chksums for the model's accelerations.
 ! Arguments: mesg - A message that appears on the chksum lines.
@@ -205,13 +236,24 @@ end subroutine MOM_accel_chksum
 ! =============================================================================
 
 subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, allowChange, permitDiminishing)
-  type(ocean_grid_type),                     intent(in) :: G    !< The ocean's grid structure
-  character(len=*),                          intent(in) :: mesg
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  real, pointer, dimension(:,:,:),           intent(in) :: Temp, Salt
-  logical, optional,                         intent(in) :: allowChange, permitDiminishing
+  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
+  character(len=*),        intent(in) :: mesg !< A message that appears on the chksum lines.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                           intent(in) :: u    !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                           intent(in) :: v    !< The meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                           intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2).
+  real, pointer, dimension(:,:,:),           &
+                           intent(in) :: Temp !< Temperature in degree C.
+  real, pointer, dimension(:,:,:),           &
+                           intent(in) :: Salt !< Salinity, in ppt.
+
+  logical, optional,       intent(in) :: allowChange !< do not flag an error
+                                                                       !! if the statistics change.
+  logical, optional,                         &
+                           intent(in) :: permitDiminishing !< do not flag error
+                                                           !!if the extrema are diminishing.
 !   This subroutine monitors statistics for the model's state variables.
 ! Arguments: mesg - A message that appears on the chksum lines.
 !  (in) u - Zonal velocity, in m s-1.
diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90
index 9332e5cc5a..b5a8409804 100644
--- a/src/core/MOM_dynamics_legacy_split.F90
+++ b/src/core/MOM_dynamics_legacy_split.F90
@@ -272,27 +272,60 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
                  Time_local, dt, fluxes, p_surf_begin, p_surf_end, &
                  dt_since_flux, dt_therm, uh, vh, uhtr, vhtr, eta_av, &
                  G, GV, CS, calc_dtbt, VarMix, MEKE)
-  type(ocean_grid_type),                     intent(inout) :: 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_(G)), target, intent(inout) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  type(thermo_var_ptrs),                     intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
-  type(vertvisc_type),                       intent(inout) :: visc
-  type(time_type),                           intent(in)    :: Time_local
-  real,                                      intent(in)    :: dt   !< The baroclinic dynamics time step, in s
-  type(forcing),                             intent(in)    :: fluxes
-  real, dimension(:,:),                      pointer       :: p_surf_begin, p_surf_end
-  real,                                      intent(in)    :: dt_since_flux, dt_therm
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: uh
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: vh
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr
-  real, dimension(SZI_(G),SZJ_(G)),          intent(out)   :: eta_av
-  type(MOM_dyn_legacy_split_CS),             pointer       :: CS
-  logical,                                   intent(in)    :: calc_dtbt
-  type(VarMix_CS),                           pointer       :: VarMix
-  type(MEKE_type),                           pointer       :: MEKE
+  type(ocean_grid_type),      intent(inout) :: 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_(G)), &
+                      target, intent(inout) :: u    !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                      target, intent(inout) :: v    !< The meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                              intent(inout) :: h   !< Layer thicknesses, in H (usually m or kg m-2).
+  type(thermo_var_ptrs),      intent(in)    :: tv   !< A structure pointing to various.
+                                                    !! thermodynamic variables.
+  type(vertvisc_type),        intent(inout) :: visc !< A structure containing vertical viscosities,
+                                                    !! bottom drag viscosities, and related fields.
+  type(time_type),            intent(in)    :: Time_local !< The model time at the end
+                                                          !! of the time step.
+  real,                       intent(in)    :: dt   !< The baroclinic dynamics time step, in s
+  type(forcing),              intent(in)    :: fluxes !< A structure containing pointers to any
+                                                      !! possible forcing fields.  Unused fields
+                                                      !! have NULL ptrs.
+  real, dimension(:,:),       pointer       :: p_surf_begin !< A pointer (perhaps NULL) to the
+                                                            !! surface pressure at the beginning
+                                                            !! of this dynamic step, in Pa.
+  real, dimension(:,:),       pointer       :: p_surf_end !< A pointer (perhaps NULL) to the
+                                                          !! surface pressure at the end of
+                                                          !! this dynamic step, in Pa.
+  real,                       intent(in)    :: dt_since_flux !< The elapsed time since fluxes
+                                                             !! were applied, in s.
+  real,                       intent(in)    :: dt_therm !< The thermodynamic time step, in s.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                      target, intent(inout) :: uh !< The zonal volume or mass transport,
+                                                  !! in m3 s-1 or kg s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                      target, intent(inout) :: vh !< The meridional volume or mass transport,
+                                                  !! in m3 s-1 or kg s-1.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                              intent(inout) :: uhtr !< The accumulated zonal volume or mass
+                                                    !! transport since the last tracer advection,
+                                                    !! in m3 or kg.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                              intent(inout) :: vhtr !< The accumulated meridional volume or mass
+                                                    !! transport since the last tracer advection,
+                                                    !! in m3 or kg.
+  real, dimension(SZI_(G),SZJ_(G)),          &
+                              intent(out)   :: eta_av !< The free surface height or column mass
+                                                      !! time-averaged over a time step,
+                                                      !! in m or kg m-2.
+  type(MOM_dyn_legacy_split_CS),             &
+                              pointer       :: CS !< The control structure set up by
+                                                  !! initialize_dyn_legacy_split.
+  logical,                    intent(in)    :: calc_dtbt !< If true, recalculate the
+                                                         !! barotropic time step.
+  type(VarMix_CS),            pointer       :: VarMix !<A pointer to a structure with fields that
+                                                      !! specify the spatially variable viscosities.
+  type(MEKE_type),            pointer       :: MEKE   !< A pointer to a structure containing fields
+                                                    !! related to the Mesoscale Eddy Kinetic Energy.
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (inout)   v - The meridional velocity, in m s-1.
 !  (inout)   h - The layer thicknesses, in m or kg m-2, depending on
@@ -1008,13 +1041,18 @@ end subroutine step_MOM_dyn_legacy_split
 ! =============================================================================
 
 subroutine adjustments_dyn_legacy_split(u, v, h, dt, G, GV, CS)
-  type(ocean_grid_type),                     intent(inout) :: 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_(G)), intent(inout) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  real,                                      intent(in)    :: dt   !< The baroclinic dynamics time step, in s
-  type(MOM_dyn_legacy_split_CS),             pointer       :: CS
+  type(ocean_grid_type),          intent(inout) :: 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_(G)), &
+                                  intent(inout) :: u    !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                  intent(inout) :: v    !< The meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                                  intent(inout) :: h    !< Layer thicknesses, in H
+                                                        !! (usually m or kg m-2).
+  real,                           intent(in)    :: dt   !< The baroclinic dynamics time step, in s.
+  type(MOM_dyn_legacy_split_CS),  pointer       :: CS   !< The control structure set up by
+                                                        !! initialize_dyn_legacy_split.
 
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (in)      v - The meridional velocity, in m s-1.
@@ -1062,13 +1100,20 @@ end subroutine adjustments_dyn_legacy_split
 ! =============================================================================
 
 subroutine register_restarts_dyn_legacy_split(HI, GV, param_file, CS, restart_CS, uh, vh)
-  type(hor_index_type),          intent(in)    :: HI
-  type(verticalGrid_type),       intent(in)    :: GV   !< The ocean's vertical grid structure
-  type(param_file_type),         intent(in)    :: param_file !< A structure to parse for run-time parameters
-  type(MOM_dyn_legacy_split_CS), pointer       :: CS
-  type(MOM_restart_CS),          pointer       :: restart_CS
-  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), target, intent(inout) :: uh
-  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), target, intent(inout) :: vh
+  type(hor_index_type),          intent(in)   :: HI         !< A horizontal index type structure.
+  type(verticalGrid_type),       intent(in)   :: GV         !< The ocean's vertical grid structure.
+  type(param_file_type),         intent(in)   :: param_file !< A structure to parse for run-time
+                                                            !! parameters.
+  type(MOM_dyn_legacy_split_CS), pointer      :: CS         !< The control structure set up by
+                                                            !! initialize_dyn_legacy_split.
+  type(MOM_restart_CS),          pointer      :: restart_CS !< A pointer to the restart
+                                                            !! control structure.
+  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
+                        target, intent(inout) :: uh         !< The zonal volume or mass transport,
+                                                            !! in m3 s-1 or kg s-1.
+  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
+                        target, intent(inout) :: vh         !< The meridional volume or mass
+                                                            !! transport, in m3 s-1 or kg s-1.
 !   This subroutine sets up any auxiliary restart variables that are specific
 ! to the unsplit time stepping scheme.  All variables registered here should
 ! have the ability to be recreated if they are not present in a restart file.
@@ -1166,32 +1211,64 @@ subroutine initialize_dyn_legacy_split(u, v, h, uh, vh, eta, Time, G, GV, param_
                       diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
                       VarMix, MEKE, OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, visc, &
                       dirs, ntrunc)
-  type(ocean_grid_type),                     intent(inout) :: 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_(G)), intent(inout) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: uh
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: vh
-  real, dimension(SZI_(G),SZJ_(G)),          intent(inout) :: eta
-  type(time_type),                   target, intent(in)    :: Time
-  type(param_file_type),                     intent(in)    :: param_file !< A structure to parse for run-time parameters
-  type(diag_ctrl),                   target, intent(inout) :: diag
-  type(MOM_dyn_legacy_split_CS),             pointer       :: CS
-  type(MOM_restart_CS),                      pointer       :: restart_CS
-  real,                                      intent(in)    :: dt   !< The baroclinic dynamics time step, in s
-  type(accel_diag_ptrs),             target, intent(inout) :: Accel_diag
-  type(cont_diag_ptrs),              target, intent(inout) :: Cont_diag
-  type(ocean_internal_state),                intent(inout) :: MIS
-  type(VarMix_CS),                           pointer       :: VarMix
-  type(MEKE_type),                           pointer       :: MEKE
-  type(ocean_OBC_type),                      pointer       :: OBC
-  type(update_OBC_CS),                       pointer       :: update_OBC_CSp
-  type(ALE_CS),                              pointer       :: ALE_CSp
-  type(set_visc_CS),                         pointer       :: setVisc_CSp
-  type(vertvisc_type),                       intent(inout) :: visc
-  type(directories),                         intent(in)    :: dirs
-  integer, target,                           intent(inout) :: ntrunc
+  type(ocean_grid_type),         intent(inout) :: 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_(G)), &
+                                 intent(inout) :: u    !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                 intent(inout) :: v    !< The meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , &
+                                 intent(inout) :: h    !< Layer thicknesses, in H
+                                                       !! (usually m or kg m-2).
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                         target, intent(inout) :: uh   !< The zonal volume or mass transport,
+                                                       !! in m3 s-1 or kg s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                         target, intent(inout) :: vh   !< The meridional volume or mass transport,
+                                                       !! in m3 s-1 or kg s-1.
+  real, dimension(SZI_(G),SZJ_(G)),          &
+                                 intent(inout) :: eta  !< The free surface height or column mass,
+                                                       !! in m or kg m-2.
+  type(time_type),       target, intent(in)    :: Time        !< The current model time.
+  type(param_file_type),         intent(in)    :: param_file  !< A structure to parse for run-time
+                                                              !! parameters.
+  type(diag_ctrl),       target, intent(inout) :: diag        !< A structure that is used to
+                                                              !! regulate diagnostic output.
+  type(MOM_dyn_legacy_split_CS), pointer       :: CS          !< The control structure set up by
+                                                              !! initialize_dyn_legacy_split.
+  type(MOM_restart_CS),          pointer       :: restart_CS  !< A pointer to the restart control
+                                                              !! structure.
+  real,                          intent(in)    :: dt          !< The baroclinic dynamics time step,
+                                                              !! in s.
+  type(accel_diag_ptrs), target, intent(inout) :: Accel_diag  !<A set of pointers to the various
+                                        !! accelerations in the momentum equations, which can be
+                                        !! used for later derived diagnostics, like energy budgets.
+  type(cont_diag_ptrs),  target, intent(inout) :: Cont_diag  !< A structure with pointers to various
+                                                              !! terms in the continuity equations.
+  type(ocean_internal_state),    intent(inout) :: MIS         !< The "MOM6 Internal State"
+                                              !! structure, used to pass around pointers to
+                                              !! various arrays for diagnostic purposes.
+  type(VarMix_CS),               pointer       :: VarMix      !< A pointer to a structure with
+                                     !! fields that specify the spatially variable viscosities.
+  type(MEKE_type),               pointer       :: MEKE        !< A pointer to a structure containing
+                                              !! fields related to the Mesoscale Eddy Kinetic Energy.
+  type(ocean_OBC_type),          pointer       :: OBC         !< If open boundary conditions are
+                                                          !! used, this points to the ocean_OBC_type
+                                                          !! that was set up in MOM_initialization.
+  type(update_OBC_CS),           pointer       :: update_OBC_CSp !< If open boundary condition
+                                                              !! updates are used, this points to
+                                                              !! the appropriate control structure.
+  type(ALE_CS),                  pointer       :: ALE_CSp     !< This points to the ALE control
+                                                              !! structure.
+  type(set_visc_CS),             pointer       :: setVisc_CSp !< This points to the set_visc control
+                                                              !! structure.
+  type(vertvisc_type),           intent(inout) :: visc        !< A structure containing vertical
+                                         !! viscosities, bottom drag viscosities, and related fields.
+  type(directories),             intent(in)    :: dirs        !< A structure containing several
+                                                              !! relevant directory paths.
+  integer, target,               intent(inout) :: ntrunc      !< A target for the variable that
+                        !! records the number of times the velocity is truncated (this should be 0).
+
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (inout)   v - The meridional velocity, in m s-1.
 !  (inout)   h - The layer thicknesses, in m or kg m-2, depending on whether
diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90
index b9c9d30bd1..422655a89a 100644
--- a/src/core/MOM_dynamics_unsplit.F90
+++ b/src/core/MOM_dynamics_unsplit.F90
@@ -183,25 +183,48 @@ module MOM_dynamics_unsplit
 subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, &
                   p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, CS, &
                   VarMix, MEKE)
-  type(ocean_grid_type),                     intent(inout) :: 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_(G)), intent(inout) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  type(thermo_var_ptrs),                     intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
-  type(vertvisc_type),                       intent(inout) :: visc
-  type(time_type),                           intent(in)    :: Time_local
-  real,                                      intent(in)    :: dt   !< The dynamics time step, in s.
-  type(forcing),                             intent(in)    :: fluxes
-  real, dimension(:,:),                      pointer       :: p_surf_begin, p_surf_end
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uh
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vh
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr
-  real, dimension(SZI_(G),SZJ_(G)),          intent(out)   :: eta_av
-  type(MOM_dyn_unsplit_CS),                  pointer       :: CS
-  type(VarMix_CS),                           pointer       :: VarMix
-  type(MEKE_type),                           pointer       :: MEKE
+  type(ocean_grid_type),            intent(inout) :: 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_(G)), &
+                                    intent(inout) :: u      !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                    intent(inout) :: v      !< The meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                                    intent(inout) :: h      !< Layer thicknesses, in H.
+                                                            !! (usually m or kg m-2).
+  type(thermo_var_ptrs),            intent(in)    :: tv     !< A structure pointing to various
+                                                            !! thermodynamic variables.
+  type(vertvisc_type),              intent(inout) :: visc   !< A structure containing vertical
+                                 !! viscosities, bottom drag viscosities, and related fields.
+  type(time_type),                  intent(in)    :: Time_local   !< The model time at the end
+                                                                  !! of the time step.
+  real,                             intent(in)    :: dt     !< The dynamics time step, in s.
+  type(forcing),                    intent(in)    :: fluxes !< A structure containing pointers to
+                                 !! any possible forcing fields.  Unused fields have NULL ptrs.
+  real, dimension(:,:),             pointer       :: p_surf_begin !< A pointer (perhaps NULL) to the
+                                 !! surface pressure at the beginning of this dynamic step, in Pa.
+  real, dimension(:,:),             pointer       :: p_surf_end   !< A pointer (perhaps NULL) to the
+                                 !! surface pressure at the end of this dynamic step, in Pa.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                                    intent(inout) :: uh     !< The zonal volume or mass transport,
+                                                            !! in m3 s-1 or kg s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                    intent(inout) :: vh     !< The meridional volume or mass
+                                                            !! transport, in m3 s-1 or kg s-1.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                                    intent(inout) :: uhtr   !< he accumulated zonal volume or mass
+                                        !! transport since the last tracer advection, in m3 or kg.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                    intent(inout) :: vhtr   !< The accumulated meridional volume or
+                                  !! mass transport since the last tracer advection, in m3 or kg.
+  real, dimension(SZI_(G),SZJ_(G)), intent(out)   :: eta_av !< The time-mean free surface height or
+                                                            !! column mass, in m or kg m-2.
+  type(MOM_dyn_unsplit_CS),         pointer       :: CS     !< The control structure set up by
+                                                            !! initialize_dyn_unsplit.
+  type(VarMix_CS),                  pointer       :: VarMix !< A pointer to a structure with fields
+                                  !! that specify the spatially variable viscosities.
+  type(MEKE_type),                  pointer       :: MEKE   !< A pointer to a structure containing
+                                  !! fields related to the Mesoscale Eddy Kinetic Energy.
 ! Arguments: u - The input and output zonal velocity, in m s-1.
 !  (inout)   v - The input and output meridional velocity, in m s-1.
 !  (inout)   h - The input and output layer thicknesses, in m or kg m-2,
@@ -540,11 +563,14 @@ end subroutine step_MOM_dyn_unsplit
 ! =============================================================================
 
 subroutine register_restarts_dyn_unsplit(HI, GV, param_file, CS, restart_CS)
-  type(hor_index_type),         intent(in)    :: HI
-  type(verticalGrid_type),      intent(in)    :: GV   !< The ocean's vertical grid structure
-  type(param_file_type),        intent(in)    :: param_file !< A structure to parse for run-time parameters
-  type(MOM_dyn_unsplit_CS),     pointer       :: CS
-  type(MOM_restart_CS),         pointer       :: restart_CS
+  type(hor_index_type),      intent(in) :: HI         !< A horizontal index type structure.
+  type(verticalGrid_type),   intent(in) :: GV         !< The ocean's vertical grid structure.
+  type(param_file_type),     intent(in) :: param_file !< A structure to parse for
+                                                      !! run-time parameters.
+  type(MOM_dyn_unsplit_CS),  pointer    :: CS         !< The control structure set up by
+                                                      !! initialize_dyn_unsplit.
+  type(MOM_restart_CS),      pointer    :: restart_CS !< A pointer to the restart control structure.
+
 !   This subroutine sets up any auxiliary restart variables that are specific
 ! to the unsplit time stepping scheme.  All variables registered here should
 ! have the ability to be recreated if they are not present in a restart file.
@@ -589,26 +615,52 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, param_file, diag, CS, &
                                   restart_CS, Accel_diag, Cont_diag, MIS, &
                                   OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
                                   visc, dirs, ntrunc)
-  type(ocean_grid_type),                     intent(inout) :: 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_(G)), intent(inout) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  type(time_type),                   target, intent(in)    :: Time
-  type(param_file_type),                     intent(in)    :: param_file !< A structure to parse for run-time parameters
-  type(diag_ctrl),                   target, intent(inout) :: diag
-  type(MOM_dyn_unsplit_CS),                  pointer       :: CS
-  type(MOM_restart_CS),                      pointer       :: restart_CS
-  type(accel_diag_ptrs),             target, intent(inout) :: Accel_diag
-  type(cont_diag_ptrs),              target, intent(inout) :: Cont_diag
-  type(ocean_internal_state),                intent(inout) :: MIS
-  type(ocean_OBC_type),                      pointer       :: OBC
-  type(update_OBC_CS),                       pointer       :: update_OBC_CSp
-  type(ALE_CS),                              pointer       :: ALE_CSp
-  type(set_visc_CS),                         pointer       :: setVisc_CSp
-  type(vertvisc_type),                       intent(inout) :: visc
-  type(directories),                         intent(in)    :: dirs
-  integer, target,                           intent(inout) :: ntrunc
+  type(ocean_grid_type),          intent(inout) :: 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_(G)), &
+                                  intent(inout) :: u          !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                  intent(inout) :: v          !< The meridional velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , &
+                                  intent(inout) :: h          !< Layer thicknesses, in H
+                                                              !! (usually m or kg m-2).
+  type(time_type),        target, intent(in)    :: Time       !< The current model time.
+  type(param_file_type),          intent(in)    :: param_file !< A structure to parse
+                                                              !! for run-time parameters.
+  type(diag_ctrl),        target, intent(inout) :: diag       !< A structure that is used to
+                                                              !! regulate diagnostic output.
+  type(MOM_dyn_unsplit_CS),       pointer       :: CS         !< The control structure set up
+                                                              !! by initialize_dyn_unsplit.
+  type(MOM_restart_CS),           pointer       :: restart_CS !< A pointer to the restart control
+                                                              !!structure.
+  type(accel_diag_ptrs),  target, intent(inout) :: Accel_diag !< A set of pointers to the various
+                                     !! accelerations in the momentum equations, which can be used
+                                     !! for later derived diagnostics, like energy budgets.
+  type(cont_diag_ptrs),   target, intent(inout) :: Cont_diag  !< A structure with pointers to
+                                                              !! various terms in the continuity
+                                                              !! equations.
+  type(ocean_internal_state),     intent(inout) :: MIS        !< The "MOM6 Internal State"
+                                                   !! structure, used to pass around pointers
+                                                   !! to various arrays for diagnostic purposes.
+  type(ocean_OBC_type),           pointer       :: OBC        !< If open boundary conditions are
+                                                       !! used, this points to the ocean_OBC_type
+                                                       !! that was set up in MOM_initialization.
+  type(update_OBC_CS),            pointer       :: update_OBC_CSp !< If open boundary condition
+                                                            !! updates are used, this points to
+                                                            !! the appropriate control structure.
+  type(ALE_CS),                   pointer       :: ALE_CSp    !< This points to the ALE control
+                                                              !! structure.
+  type(set_visc_CS),              pointer       :: setVisc_CSp !< This points to the set_visc
+                                                               !! control structure.
+  type(vertvisc_type),            intent(inout) :: visc       !< A structure containing vertical
+                                                              !! viscosities, bottom drag
+                                                              !! viscosities, and related fields.
+  type(directories),              intent(in)    :: dirs       !< A structure containing several
+                                                              !! relevant directory paths.
+  integer, target,                intent(inout) :: ntrunc     !< A target for the variable that
+                                                        !! records the number of times the velocity
+                                                        !! is truncated (this should be 0).
 
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (inout)   v - The meridional velocity, in m s-1.
diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90
index 2fbd5356d4..0275c539d3 100644
--- a/src/core/MOM_dynamics_unsplit_RK2.F90
+++ b/src/core/MOM_dynamics_unsplit_RK2.F90
@@ -189,25 +189,61 @@ module MOM_dynamics_unsplit_RK2
 subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, fluxes, &
                   p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, CS, &
                   VarMix, MEKE)
-  type(ocean_grid_type),                     intent(inout) :: 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_(G)), intent(inout) :: u_in
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_in
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h_in
-  type(thermo_var_ptrs),                     intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
-  type(vertvisc_type),                       intent(inout) :: visc
-  type(time_type),                           intent(in)    :: Time_local
-  real,                                      intent(in)    :: dt   !< The baroclinic dynamics time step, in s.
-  type(forcing),                             intent(in)    :: fluxes
-  real, dimension(:,:),                      pointer       :: p_surf_begin, p_surf_end
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uh
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vh
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr
-  real, dimension(SZI_(G),SZJ_(G)),          intent(out)   :: eta_av
-  type(MOM_dyn_unsplit_RK2_CS),              pointer       :: CS
-  type(VarMix_CS),                           pointer       :: VarMix
-  type(MEKE_type),                           pointer       :: MEKE
+  type(ocean_grid_type),             intent(inout) :: 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_(G)), &
+                                     intent(inout) :: u_in    !< The input and output zonal
+                                                              !! velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                     intent(inout) :: v_in    !< The input and output meridional
+                                                              !! velocity, in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
+                                     intent(inout) :: h_in    !< The input and output layer
+                                              !! thicknesses, in m or kg m-2, depending on
+                                              !! whether the Boussinesq approximation is made.
+  type(thermo_var_ptrs),             intent(in)    :: tv      !< A structure pointing to various
+                                                              !! thermodynamic variables.
+  type(vertvisc_type),               intent(inout) :: visc    !< A structure containing vertical
+                                                              !! viscosities, bottom drag
+                                                              !! viscosities, and related fields.
+  type(time_type),                   intent(in)    :: Time_local   !< The model time at the end of
+                                                              !! the time step.
+  real,                              intent(in)    :: dt      !< The baroclinic dynamics time step,
+                                                              !! in s.
+  type(forcing),                     intent(in)    :: fluxes  !< A structure containing pointers to
+                                                              !! any possible forcing fields. Unused
+                                                              !! fields have NULL ptrs.
+  real, dimension(:,:),              pointer       :: p_surf_begin !< A pointer (perhaps NULL) to
+                                                            !! the surface pressure at the beginning
+                                                            !! of this dynamic step, in Pa.
+  real, dimension(:,:),              pointer       :: p_surf_end   !< A pointer (perhaps NULL) to
+                                                              !! the surface pressure at the end of
+                                                              !! this dynamic step, in Pa.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                                     intent(inout) :: uh      !< The zonal volume or mass transport,
+                                                              !! in m3 s-1 or kg s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                     intent(inout) :: vh      !< The meridional volume or mass
+                                                              !! transport, in m3 s-1 or kg s-1.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
+                                     intent(inout) :: uhtr    !< The accumulated zonal volume or
+                                                              !! mass transport since the last
+                                                              !! tracer advection, in m3 or kg.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
+                                     intent(inout) :: vhtr    !< The accumulated meridional volume
+                                                              !! or mass transport since the last
+                                                              !! tracer advection, in m3 or kg.
+  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: eta_av  !< The time-mean free surface height
+                                                              !! or column mass, in m or kg m-2.
+  type(MOM_dyn_unsplit_RK2_CS),      pointer       :: CS      !< The control structure set up by
+                                                              !! initialize_dyn_unsplit_RK2.
+  type(VarMix_CS),                   pointer       :: VarMix  !< A pointer to a structure with
+                                                              !! fields that specify the spatially
+                                                              !! variable viscosities.
+  type(MEKE_type),                   pointer       :: MEKE    !< A pointer to a structure containing
+                                                              !! fields related to the Mesoscale
+                                                              !! Eddy Kinetic Energy.
 ! Arguments: u_in - The input and output zonal velocity, in m s-1.
 !  (inout)   v_in - The input and output meridional velocity, in m s-1.
 !  (inout)   h_in - The input and output layer thicknesses, in m or kg m-2,
@@ -478,11 +514,14 @@ end subroutine step_MOM_dyn_unsplit_RK2
 ! =============================================================================
 
 subroutine register_restarts_dyn_unsplit_RK2(HI, GV, param_file, CS, restart_CS)
-  type(hor_index_type),         intent(in)    :: HI
-  type(verticalGrid_type),      intent(in)    :: GV   !< The ocean's vertical grid structure
-  type(param_file_type),        intent(in)    :: param_file !< A structure to parse for run-time parameters
-  type(MOM_dyn_unsplit_RK2_CS), pointer       :: CS
-  type(MOM_restart_CS),         pointer       :: restart_CS
+  type(hor_index_type),         intent(in)    :: HI         !< A horizontal index type structure.
+  type(verticalGrid_type),      intent(in)    :: GV         !< The ocean's vertical grid structure.
+  type(param_file_type),        intent(in)    :: param_file !< A structure to parse for run-time
+                                                            !! parameters.
+  type(MOM_dyn_unsplit_RK2_CS), pointer       :: CS         !< The control structure set up by
+                                                            !! initialize_dyn_unsplit_RK2.
+  type(MOM_restart_CS),         pointer       :: restart_CS !< A pointer to the restart control
+                                                            !! structure.
 !   This subroutine sets up any auxiliary restart variables that are specific
 ! to the unsplit time stepping scheme.  All variables registered here should
 ! have the ability to be recreated if they are not present in a restart file.
@@ -526,26 +565,51 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, param_file, diag, CS
                                       restart_CS, Accel_diag, Cont_diag, MIS, &
                                       OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
                                       visc, dirs, ntrunc)
-  type(ocean_grid_type),                     intent(inout) :: 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_(G)), intent(inout) :: u    !< The zonal velocity, in m s-1
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v    !< The meridional velocity, in m s-1
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  type(time_type),                   target, intent(in)    :: Time
-  type(param_file_type),                     intent(in)    :: param_file !< A structure to parse for run-time parameters
-  type(diag_ctrl),                   target, intent(inout) :: diag
-  type(MOM_dyn_unsplit_RK2_CS),              pointer       :: CS
-  type(MOM_restart_CS),                      pointer       :: restart_CS
-  type(accel_diag_ptrs),             target, intent(inout) :: Accel_diag
-  type(cont_diag_ptrs),              target, intent(inout) :: Cont_diag
-  type(ocean_internal_state),                intent(inout) :: MIS
-  type(ocean_OBC_type),                      pointer       :: OBC
-  type(update_OBC_CS),                       pointer       :: update_OBC_CSp
-  type(ALE_CS),                              pointer       :: ALE_CSp
-  type(set_visc_CS),                         pointer       :: setVisc_CSp
-  type(vertvisc_type),                       intent(inout) :: visc
-  type(directories),                         intent(in)    :: dirs
-  integer, target,                           intent(inout) :: ntrunc
+  type(ocean_grid_type),                     intent(inout) :: 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_(G)), intent(inout) :: u    !< The zonal velocity, in m s-1.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v    !< The meridional velocity,
+                                                                   !! in m s-1.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h    !< Layer thicknesses, in H
+                                                                   !! (usually m or kg m-2).
+  type(time_type),                   target, intent(in)    :: Time !< The current model time.
+  type(param_file_type),                     intent(in)    :: param_file !< A structure to parse
+                                                                         !! for run-time parameters.
+  type(diag_ctrl),                   target, intent(inout) :: diag !< A structure that is used to
+                                                                   !! regulate diagnostic output.
+  type(MOM_dyn_unsplit_RK2_CS),              pointer       :: CS   !< The control structure set up
+                                                                   !! by initialize_dyn_unsplit_RK2.
+  type(MOM_restart_CS),                      pointer       :: restart_CS !< A pointer to the restart
+                                                                         !! control structure.
+  type(accel_diag_ptrs),             target, intent(inout) :: Accel_diag !< A set of pointers to the
+                                      !! various accelerations in the momentum equations, which can
+                                      !! be used for later derived diagnostics, like energy budgets.
+  type(cont_diag_ptrs),              target, intent(inout) :: Cont_diag  !<A structure with pointers
+                                                                         !! to various terms in the
+                                                                         !! continuity equations.
+  type(ocean_internal_state),                intent(inout) :: MIS  !< The "MOM6 Internal State"
+                                                     !! structure, used to pass around pointers
+                                                     !! to various arrays for diagnostic purposes.
+  type(ocean_OBC_type),                      pointer       :: OBC  !< If open boundary conditions
+                                                    !! are used, this points to the ocean_OBC_type
+                                                    !! that was set up in MOM_initialization.
+  type(update_OBC_CS),                       pointer       :: update_OBC_CSp !< If open boundary
+                                                         !! condition updates are used, this points
+                                                         !! to the appropriate control structure.
+  type(ALE_CS),                              pointer       :: ALE_CSp     !< This points to the ALE
+                                                                          !! control structure.
+  type(set_visc_CS),                         pointer       :: setVisc_CSp !< This points to the
+                                                                          !! set_visc control
+                                                                          !! structure.
+  type(vertvisc_type),                       intent(inout) :: visc !< A structure containing
+                                                         !! vertical viscosities, bottom drag
+                                                         !! viscosities, and related fields.
+  type(directories),                         intent(in)    :: dirs !< A structure containing several
+                                                                   !! relevant directory paths.
+  integer, target,                           intent(inout) :: ntrunc !< A target for the variable
+                                                       !! that records the number of times the
+                                                       !! velocity is truncated (this should be 0).
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (inout)   v - The meridional velocity, in m s-1.
 !  (inout)   h - The layer thicknesses, in m or kg m-2, depending on whether
diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90
index d8efbe8561..aa50373755 100644
--- a/src/core/MOM_interface_heights.F90
+++ b/src/core/MOM_interface_heights.F90
@@ -67,14 +67,25 @@ module MOM_interface_heights
 contains
 
 subroutine find_eta_3d(h, tv, G_Earth, G, GV, eta, eta_bt, halo_size)
-  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(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  type(thermo_var_ptrs),                      intent(in)  :: tv   !< A structure pointing to various thermodynamic variables
-  real,                                       intent(in)  :: G_Earth
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: eta
-  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)  :: eta_bt
-  integer,                          optional, intent(in)  :: halo_size
+  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(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h         !< Layer thicknesses, in H
+                                                                       !! (usually m or kg m-2).
+  type(thermo_var_ptrs),                      intent(in)  :: tv        !< A structure pointing to
+                                                                       !! various thermodynamic
+                                                                       !! variables.
+  real,                                       intent(in)  :: G_Earth   !< Earth gravitational
+                                                                       !! acceleration (m/s2).
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: eta       !< layer interface heights
+                                                                       !! (meter).
+  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)  :: eta_bt    !< optional barotropic
+             !! variable that gives the "correct" free surface height (Boussinesq) or total water
+             !! column mass per unit aread (non-Boussinesq).  This is used to dilate the layer.
+             !! thicknesses when calculating interfaceheights, in H (m or kg m-2).
+  integer,                          optional, intent(in)  :: halo_size !< width of halo points on
+                                                                       !! which to calculate eta.
 
 !   This subroutine determines the heights of all interfaces between layers,
 ! using the appropriate form for consistency with the calculation of the
@@ -186,14 +197,24 @@ end subroutine find_eta_3d
 
 
 subroutine find_eta_2d(h, tv, G_Earth, G, GV, eta, eta_bt, halo_size)
-  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(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h    !< Layer thicknesses, in H (usually m or kg m-2)
-  type(thermo_var_ptrs),                      intent(in)  :: tv   !< A structure pointing to various thermodynamic variables
-  real,                                       intent(in)  :: G_Earth
-  real, dimension(SZI_(G),SZJ_(G)),           intent(out) :: eta
-  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)  :: eta_bt
-  integer,                          optional, intent(in)  :: halo_size
+  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(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h        !< Layer thicknesses, in H
+                                                                      !! (usually m or kg m-2).
+  type(thermo_var_ptrs),                      intent(in)  :: tv       !< A structure pointing to
+                                                                      !! various thermodynamic
+                                                                      !! variables.
+  real,                                       intent(in)  :: G_Earth  !< Earth gravitational
+                                                                      !! acceleration (m/s2).
+  real, dimension(SZI_(G),SZJ_(G)),           intent(out) :: eta      !< free surface height
+                                                                      !! relative to mean sea
+                                                                      !! level (z=0) (m).
+  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)  :: eta_bt   !< optional barotropic
+                   !! variable that gives the "correct" free surface height (Boussinesq) or total
+                   !! water column mass per unit aread (non-Boussinesq), in H (m or kg m-2).
+  integer,                          optional, intent(in)  :: halo_size !< width of halo points on
+                                                                       !! which to calculate eta.
 
 !   This subroutine determines the free surface height, using the appropriate
 ! form for consistency with the calculation of the pressure gradient forces.
diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90
index d2be2394d1..34965659f6 100644
--- a/src/core/MOM_verticalGrid.F90
+++ b/src/core/MOM_verticalGrid.F90
@@ -201,11 +201,21 @@ end function get_flux_units
 
 !> Returns the model's tracer flux units.
 function get_tr_flux_units(GV, tr_units, tr_vol_conc_units, tr_mass_conc_units)
-  character(len=48)                      :: get_tr_flux_units
-  type(verticalGrid_type),    intent(in) :: GV   !< The ocean's vertical grid structure
-  character(len=*), optional, intent(in) :: tr_units
-  character(len=*), optional, intent(in) :: tr_vol_conc_units
-  character(len=*), optional, intent(in) :: tr_mass_conc_units
+  character(len=48)                      :: get_tr_flux_units !< The model's flux units
+                                                              !! for a tracer.
+  type(verticalGrid_type),    intent(in) :: GV                !< The ocean's vertical
+                                                              !! grid structure.
+  character(len=*), optional, intent(in) :: tr_units          !< Units for a tracer, for example
+                                                              !! Celsius or PSU.
+  character(len=*), optional, intent(in) :: tr_vol_conc_units !< The concentration units per unit
+                                                              !! volume, forexample if the units are
+                                                              !! umol m-3, tr_vol_conc_units would
+                                                              !! be umol.
+  character(len=*), optional, intent(in) :: tr_mass_conc_units !< The concentration units per unit
+                                                              !! mass of sea water, for example if
+                                                              !! the units are mol kg-1,
+                                                              !! tr_vol_conc_units would be mol.
+
 !   This subroutine returns the appropriate units for thicknesses and fluxes,
 ! depending on whether the model is Boussinesq or not and the scaling for
 ! the vertical thickness.