diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90
index db97286f95..8a3f0ca29a 100644
--- a/src/ALE/MOM_regridding.F90
+++ b/src/ALE/MOM_regridding.F90
@@ -87,7 +87,7 @@ module MOM_regridding
   !> Reference pressure for potential density calculations (Pa)
   real :: ref_pressure = 2.e7
 
-  !> Weight given to old coordinate when blending between new and old grids (nondim)
+  !> Weight given to old coordinate when blending between new and old grids [nondim]
   !! Used only below depth_of_time_filter_shallow, with a cubic variation
   !! from zero to full effect between depth_of_time_filter_shallow and
   !! depth_of_time_filter_deep.
@@ -100,7 +100,7 @@ module MOM_regridding
   real :: depth_of_time_filter_deep = 0.
 
   !> Fraction (between 0 and 1) of compressibility to add to potential density
-  !! profiles when interpolating for target grid positions. (nondim)
+  !! profiles when interpolating for target grid positions. [nondim]
   real :: compressibility_fraction = 0.
 
   !> If true, each interface is given a maximum depth based on a rescaling of
@@ -184,6 +184,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
   character(len=*),           intent(in)    :: param_prefix !< String to prefix to parameter names.
                                                             !! If empty, causes main model parameters to be used.
   character(len=*),           intent(in)    :: param_suffix !< String to append to parameter names.
+
   ! Local variables
   integer :: ke ! Number of levels
   character(len=80)  :: string, string2, varName ! Temporary strings
@@ -194,18 +195,19 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
   logical :: tmpLogical, fix_haloclines, set_max, do_sum, main_parameters
   logical :: coord_is_state_dependent, ierr
   real :: filt_len, strat_tol, index_scale, tmpReal
-  real :: maximum_depth !< The maximum depth of the ocean [m] (not in Z).
+  real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
   real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
   real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
   integer :: nz_fixed_sfc, k, nzf(4)
-  real, dimension(:), allocatable :: dz     ! Resolution (thickness) in units of coordinate
+  real, dimension(:), allocatable :: dz     ! Resolution (thickness) in units of coordinate, which may be
+                                            ! [m] or [Z ~> m] or [H ~> m or kg m-2] or [kg m-3] or other units.
   real, dimension(:), allocatable :: h_max  ! Maximum layer thicknesses [H ~> m or kg m-2]
   real, dimension(:), allocatable :: z_max  ! Maximum interface depths [H ~> m or kg m-2] or other
                                             ! units depending on the coordinate
   real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths
                                             ! [H ~> m or kg m-2] or other units
-  real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode
-  ! Thicknesses that give level centers corresponding to table 2 of WOA09
+  real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3]
+  ! Thicknesses [m] that give level centers corresponding to table 2 of WOA09
   real, dimension(40) :: woa09_dz = (/ 5.,  10.,  10.,  15.,  22.5, 25., 25.,  25.,  &
                                       37.5, 50.,  50.,  75., 100., 100., 100., 100., &
                                      100., 100., 100., 100., 100., 100., 100., 175., &
@@ -918,10 +920,11 @@ end subroutine calc_h_new_by_dz
 
 !> Check that the total thickness of two grids match
 subroutine check_remapping_grid( G, GV, h, dzInterface, msg )
-  type(ocean_grid_type),                       intent(in) :: G !< Grid structure
-  type(verticalGrid_type),                     intent(in) :: GV !< Ocean vertical grid structure
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in) :: h !< Layer thicknesses (m)
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzInterface !< Change in interface positions (m)
+  type(ocean_grid_type),                       intent(in) :: G   !< Grid structure
+  type(verticalGrid_type),                     intent(in) :: GV  !< Ocean vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in) :: h   !< Layer thicknesses [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzInterface !< Change in interface positions
+                                                                 !! [H ~> m or kg m-2]
   character(len=*),                            intent(in) :: msg !< Message to append to errors
   ! Local variables
   integer :: i, j
@@ -936,8 +939,8 @@ end subroutine check_remapping_grid
 !> Check that the total thickness of new and old grids are consistent
 subroutine check_grid_column( nk, depth, h, dzInterface, msg )
   integer,               intent(in) :: nk !< Number of cells
-  real,                  intent(in) :: depth !< Depth of bottom (m or arbitrary units)
-  real, dimension(nk),   intent(in) :: h  !< Cell thicknesses (m or arbitrary units)
+  real,                  intent(in) :: depth !< Depth of bottom [Z ~> m] or arbitrary units
+  real, dimension(nk),   intent(in) :: h  !< Cell thicknesses [Z ~> m] or arbitrary units
   real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (same units as h)
   character(len=*),      intent(in) :: msg !< Message to append to errors
   ! Local variables
@@ -998,9 +1001,9 @@ end subroutine check_grid_column
 subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
   type(regridding_CS),      intent(in)    :: CS !< Regridding control structure
   integer,                  intent(in)    :: nk !< Number of cells in source grid
-  real, dimension(nk+1),    intent(in)    :: z_old !< Old grid position (m)
-  real, dimension(CS%nk+1), intent(in)    :: z_new !< New grid position (m)
-  real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions (m)
+  real, dimension(nk+1),    intent(in)    :: z_old !< Old grid position [m]
+  real, dimension(CS%nk+1), intent(in)    :: z_new !< New grid position [m]
+  real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [m]
   ! Local variables
   real :: sgn  ! The sign convention for downward.
   real :: dz_tgt, zr1, z_old_k
@@ -1475,7 +1478,7 @@ subroutine build_grid_HyCOM1( G, GV, h, tv, h_new, dzInterface, CS )
 
       z_col(1) = 0. ! Work downward rather than bottom up
       do K = 1, GV%ke
-        z_col(K+1) = z_col(K) + h(i,j,k) ! Work in units of [H ~> m or kg m-2]
+        z_col(K+1) = z_col(K) + h(i,j,k)
         p_col(k) = CS%ref_pressure + CS%compressibility_fraction * &
              ( 0.5 * ( z_col(K) + z_col(K+1) ) * GV%H_to_Pa - CS%ref_pressure )
       enddo
@@ -1607,7 +1610,7 @@ subroutine build_grid_SLight(G, GV, h, tv, dzInterface, CS)
       depth = G%bathyT(i,j) * GV%Z_to_H
       z_col(1) = 0. ! Work downward rather than bottom up
       do K=1,nz
-        z_col(K+1) = z_col(K) + h(i, j, k) ! Work in units of [H ~> m or kg m-2]
+        z_col(K+1) = z_col(K) + h(i,j,k)
         p_col(k) = CS%ref_pressure + CS%compressibility_fraction * &
                     ( 0.5 * ( z_col(K) + z_col(K+1) ) * GV%H_to_Pa - CS%ref_pressure )
       enddo
@@ -2221,19 +2224,19 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
   real,    optional, intent(in) :: Rho_ml_avg_depth !< Averaging depth over which to determine mixed layer potential
                                                     !! density [H ~> m or kg m-2]
   real,    optional, intent(in) :: nlay_ML_to_interior !< Number of layers to offset the mixed layer density to find
-                                                    !! resolved stratification (nondim)
+                                                    !! resolved stratification [nondim]
   logical, optional, intent(in) :: fix_haloclines   !< Detect regions with much weaker stratification in the coordinate
   real,    optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for
-                                                    !! spuriously unstable water mass profiles (m)
+                                                    !! spuriously unstable water mass profiles [m]
   real,    optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic
                                                     !! halocline region.
   logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward
                                                     !! from the top.
-  real,    optional, intent(in) :: adaptTimeRatio   !< Ratio of the ALE timestep to the grid timescale, ND.
+  real,    optional, intent(in) :: adaptTimeRatio   !< Ratio of the ALE timestep to the grid timescale [nondim].
   real,    optional, intent(in) :: adaptZoom        !< Depth of near-surface zooming region [H ~> m or kg m-2].
-  real,    optional, intent(in) :: adaptZoomCoeff   !< Coefficient of near-surface zooming diffusivity, ND.
-  real,    optional, intent(in) :: adaptBuoyCoeff   !< Coefficient of buoyancy diffusivity, ND.
-  real,    optional, intent(in) :: adaptAlpha       !< Scaling factor on optimization tendency, ND.
+  real,    optional, intent(in) :: adaptZoomCoeff   !< Coefficient of near-surface zooming diffusivity [nondim].
+  real,    optional, intent(in) :: adaptBuoyCoeff   !< Coefficient of buoyancy diffusivity [nondim].
+  real,    optional, intent(in) :: adaptAlpha       !< Scaling factor on optimization tendency [nondim].
   logical, optional, intent(in) :: adaptDoMin       !< If true, make a HyCOM-like mixed layer by
                                                     !! preventing interfaces from being shallower than
                                                     !! the depths specified by the regridding coordinate.
diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90
index c0620122c1..f399aa2c0f 100644
--- a/src/ALE/MOM_remapping.F90
+++ b/src/ALE/MOM_remapping.F90
@@ -71,7 +71,7 @@ module MOM_remapping
 ! outside of the range 0 to 1.
 #define __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
 
-real, parameter :: hNeglect_dflt = 1.E-30 !< A dimensional (H units) number that can be
+real, parameter :: hNeglect_dflt = 1.E-30 !< A thickness [H ~> m or kg m-2] that can be
                                       !! added to thicknesses in a denominator without
                                       !! changing the numerical result, except where
                                       !! a division by zero would otherwise occur.
diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90
index 4b38c683c7..6928425e33 100644
--- a/src/ALE/coord_hycom.F90
+++ b/src/ALE/coord_hycom.F90
@@ -39,8 +39,8 @@ module coord_hycom
 subroutine init_coord_hycom(CS, nk, coordinateResolution, target_density, interp_CS)
   type(hycom_CS),       pointer    :: CS !< Unassociated pointer to hold the control structure
   integer,              intent(in) :: nk !< Number of layers in generated grid
-  real, dimension(nk),  intent(in) :: coordinateResolution !< Nominal near-surface resolution (m)
-  real, dimension(nk+1),intent(in) :: target_density !< Interface target densities (kg/m3)
+  real, dimension(nk),  intent(in) :: coordinateResolution !< Nominal near-surface resolution [m]
+  real, dimension(nk+1),intent(in) :: target_density !< Interface target densities [kg m-3]
   type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation
 
   if (associated(CS)) call MOM_error(FATAL, "init_coord_hycom: CS already associated!")
@@ -100,13 +100,13 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
   type(EOS_type),        pointer       :: eqn_of_state !< Equation of state structure
   integer,               intent(in)    :: nz !< Number of levels
   real,                  intent(in)    :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2])
-  real, dimension(nz),   intent(in)    :: T !< Temperature of column (degC)
-  real, dimension(nz),   intent(in)    :: S !< Salinity of column (psu)
-  real, dimension(nz),   intent(in)    :: h  !< Layer thicknesses, (in m or H)
-  real, dimension(nz),   intent(in)    :: p_col !< Layer pressure in Pa
+  real, dimension(nz),   intent(in)    :: T !< Temperature of column [degC]
+  real, dimension(nz),   intent(in)    :: S !< Salinity of column [ppt]
+  real, dimension(nz),   intent(in)    :: h  !< Layer thicknesses, in [m] or [H ~> m or kg m-2]
+  real, dimension(nz),   intent(in)    :: p_col !< Layer pressure [Pa]
   real, dimension(nz+1), intent(in)    :: z_col !< Interface positions relative to the surface [H ~> m or kg m-2]
   real, dimension(CS%nk+1), intent(inout) :: z_col_new !< Absolute positions of interfaces
-  real, optional,        intent(in)    :: zScale !< Scaling factor from the input thicknesses in m
+  real, optional,        intent(in)    :: zScale !< Scaling factor from the input thicknesses in [m]
                                                  !! to desired units for zInterface, perhaps m_to_H.
   real,        optional, intent(in)    :: h_neglect !< A negligibly small width for the
                                              !! purpose of cell reconstructions
@@ -121,7 +121,7 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
   real, dimension(CS%nk) :: h_col_new ! New layer thicknesses
   real :: z_scale
   real :: stretching ! z* stretching, converts z* to z.
-  real :: nominal_z ! Nominal depth of interface is using z* (m or Pa)
+  real :: nominal_z ! Nominal depth of interface when using z* [Z ~> m]
   real :: hNew
   logical :: maximum_depths_set ! If true, the maximum depths of interface have been set.
   logical :: maximum_h_set      ! If true, the maximum layer thicknesses have been set.
diff --git a/src/ALE/coord_sigma.F90 b/src/ALE/coord_sigma.F90
index 0abc0c8e2b..3bf666ec52 100644
--- a/src/ALE/coord_sigma.F90
+++ b/src/ALE/coord_sigma.F90
@@ -28,7 +28,7 @@ module coord_sigma
 subroutine init_coord_sigma(CS, nk, coordinateResolution)
   type(sigma_CS),     pointer    :: CS !< Unassociated pointer to hold the control structure
   integer,            intent(in) :: nk !< Number of layers in the grid
-  real, dimension(:), intent(in) :: coordinateResolution !< Nominal coordinate resolution (nondim)
+  real, dimension(:), intent(in) :: coordinateResolution !< Nominal coordinate resolution [nondim]
 
   if (associated(CS)) call MOM_error(FATAL, "init_coord_sigma: CS already associated!")
   allocate(CS)
diff --git a/src/ALE/coord_slight.F90 b/src/ALE/coord_slight.F90
index d510388414..fa726548b8 100644
--- a/src/ALE/coord_slight.F90
+++ b/src/ALE/coord_slight.F90
@@ -24,14 +24,14 @@ module coord_slight
   real :: ref_pressure
 
   !> Fraction (between 0 and 1) of compressibility to add to potential density
-  !! profiles when interpolating for target grid positions. (nondim)
+  !! profiles when interpolating for target grid positions. [nondim]
   real :: compressibility_fraction
 
   ! The following 4 parameters were introduced for use with the SLight coordinate:
   !> Depth over which to average to determine the mixed layer potential density [H ~> m or kg m-2]
   real :: Rho_ML_avg_depth
 
-  !> Number of layers to offset the mixed layer density to find resolved stratification (nondim)
+  !> Number of layers to offset the mixed layer density to find resolved stratification [nondim]
   real :: nlay_ml_offset
 
   !> The number of fixed-thickness layers at the top of the model
@@ -48,7 +48,7 @@ module coord_slight
   !! unstable water mass profiles [H ~> m or kg m-2].
   real :: halocline_filter_length
 
-  !> A value of the stratification ratio that defines a problematic halocline region (nondim).
+  !> A value of the stratification ratio that defines a problematic halocline region [nondim].
   real :: halocline_strat_tol
 
   !> Nominal density of interfaces [kg m-3].
@@ -124,7 +124,7 @@ subroutine set_slight_params(CS, max_interface_depths, max_layer_thickness, &
                                       !! new grid through regridding [H ~> m or kg m-2]
   real,    optional, intent(in) :: compressibility_fraction !< Fraction (between 0 and 1) of
                                       !! compressibility to add to potential density profiles when
-                                      !! interpolating for target grid positions. (nondim)
+                                      !! interpolating for target grid positions. [nondim]
   real,    optional, intent(in) :: dz_ml_min        !< The fixed resolution in the topmost
                                       !! SLight_nkml_min layers [H ~> m or kg m-2]
   integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickness layers at the
@@ -132,13 +132,13 @@ subroutine set_slight_params(CS, max_interface_depths, max_layer_thickness, &
   real,    optional, intent(in) :: Rho_ML_avg_depth !< Depth over which to average to determine
                                       !! the mixed layer potential density [H ~> m or kg m-2]
   real,    optional, intent(in) :: nlay_ML_offset   !< Number of layers to offset the mixed layer
-                                      !! density to find resolved stratification (nondim)
+                                      !! density to find resolved stratification [nondim]
   logical, optional, intent(in) :: fix_haloclines   !< If true, detect regions with much weaker than
                                       !! based on in-situ density, and use a stretched coordinate there.
   real,    optional, intent(in) :: halocline_filter_length !< A length scale over which to filter T & S
                                       !! when looking for spuriously unstable water mass profiles [H ~> m or kg m-2].
   real,    optional, intent(in) :: halocline_strat_tol !< A value of the stratification ratio that
-                                      !! defines a problematic halocline region (nondim).
+                                      !! defines a problematic halocline region [nondim].
   type(interp_CS_type), &
            optional, intent(in) :: interp_CS !< Controls for interpolation
 
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index b5490ab497..1b364d7ef0 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -466,7 +466,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, &
   logical :: therm_reset ! If true, reset running sums of thermodynamic quantities.
   real :: cycle_time    ! The length of the coupled time-stepping cycle [s].
   real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
-    ssh         ! sea surface height, which may be based on eta_av (meter)
+    ssh         ! sea surface height, which may be based on eta_av [m]
 
   real, dimension(:,:,:), pointer :: &
     u => NULL(), & ! u : zonal velocity component [m s-1]
@@ -2488,7 +2488,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp)
   type(unit_scale_type),   pointer :: US => NULL() ! Pointer to a structure containing
                                                    ! various unit conversion factors
   type(MOM_restart_CS),    pointer :: restart_CSp_tmp => NULL()
-  real, allocatable :: z_interface(:,:,:) ! Interface heights (meter)
+  real, allocatable :: z_interface(:,:,:) ! Interface heights [m]
   type(vardesc) :: vd
 
   call cpu_clock_begin(id_clock_init)
diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90
index 97bc80ef7f..1bb1fe9c5a 100644
--- a/src/core/MOM_PressureForce_analytic_FV.F90
+++ b/src/core/MOM_PressureForce_analytic_FV.F90
@@ -175,7 +175,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
                     ! to reduce the impact of truncation errors.
   real :: rho_in_situ(SZI_(G)) ! The in situ density [kg m-3].
   real :: Pa_to_H   ! A factor to convert from Pa to the thicknesss units (H).
-!  real :: oneatm = 101325.0  ! 1 atm in Pa (kg/ms2)
+!  real :: oneatm = 101325.0  ! 1 atm in [Pa] = [kg m-1 s-2]
   real, parameter :: C1_6 = 1.0/6.0
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
   integer :: i, j, k
diff --git a/src/core/MOM_PressureForce_blocked_AFV.F90 b/src/core/MOM_PressureForce_blocked_AFV.F90
index d94adf7872..5c6cbd3dcb 100644
--- a/src/core/MOM_PressureForce_blocked_AFV.F90
+++ b/src/core/MOM_PressureForce_blocked_AFV.F90
@@ -172,7 +172,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
                     ! to reduce the impact of truncation errors.
   real :: rho_in_situ(SZI_(G)) ! The in situ density [kg m-3].
   real :: Pa_to_H   ! A factor to convert from Pa to the thicknesss units (H).
-!  real :: oneatm = 101325.0  ! 1 atm in Pa (kg/ms2)
+!  real :: oneatm = 101325.0  ! 1 atm in [Pa] = [kg m-1 s-2]
   real, parameter :: C1_6 = 1.0/6.0
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
   integer :: is_bk, ie_bk, js_bk, je_bk, Isq_bk, Ieq_bk, Jsq_bk, Jeq_bk
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index f6fc793946..4cc413926d 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -56,60 +56,60 @@ module MOM_forcing_type
 
   ! surface buoyancy force, used when temperature is not a state variable
   real, pointer, dimension(:,:) :: &
-    buoy          => NULL()  !< buoyancy flux (m^2/s^3)
+    buoy          => NULL()  !< buoyancy flux [m2 s-3]
 
-  ! radiative heat fluxes into the ocean (W/m^2)
+  ! radiative heat fluxes into the ocean [W m-2]
   real, pointer, dimension(:,:) :: &
-    sw         => NULL(), & !< shortwave (W/m^2)
-    sw_vis_dir => NULL(), & !< visible, direct shortwave (W/m^2)
-    sw_vis_dif => NULL(), & !< visible, diffuse shortwave (W/m^2)
-    sw_nir_dir => NULL(), & !< near-IR, direct shortwave (W/m^2)
-    sw_nir_dif => NULL(), & !< near-IR, diffuse shortwave (W/m^2)
-    lw         => NULL()    !< longwave (W/m^2) (typically negative)
-
-  ! turbulent heat fluxes into the ocean (W/m^2)
+    sw         => NULL(), & !< shortwave [W m-2]
+    sw_vis_dir => NULL(), & !< visible, direct shortwave [W m-2]
+    sw_vis_dif => NULL(), & !< visible, diffuse shortwave [W m-2]
+    sw_nir_dir => NULL(), & !< near-IR, direct shortwave [W m-2]
+    sw_nir_dif => NULL(), & !< near-IR, diffuse shortwave [W m-2]
+    lw         => NULL()    !< longwave [W m-2] (typically negative)
+
+  ! turbulent heat fluxes into the ocean [W m-2]
   real, pointer, dimension(:,:) :: &
-    latent     => NULL(), & !< latent (W/m^2)   (typically < 0)
-    sens       => NULL(), & !< sensible (W/m^2) (typically negative)
-    heat_added => NULL()    !< additional heat flux from SST restoring or flux adjustments (W/m^2)
+    latent     => NULL(), & !< latent [W m-2]   (typically < 0)
+    sens       => NULL(), & !< sensible [W m-2] (typically negative)
+    heat_added => NULL()    !< additional heat flux from SST restoring or flux adjustments [W m-2]
 
   ! components of latent heat fluxes used for diagnostic purposes
   real, pointer, dimension(:,:) :: &
-    latent_evap_diag    => NULL(), & !< latent (W/m^2) from evaporating liquid water (typically < 0)
-    latent_fprec_diag   => NULL(), & !< latent (W/m^2) from melting fprec  (typically < 0)
-    latent_frunoff_diag => NULL()    !< latent (W/m^2) from melting frunoff (calving) (typically < 0)
+    latent_evap_diag    => NULL(), & !< latent [W m-2] from evaporating liquid water (typically < 0)
+    latent_fprec_diag   => NULL(), & !< latent [W m-2] from melting fprec  (typically < 0)
+    latent_frunoff_diag => NULL()    !< latent [W m-2] from melting frunoff (calving) (typically < 0)
 
-  ! water mass fluxes into the ocean ( kg/(m^2 s) ); these fluxes impact the ocean mass
+  ! water mass fluxes into the ocean [kg m-2 s-1]; these fluxes impact the ocean mass
   real, pointer, dimension(:,:) :: &
-    evap        => NULL(), & !< (-1)*fresh water flux evaporated out of the ocean ( kg/(m^2 s) )
-    lprec       => NULL(), & !< precipitating liquid water into the ocean ( kg/(m^2 s) )
-    fprec       => NULL(), & !< precipitating frozen water into the ocean ( kg/(m^2 s) )
-    vprec       => NULL(), & !< virtual liquid precip associated w/ SSS restoring ( kg/(m^2 s) )
-    lrunoff     => NULL(), & !< liquid river runoff entering ocean ( kg/(m^2 s) )
-    frunoff     => NULL(), & !< frozen river runoff (calving) entering ocean ( kg/(m^2 s) )
-    seaice_melt => NULL(), & !< seaice melt (positive) or formation (negative) ( kg/(m^2 s) )
-    netMassIn   => NULL(), & !< Sum of water mass flux out of the ocean ( kg/(m^2 s) )
-    netMassOut  => NULL(), & !< Net water mass flux into of the ocean ( kg/(m^2 s) )
-    netSalt     => NULL()    !< Net salt entering the ocean
+    evap        => NULL(), & !< (-1)*fresh water flux evaporated out of the ocean [kg m-2 s-1]
+    lprec       => NULL(), & !< precipitating liquid water into the ocean [kg m-2 s-1]
+    fprec       => NULL(), & !< precipitating frozen water into the ocean [kg m-2 s-1]
+    vprec       => NULL(), & !< virtual liquid precip associated w/ SSS restoring [kg m-2 s-1]
+    lrunoff     => NULL(), & !< liquid river runoff entering ocean [kg m-2 s-1]
+    frunoff     => NULL(), & !< frozen river runoff (calving) entering ocean [kg m-2 s-1]
+    seaice_melt => NULL(), & !< seaice melt (positive) or formation (negative) [kg m-2 s-1]
+    netMassIn   => NULL(), & !< Sum of water mass flux out of the ocean [kg m-2 s-1]
+    netMassOut  => NULL(), & !< Net water mass flux into of the ocean [kg m-2 s-1]
+    netSalt     => NULL()    !< Net salt entering the ocean [kgSalt m-2 s-1]
 
   ! heat associated with water crossing ocean surface
   real, pointer, dimension(:,:) :: &
-    heat_content_cond    => NULL(), & !< heat content associated with condensating water (W/m^2)
-    heat_content_lprec   => NULL(), & !< heat content associated with liquid >0 precip   (W/m^2) (diagnostic)
-    heat_content_fprec   => NULL(), & !< heat content associated with frozen precip      (W/m^2)
-    heat_content_vprec   => NULL(), & !< heat content associated with virtual >0 precip  (W/m^2)
-    heat_content_lrunoff => NULL(), & !< heat content associated with liquid runoff      (W/m^2)
-    heat_content_frunoff => NULL(), & !< heat content associated with frozen runoff      (W/m^2)
-    heat_content_icemelt => NULL(), & !< heat content associated with liquid sea ice     (W/m^2)
-    heat_content_massout => NULL(), & !< heat content associated with mass leaving ocean (W/m^2)
-    heat_content_massin  => NULL()    !< heat content associated with mass entering ocean (W/m^2)
+    heat_content_cond    => NULL(), & !< heat content associated with condensating water [W m-2]
+    heat_content_lprec   => NULL(), & !< heat content associated with liquid >0 precip   [W m-2] (diagnostic)
+    heat_content_fprec   => NULL(), & !< heat content associated with frozen precip      [W m-2]
+    heat_content_vprec   => NULL(), & !< heat content associated with virtual >0 precip  [W m-2]
+    heat_content_lrunoff => NULL(), & !< heat content associated with liquid runoff      [W m-2]
+    heat_content_frunoff => NULL(), & !< heat content associated with frozen runoff      [W m-2]
+    heat_content_icemelt => NULL(), & !< heat content associated with liquid sea ice     [W m-2]
+    heat_content_massout => NULL(), & !< heat content associated with mass leaving ocean [W m-2]
+    heat_content_massin  => NULL()    !< heat content associated with mass entering ocean [W m-2]
 
   ! salt mass flux (contributes to ocean mass only if non-Bouss )
   real, pointer, dimension(:,:) :: &
-    salt_flux       => NULL(), & !< net salt flux into the ocean ( kg salt/(m^2 s) )
-    salt_flux_in    => NULL(), & !< salt flux provided to the ocean from coupler ( kg salt/(m^2 s) )
+    salt_flux       => NULL(), & !< net salt flux into the ocean [kgSalt m-2 s-1]
+    salt_flux_in    => NULL(), & !< salt flux provided to the ocean from coupler [kgSalt m-2 s-1]
     salt_flux_added => NULL()    !< additional salt flux from restoring or flux adjustment before adjustment
-                                 !! to net zero ( kg salt/(m^2 s) )
+                                 !! to net zero [kgSalt m-2 s-1]
 
   ! applied surface pressure from other component models (e.g., atmos, sea ice, land ice)
   real, pointer, dimension(:,:) :: p_surf_full => NULL()
@@ -119,7 +119,7 @@ module MOM_forcing_type
                 !< Pressure at the top ocean interface [Pa] as used to drive the ocean model.
                 !! If p_surf is limited, p_surf may be smaller than p_surf_full, otherwise they are the same.
   real, pointer, dimension(:,:) :: p_surf_SSH => NULL()
-                !< Pressure at the top ocean interface that is used in corrections to the sea surface
+                !< Pressure at the top ocean interface [Pa] that is used in corrections to the sea surface
                 !! height field that is passed back to the calling routines.
                 !! p_surf_SSH may point to p_surf or to p_surf_full.
   logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere
@@ -129,14 +129,14 @@ module MOM_forcing_type
 
   ! tide related inputs
   real, pointer, dimension(:,:) :: &
-    TKE_tidal     => NULL(), & !< tidal energy source driving mixing in bottom boundary layer (W/m^2)
+    TKE_tidal     => NULL(), & !< tidal energy source driving mixing in bottom boundary layer [W m-2]
     ustar_tidal   => NULL()    !< tidal contribution to bottom ustar [m s-1]
 
   ! iceberg related inputs
   real, pointer, dimension(:,:) :: &
     ustar_berg => NULL(), &   !< iceberg contribution to top ustar [Z s-1 ~> m s-1].
-    area_berg  => NULL(), &   !< area of ocean surface covered by icebergs (m2/m2)
-    mass_berg  => NULL()      !< mass of icebergs (kg/m2)
+    area_berg  => NULL(), &   !< area of ocean surface covered by icebergs [m2 m-2]
+    mass_berg  => NULL()      !< mass of icebergs [kg m-2]
 
   ! land ice-shelf related inputs
   real, pointer, dimension(:,:) :: ustar_shelf => NULL()  !< Friction velocity under ice-shelves [Z s-1 ~> m s-1].
@@ -146,15 +146,15 @@ module MOM_forcing_type
                                  !! associated if ice shelves are enabled, and are
                                  !! exactly 0 away from shelves or on land.
   real, pointer, dimension(:,:) :: iceshelf_melt => NULL() !< Ice shelf melt rate (positive)
-                                 !! or freezing (negative) (in m/year)
+                                 !! or freezing (negative) [m year-1]
 
   ! Scalars set by surface forcing modules
-  real :: vPrecGlobalAdj     !< adjustment to restoring vprec to zero out global net ( kg/(m^2 s) )
-  real :: saltFluxGlobalAdj  !< adjustment to restoring salt flux to zero out global net ( kg salt/(m^2 s) )
-  real :: netFWGlobalAdj     !< adjustment to net fresh water to zero out global net ( kg/(m^2 s) )
-  real :: vPrecGlobalScl     !< scaling of restoring vprec to zero out global net ( -1..1 )
-  real :: saltFluxGlobalScl  !< scaling of restoring salt flux to zero out global net ( -1..1 )
-  real :: netFWGlobalScl     !< scaling of net fresh water to zero out global net ( -1..1 )
+  real :: vPrecGlobalAdj     !< adjustment to restoring vprec to zero out global net [kg m-2 s-1]
+  real :: saltFluxGlobalAdj  !< adjustment to restoring salt flux to zero out global net [kgSalt m-2 s-1]
+  real :: netFWGlobalAdj     !< adjustment to net fresh water to zero out global net [kg m-2 s-1]
+  real :: vPrecGlobalScl     !< scaling of restoring vprec to zero out global net ( -1..1 ) [nondim]
+  real :: saltFluxGlobalScl  !< scaling of restoring salt flux to zero out global net ( -1..1 ) [nondim]
+  real :: netFWGlobalScl     !< scaling of net fresh water to zero out global net ( -1..1 ) [nondim]
 
   logical :: fluxes_used = .true. !< If true, all of the heat, salt, and mass
                                   !! fluxes have been applied to the ocean.
@@ -162,7 +162,7 @@ module MOM_forcing_type
                                   !! should be applied [s].  If negative, this forcing
                                   !! type variable has not yet been inialized.
 
-  real :: C_p                !< heat capacity of seawater ( J/(K kg) ).
+  real :: C_p                !< heat capacity of seawater [J kg-1 degC-1].
                              !! C_p is is the same value as in thermovar_ptrs_type.
 
   ! passive tracer surface fluxes
@@ -884,9 +884,9 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, h, Temp, Salt, tv,
 
   ! The surface forcing is contained in the fluxes type.
   ! We aggregate the thermodynamic forcing for a time step into the following:
-  ! netH       = water (H units/s) added/removed via surface fluxes
-  ! netHeat    = heat (degC * H/s) via surface fluxes
-  ! netSalt    = salt ( g(salt)/m2 for non-Bouss and ppt*m for Bouss /s) via surface fluxes
+  ! netH       = water added/removed via surface fluxes [H s-1 ~> m s-1 or kg m-2 s-1]
+  ! netHeat    = heat via surface fluxes [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
+  ! netSalt    = salt via surface fluxes [ppt H s-1 ~> ppt m s-1 or gSalt m-2 s-1]
   ! Note that unlike other calls to extractFLuxes1d() that return the time-integrated flux
   ! this call returns the rate because dt=1
   call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,                                 &
diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90
index e5e2045c7a..f5775fc2c4 100644
--- a/src/core/MOM_isopycnal_slopes.F90
+++ b/src/core/MOM_isopycnal_slopes.F90
@@ -35,14 +35,14 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
                                                                      !! thermodynamic variables
   real,                                        intent(in)    :: dt_kappa_smooth !< A smoothing vertical diffusivity
                                                                      !! times a smoothing timescale [Z2 ~> m2].
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: slope_x !< Isopycnal slope in i-direction (nondim)
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: slope_y !< Isopycnal slope in j-direction (nondim)
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: slope_x !< Isopycnal slope in i-direction [nondim]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: slope_y !< Isopycnal slope in j-direction [nondim]
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), &
                                      optional, intent(inout) :: N2_u !< Brunt-Vaisala frequency squared at
-                                                                     !! interfaces between u-points (s-2)
+                                                                     !! interfaces between u-points [s-2]
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), &
                                      optional, intent(inout) :: N2_v !< Brunt-Vaisala frequency squared at
-                                                                     !! interfaces between u-points (s-2)
+                                                                     !! interfaces between u-points [s-2]
   integer,                           optional, intent(in)    :: halo !< Halo width over which to compute
 
   ! real,                              optional, intent(in)    :: eta_to_m !< The conversion factor from the units
@@ -234,7 +234,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
           slope_x(I,j,K) = 0.0
         endif
 
-        if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of Brunt-Vaisala frequency (s-2)
+        if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of Brunt-Vaisala frequency [s-2]
 
       else ! With .not.use_EOS, the layers are constant density.
         slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j)
@@ -318,7 +318,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
           slope_y(i,J,K) = 0.0
         endif
 
-        if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of Brunt-Vaisala frequency (s-2)
+        if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of Brunt-Vaisala frequency [s-2]
 
       else ! With .not.use_EOS, the layers are constant density.
         slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J)
diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90
index c78c08933e..c59eafc4c2 100644
--- a/src/core/MOM_open_boundary.F90
+++ b/src/core/MOM_open_boundary.F90
@@ -151,7 +151,7 @@ module MOM_open_boundary
   real, pointer, dimension(:,:,:) :: tangential_grad=>NULL() !< The gradient of the velocity tangential
                                                             !! to the OB segment [m s-1].
   real, pointer, dimension(:,:,:) :: normal_trans=>NULL()   !< The layer transport normal to the OB
-                                                            !! segment (m3 s-1).
+                                                            !! segment [m3 s-1].
   real, pointer, dimension(:,:)   :: normal_vel_bt=>NULL()  !< The barotropic velocity normal to
                                                             !! the OB segment [m s-1].
   real, pointer, dimension(:,:)   :: eta=>NULL()            !< The sea-surface elevation along the segment [m].
@@ -1412,8 +1412,8 @@ end subroutine open_boundary_impose_normal_slope
 subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv)
   type(ocean_OBC_type),              pointer       :: OBC !< Open boundary control structure
   type(dyn_horgrid_type),            intent(inout) :: G !< Ocean grid structure
-  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: areaCu !< Area of a u-cell (m2)
-  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: areaCv !< Area of a u-cell (m2)
+  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: areaCu !< Area of a u-cell [m2]
+  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: areaCv !< Area of a u-cell [m2]
   ! Local variables
   integer :: i, j, n
   type(OBC_segment_type), pointer :: segment => NULL()
@@ -2921,7 +2921,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time)
   type(verticalGrid_type),                   intent(in)    :: GV   !<  Ocean vertical grid structure
   type(ocean_OBC_type),                      pointer       :: OBC  !< Open boundary structure
   type(thermo_var_ptrs),                     intent(in)    :: tv   !< Thermodynamics structure
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h    !< Thickness
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h    !< Thickness [m]
   type(time_type),                           intent(in)    :: Time !< Model time
   ! Local variables
   integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed
@@ -2930,7 +2930,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time)
   character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path
   type(OBC_segment_type), pointer :: segment => NULL()
   integer, dimension(4) :: siz,siz2
-  real :: sumh ! column sum of thicknesses (m)
+  real :: sumh ! column sum of thicknesses [m]
   integer :: ni_seg, nj_seg  ! number of src gridpoints along the segments
   integer :: i2, j2          ! indices for referencing local domain array
   integer :: is_obc, ie_obc, js_obc, je_obc  ! segment indices within local domain
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 92b5f4f918..c623848c15 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -45,8 +45,8 @@ module MOM_variables
     v, &        !< The mixed layer meridional velocity [m s-1].
     sea_lev, &  !< The sea level [m].  If a reduced surface gravity is
                 !! used, that is compensated for in sea_lev.
-    melt_potential, & !< instantaneous amount of heat that can be used to melt sea ice,
-                      !! in J m-2. This is computed w.r.t. surface freezing temperature.
+    melt_potential, & !< Instantaneous amount of heat that can be used to melt sea ice [J m-2].
+                      !! This is computed w.r.t. surface freezing temperature.
     ocean_mass, &  !< The total mass of the ocean [kg m-2].
     ocean_heat, &  !< The total heat content of the ocean in [degC kg m-2].
     ocean_salt, &  !< The total salt content of the ocean in [kgSalt m-2].
@@ -219,7 +219,7 @@ module MOM_variables
   real, pointer, dimension(:,:) :: kv_tbl_shelf_v => NULL()
                 !< Viscosity in the viscous top boundary layer under ice shelves at v-points [Z2 s-1 ~> m2 s-1].
   real, pointer, dimension(:,:) :: nkml_visc_u => NULL()
-                !< The number of layers in the viscous surface mixed layer at u-points (nondimensional).
+                !< The number of layers in the viscous surface mixed layer at u-points [nondim].
                 !! This is not an integer because there may be fractional layers, and it is stored in
                 !! terms of layers, not depth, to facilitate the movement of the viscous boundary layer
                 !! with the flow.
diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90
index 07769f6077..3c50f00061 100644
--- a/src/diagnostics/MOM_diag_to_Z.F90
+++ b/src/diagnostics/MOM_diag_to_Z.F90
@@ -177,11 +177,11 @@ subroutine calculate_Z_diag_fields(u, v, h, ssh_in, frac_shelf_h, G, GV, US, CS)
   real :: htot           ! summed layer thicknesses [H ~> m or kg m-2]
   real :: dilate         ! proportion by which to dilate every layer
   real :: wt(SZK_(G)+1)  ! fractional weight for each layer in the
-                         ! range between k_top and k_bot (nondim)
+                         ! range between k_top and k_bot [nondim]
   real :: z1(SZK_(G)+1)  ! z1 and z2 are the depths of the top and bottom
   real :: z2(SZK_(G)+1)  ! limits of the part of a layer that contributes
                          ! to a depth level, relative to the cell center
-                         ! and normalized by the cell thickness (nondim)
+                         ! and normalized by the cell thickness [nondim]
                          ! Note that -1/2 <= z1 < z2 <= 1/2.
   real :: sl_tr(max(CS%num_tr_used,1)) ! normalized slope of the tracer
                                        ! within the cell, in tracer units
@@ -674,10 +674,10 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z
   real, dimension(:), intent(out)   :: wt     !< Relative weights of each layer from k_top to k_bot.
   real, dimension(:), intent(out)   :: z1     !< Depth of the top limits of the part of
        !! a layer that contributes to a depth level, relative to the cell center and normalized
-       !! by the cell thickness (nondim).  Note that -1/2 <= z1 < z2 <= 1/2.
+       !! by the cell thickness [nondim].  Note that -1/2 <= z1 < z2 <= 1/2.
   real, dimension(:), intent(out)   :: z2     !< Depths of the bottom limit of the part of
        !! a layer that contributes to a depth level, relative to the cell center and normalized
-       !! by the cell thickness (nondim).  Note that -1/2 <= z1 < z2 <= 1/2.
+       !! by the cell thickness [nondim].  Note that -1/2 <= z1 < z2 <= 1/2.
   ! Local variables
   real    :: Ih, e_c, tot_wt, I_totwt
   integer :: k
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index caa0f3cd00..cd3c87b922 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -233,13 +233,13 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
   real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS
   real :: wt, wt_p
 
-  ! squared Coriolis parameter at to h-points (1/s2)
+  ! squared Coriolis parameter at to h-points [s-2]
   real :: f2_h
 
-  ! magnitude of the gradient of f (1/(m*s))
+  ! magnitude of the gradient of f [s-1 m-1]
   real :: mag_beta
 
-  ! frequency squared used to avoid division by 0 (1/s2)
+  ! frequency squared used to avoid division by 0 [s-2]
   ! value is roughly (pi / (the age of the universe) )^2.
   real, parameter :: absurdly_small_freq2 = 1e-34
 
@@ -327,7 +327,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     call post_data(CS%id_masso, masso, CS%diag)
   endif
 
-  ! diagnose thickness/volumes of grid cells (meter)
+  ! diagnose thickness/volumes of grid cells [m]
   if (CS%id_thkcello>0 .or. CS%id_volcello>0) then
     if (GV%Boussinesq) then ! thkcello = h for Boussinesq
       if (CS%id_thkcello > 0) then ; if (GV%H_to_m == 1.0) then
@@ -359,7 +359,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
           do i=is,ie ! Pressure for EOS at the layer center [Pa]
             pressure_1d(i) = pressure_1d(i) + 0.5*GV%H_to_Pa*h(i,j,k)
           enddo
-          ! Store in-situ density (kg/m3) in work_3d
+          ! Store in-situ density [kg m-3] in work_3d
           call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, &
                                  work_3d(:,j,k), is, ie-is+1, tv%eqn_of_state)
           do i=is,ie ! Cell thickness = dz = dp/(g*rho) (meter); store in work_3d
@@ -1140,7 +1140,7 @@ subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh)
   type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
   type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
   real, dimension(SZI_(G),SZJ_(G)), &
-                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement (m)
+                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
 
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G)) :: work_2d  ! A 2-d work array
@@ -1181,13 +1181,13 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
   type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
   type(thermo_var_ptrs),    intent(in) :: tv  !< A structure pointing to various thermodynamic variables
   real, dimension(SZI_(G),SZJ_(G)), &
-                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement (m)
+                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
   real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections
-                                                  !! for ice displacement and the inverse barometer (m)
+                                                  !! for ice displacement and the inverse barometer [m]
 
   real, dimension(SZI_(G),SZJ_(G)) :: work_2d  ! A 2-d work array
   real, dimension(SZI_(G),SZJ_(G)) :: &
-    zos  ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh (meter)
+    zos  ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [m]
   real :: I_time_int    ! The inverse of the time interval [s-1].
   real :: zos_area_mean, volo, ssh_ga
   integer :: i, j, is, ie, js, je
diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90
index 872c6f4783..0c4b0386a4 100644
--- a/src/diagnostics/MOM_wave_speed.F90
+++ b/src/diagnostics/MOM_wave_speed.F90
@@ -67,7 +67,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
                                           !! monotonic for the purposes of calculating vertical
                                           !! modal structure [m].
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
-        optional,                   intent(out) :: modal_structure !< Normalized model structure (non-dim)
+        optional,                   intent(out) :: modal_structure !< Normalized model structure [nondim]
 
   ! Local variables
   real, dimension(SZK_(G)+1) :: &
diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90
index 64592e7303..28ad4c6bfc 100644
--- a/src/diagnostics/MOM_wave_structure.F90
+++ b/src/diagnostics/MOM_wave_structure.F90
@@ -99,15 +99,12 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
                                                               !! gravity wave speed [m s-1].
   integer,                                  intent(in)  :: ModeNum !< Mode number
   real,                                     intent(in)  :: freq !< Intrinsic wave frequency [s-1].
-  type(wave_structure_CS),                  pointer     :: CS !< The control structure returned
-                                                              !! by a previous call to
-                                                              !! wave_structure_init.
+  type(wave_structure_CS),                  pointer     :: CS !< The control structure returned by a
+                                                              !! previous call to wave_structure_init.
   real, dimension(SZI_(G),SZJ_(G)), &
-                                  optional, intent(in)  :: En !< Internal wave energy density,
-                                                              !! in Jm-2.
+                                  optional, intent(in)  :: En !< Internal wave energy density [J m-2].
   logical,optional,                         intent(in)  :: full_halos !< If true, do the calculation
-                                                              !! over the entire computational
-                                                              !! domain.
+                                                              !! over the entire computational domain.
   ! Local variables
   real, dimension(SZK_(G)+1) :: &
     dRho_dT, dRho_dS, &
@@ -115,8 +112,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
     gprime        ! The reduced gravity across each interface [m2 Z-1 s-2 ~> m s-2].
   real, dimension(SZK_(G)) :: &
     Igl, Igu      ! The inverse of the reduced gravity across an interface times
-                  ! the thickness of the layer below (Igl) or above (Igu) it,
-                  ! [s2 m-2].
+                  ! the thickness of the layer below (Igl) or above (Igu) it [s2 m-2].
   real, dimension(SZK_(G),SZI_(G)) :: &
     Hf, Tf, Sf, Rf
   real, dimension(SZK_(G)) :: &
diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90
index e4435a60df..4538f9c8ae 100644
--- a/src/equation_of_state/MOM_EOS.F90
+++ b/src/equation_of_state/MOM_EOS.F90
@@ -200,7 +200,7 @@ subroutine calculate_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref)
   real,           intent(in)  :: T !< Potential temperature referenced to the surface [degC]
   real,           intent(in)  :: S !< Salinity [ppt]
   real,           intent(in)  :: pressure !< Pressure [Pa]
-  real,           intent(out) :: specvol  !< specific volume (in-situ if pressure is local) (m3 kg-1)
+  real,           intent(out) :: specvol  !< specific volume (in-situ if pressure is local) [m3 kg-1]
   type(EOS_type), pointer     :: EOS      !< Equation of state structure
   real, optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
 
@@ -1379,7 +1379,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
   real,           intent(in)  :: P_t !< Anomalous pressure of top of cell, relative to g*rho_ref*z_t [Pa]
   real,           intent(in)  :: P_tgt !< Target pressure at height z_out, relative to g*rho_ref*z_out [Pa]
   real,           intent(in)  :: rho_ref !< Reference density with which calculation are anomalous to
-  real,           intent(in)  :: G_e !< Gravitational acceleration (m2 Z-1 s-2)
+  real,           intent(in)  :: G_e !< Gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
   type(EOS_type), pointer     :: EOS !< Equation of state structure
   real,           intent(out) :: P_b !< Pressure at the bottom of the cell [Pa]
   real,           intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m].
diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90
index db0afa3d8a..c516c96e86 100644
--- a/src/framework/MOM_io.F90
+++ b/src/framework/MOM_io.F90
@@ -97,7 +97,7 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit
   integer,               intent(in)    :: novars     !< number of fields written to filename
   type(fieldtype),       intent(inout) :: fields(:)  !< array of fieldtypes for each variable
   integer, optional,     intent(in)    :: threading  !< SINGLE_FILE or MULTIPLE
-  real, optional,        intent(in)    :: timeunit   !< length, in seconds, of the units for time. The
+  real, optional,        intent(in)    :: timeunit   !< length of the units for time [s]. The
                                                      !! default value is 86400.0, for 1 day.
   type(ocean_grid_type),   optional, intent(in) :: G !< ocean horizontal grid structure; G or dG
                                                      !! is required if the new file uses any
@@ -357,7 +357,7 @@ subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit
   integer,               intent(in)    :: novars     !< number of fields written to filename
   type(fieldtype),       intent(inout) :: fields(:)  !< array of fieldtypes for each variable
   integer, optional,     intent(in)    :: threading  !< SINGLE_FILE or MULTIPLE
-  real, optional,        intent(in)    :: timeunit   !< length, in seconds, of the units for time. The
+  real, optional,        intent(in)    :: timeunit   !< length of the units for time [s]. The
                                                      !! default value is 86400.0, for 1 day.
   type(ocean_grid_type),   optional, intent(in) :: G !< ocean horizontal grid structure; G or dG
                                                      !! is required if a new file uses any
diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90
index f3f8e761b7..5cbf1f4056 100644
--- a/src/ice_shelf/MOM_ice_shelf.F90
+++ b/src/ice_shelf/MOM_ice_shelf.F90
@@ -629,7 +629,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!End of safety checks !!!!!!!!!!!!!!!!!!!
   enddo ; enddo ! i- and j-loops
 
-  ! mass flux (kg/s), part of ISOMIP diags.
+  ! mass flux [kg s-1], part of ISOMIP diags.
   mass_flux(:,:) = 0.0
   mass_flux(:,:) = ISS%water_flux(:,:) * ISS%area_shelf_h(:,:)
 
@@ -1056,7 +1056,7 @@ subroutine add_shelf_flux(G, CS, state, fluxes)
     enddo ; enddo
 
     if (CS%DEBUG) then
-      write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux, CS%time_step
+      write(mesg,*) 'Mean melt flux [kg m-2 s-1], dt = ', mean_melt_flux, CS%time_step
       call MOM_mesg(mesg)
       call MOM_forcing_chksum("After constant sea level", fluxes, G, CS%US, haloshift=0)
     endif
diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90
index 761da880b4..c2f188bc6f 100644
--- a/src/initialization/MOM_fixed_initialization.F90
+++ b/src/initialization/MOM_fixed_initialization.F90
@@ -169,18 +169,20 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
 end subroutine MOM_initialize_fixed
 
 !> MOM_initialize_topography makes the appropriate call to set up the bathymetry.  At this
-!! point the topography is in units of m, but this can be changed later.
+!! point the topography is in units of [m], but this can be changed later.
 subroutine MOM_initialize_topography(D, max_depth, G, PF, US)
   type(dyn_horgrid_type),           intent(in)  :: G  !< The dynamic horizontal grid type
   real, dimension(G%isd:G%ied,G%jsd:G%jed), &
-                                    intent(out) :: D  !< Ocean bottom depth in m
+                                    intent(out) :: D  !< Ocean bottom depth [m]
   type(param_file_type),            intent(in)  :: PF !< Parameter file structure
-  real,                             intent(out) :: max_depth !< Maximum depth of model in m
+  real,                             intent(out) :: max_depth !< Maximum depth of model [m]
   type(unit_scale_type),  optional, intent(in)  :: US !< A dimensional unit scaling type
 
-!  This subroutine makes the appropriate call to set up the bottom depth.
-!  This is a separate subroutine so that it can be made public and shared with
-!  the ice-sheet code or other components.
+  ! This subroutine makes the appropriate call to set up the bottom depth.
+  ! This is a separate subroutine so that it can be made public and shared with
+  ! the ice-sheet code or other components.
+
+  ! Local variables
   real :: m_to_Z, Z_to_m  ! Dimensional rescaling factors
   character(len=40)  :: mdl = "MOM_initialize_topography" ! This subroutine's name.
   character(len=200) :: config
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 56c68001b8..d426550762 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -707,7 +707,7 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h)
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: h   !< Layer thicknesses [H ~> m or kg m-2]
   ! Local variables
   integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
-  real :: hTolerance = 0.1 !<  Tolerance to exceed adjustment criteria (m)
+  real :: hTolerance = 0.1 !<  Tolerance to exceed adjustment criteria [Z ~> m]
   real :: hTmp, eTmp, dilate
   character(len=100) :: mesg
 
@@ -1164,7 +1164,7 @@ subroutine cut_off_column_top(nk, tv, GV, G_earth, depth, min_thickness, &
   real, dimension(nk),   intent(in)    :: S_t !< Salinity at top of layer [ppt]
   real, dimension(nk),   intent(in)    :: S_b !< Salinity at bottom of layer [ppt]
   real,                  intent(in)    :: p_surf !< Imposed pressure on ocean at surface [Pa]
-  real, dimension(nk),   intent(inout) :: h   !< Layer thickness (H units, m or Pa)
+  real, dimension(nk),   intent(inout) :: h   !< Layer thickness [H ~> m or kg m-2]
   type(remapping_CS),    pointer       :: remap_CS !< Remapping structure for remapping T and S,
                                                    !! if associated
   real,        optional, intent(in)    :: z_tol !< The tolerance with which to find the depth
diff --git a/src/initialization/midas_vertmap.F90 b/src/initialization/midas_vertmap.F90
index 9277d42f21..a985cf2982 100644
--- a/src/initialization/midas_vertmap.F90
+++ b/src/initialization/midas_vertmap.F90
@@ -714,7 +714,7 @@ end subroutine meshgrid
 !! in each region is an approximation to del2(zi)=0 subject to
 !! boundary conditions along the valid points curve bounding this region.
 subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n)
-  real, dimension(:,:), intent(inout) :: zi !< interface positions (m)
+  real, dimension(:,:), intent(inout) :: zi !< interface positions [m] or arbitrary
   integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< points to be smoothed
   integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< ignore these points
   real, intent(in)  :: sor !< successive over-relaxation coefficient (typically 0.6)
diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90
index 010e531973..4555a52f1d 100644
--- a/src/user/BFB_surface_forcing.F90
+++ b/src/user/BFB_surface_forcing.F90
@@ -24,27 +24,20 @@ module BFB_surface_forcing
 !> Control structure for BFB_surface_forcing
 type, public :: BFB_surface_forcing_CS ; private
 
-  logical :: use_temperature !< If true, temperature and salinity are used as
-                             !! state variables.
+  logical :: use_temperature !< If true, temperature and salinity are used as state variables.
   logical :: restorebuoy     !< If true, use restoring surface buoyancy forcing.
-  real :: Rho0               !<   The density used in the Boussinesq
-                             !! approximation [kg m-3].
-  real :: G_Earth            !<   The gravitational acceleration [m s-2]
-  real :: Flux_const         !<   The restoring rate at the surface [m s-1].
-  real :: gust_const         !<   A constant unresolved background gustiness
+  real :: Rho0               !< The density used in the Boussinesq approximation [kg m-3].
+  real :: G_Earth            !< The gravitational acceleration [m s-2]
+  real :: Flux_const         !< The restoring rate at the surface [m s-1].
+  real :: gust_const         !< A constant unresolved background gustiness
                              !! that contributes to ustar [Pa].
-  real :: SST_s              !< SST at the southern edge of the linear
-                             !! forcing ramp
-  real :: SST_n              !< SST at the northern edge of the linear
-                             !! forcing ramp
-  real :: lfrslat            !< Southern latitude where the linear forcing ramp
-                             !! begins
-  real :: lfrnlat            !< Northern latitude where the linear forcing ramp
-                             !! ends
-  real :: drho_dt            !< Rate of change of density with temperature.
-                             !! Note that temperature is being used as a dummy
-                             !! variable here. All temperatures are converted
-                             !! into density.
+  real :: SST_s              !< SST at the southern edge of the linear forcing ramp [degC]
+  real :: SST_n              !< SST at the northern edge of the linear forcing ramp [degC]
+  real :: lfrslat            !< Southern latitude where the linear forcing ramp begins [degLat]
+  real :: lfrnlat            !< Northern latitude where the linear forcing ramp ends [degLat]
+  real :: drho_dt            !< Rate of change of density with temperature [kg m-3 degC-1].
+                             !!   Note that temperature is being used as a dummy variable here.
+                             !! All temperatures are converted into density.
 
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                              !! regulate the timing of diagnostic output.
@@ -102,7 +95,7 @@ subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, CS)
     ! Set whichever fluxes are to be used here.  Any fluxes that
     ! are always zero do not need to be changed here.
     do j=js,je ; do i=is,ie
-      ! Fluxes of fresh water through the surface are in units of kg m-2 s-1
+      ! Fluxes of fresh water through the surface are in units of [kg m-2 s-1]
       ! and are positive downward - i.e. evaporation should be negative.
       fluxes%evap(i,j) = -0.0 * G%mask2dT(i,j)
       fluxes%lprec(i,j) = 0.0 * G%mask2dT(i,j)
@@ -110,7 +103,7 @@ subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, CS)
       ! vprec will be set later, if it is needed for salinity restoring.
       fluxes%vprec(i,j) = 0.0
 
-      !   Heat fluxes are in units of W m-2 and are positive into the ocean.
+      ! Heat fluxes are in units of [W m-2] and are positive into the ocean.
       fluxes%lw(i,j) = 0.0 * G%mask2dT(i,j)
       fluxes%latent(i,j) = 0.0 * G%mask2dT(i,j)
       fluxes%sens(i,j) = 0.0 * G%mask2dT(i,j)
diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90
index 365e10fdbf..b81061ab29 100644
--- a/src/user/DOME2d_initialization.F90
+++ b/src/user/DOME2d_initialization.F90
@@ -225,8 +225,8 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, param_file,
                      eqn_of_state, just_read_params)
   type(ocean_grid_type),                     intent(in)  :: G !< Ocean grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< The ocean's vertical grid structure.
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC)
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (ppt)
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature [degC]
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity [ppt]
   real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in)  :: h !< Layer thickness [H ~> m or kg m-2]
   type(param_file_type),                     intent(in)  :: param_file !< Parameter file structure
   type(EOS_type),                            pointer     :: eqn_of_state !< Equation of state structure
diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90
index 260512ecfa..fedd46ab03 100644
--- a/src/user/MOM_wave_interface.F90
+++ b/src/user/MOM_wave_interface.F90
@@ -69,7 +69,7 @@ module MOM_wave_interface
 
   ! Surface Wave Dependent 1d/2d/3d vars
   real, allocatable, dimension(:), public :: &
-       WaveNum_Cen        !< Wavenumber bands for read/coupled (1/m)
+       WaveNum_Cen        !< Wavenumber bands for read/coupled [m-1]
   real, allocatable, dimension(:), public :: &
        Freq_Cen           !< Frequency bands for read/coupled [s-1]
   real, allocatable, dimension(:), public :: &
@@ -77,11 +77,11 @@ module MOM_wave_interface
   real, allocatable, dimension(:), public :: &
        PrescribedSurfStkY !< Surface Stokes drift if prescribed [m s-1]
   real, allocatable, dimension(:,:,:), public :: &
-       Us_x               !< 3d Stokes drift profile (zonal, m/s)
+       Us_x               !< 3d zonal Stokes drift profile [m s-1]
                           !! Horizontal -> U points
                           !! Vertical -> Mid-points
   real, allocatable, dimension(:,:,:), public :: &
-       Us_y               !< 3d Stokes drift profile (meridional, m/s)
+       Us_y               !< 3d meridional Stokes drift profile [m s-1]
                           !! Horizontal -> V points
                           !! Vertical -> Mid-points
   real, allocatable, dimension(:,:), public :: &
@@ -189,7 +189,7 @@ module MOM_wave_interface
 
 !> Initializes parameters related to MOM_wave_interface
 subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag )
-  type(time_type), target, intent(in)    :: Time       !< Time (s)
+  type(time_type), target, intent(in)    :: Time       !< Model time
   type(ocean_grid_type),   intent(inout) :: G          !< Grid structure
   type(verticalGrid_type), intent(in)    :: GV         !< Vertical grid structure
   type(unit_scale_type),   intent(in)    :: US         !< A dimensional unit scaling type
@@ -426,8 +426,8 @@ subroutine Update_Surface_Waves(G, GV, US, Day, dt, CS)
   type(ocean_grid_type), intent(inout) :: G   !< Grid structure
   type(verticalGrid_type), intent(in)  :: GV  !< Vertical grid structure
   type(unit_scale_type),   intent(in)  :: US   !< A dimensional unit scaling type
-  type(time_type),         intent(in)  :: Day !< Time (s)
-  type(time_type),         intent(in)  :: dt  !< Timestep (s)
+  type(time_type),         intent(in)  :: Day !< Current model time
+  type(time_type),         intent(in)  :: dt  !< Timestep as a time-type
   ! Local variables
   integer :: ii, jj, kk, b
   type(time_type) :: Day_Center
@@ -966,7 +966,7 @@ end subroutine get_Langmuir_Number
 !!
 !! Original description:
 !! - This function returns the enhancement factor, given the 10-meter
-!!   wind [m s-1], friction velocity [m s-1] and the boundary layer depth (m).
+!!   wind [m s-1], friction velocity [m s-1] and the boundary layer depth [m].
 !!
 !! Update (Jan/25):
 !! - Converted from function to subroutine, now returns Langmuir number.
@@ -1102,7 +1102,7 @@ subroutine Get_SL_Average_Band( GV, AvgDepth, NB, WaveNumbers, SurfStokes, Avera
   real, intent(in)    :: AvgDepth    !< Depth to average over [Z ~> m].
   integer, intent(in) :: NB          !< Number of bands used
   real, dimension(NB), &
-       intent(in)     :: WaveNumbers !< Wavenumber corresponding to each band (1/Z)
+       intent(in)     :: WaveNumbers !< Wavenumber corresponding to each band [Z-1 ~> m-1]
   real, dimension(NB), &
        intent(in)     :: SurfStokes  !< Surface Stokes drift for each band [m s-1]
   real, intent(out)   :: Average     !< Output average Stokes drift over depth AvgDepth [m s-1]
diff --git a/src/user/baroclinic_zone_initialization.F90 b/src/user/baroclinic_zone_initialization.F90
index 10ddf57ff2..8f3ad67ca9 100644
--- a/src/user/baroclinic_zone_initialization.F90
+++ b/src/user/baroclinic_zone_initialization.F90
@@ -33,15 +33,15 @@ subroutine bcz_params(G, GV, US, param_file, S_ref, dSdz, delta_S, dSdx, T_ref,
   type(verticalGrid_type), intent(in)  :: GV         !< The ocean's vertical grid structure.
   type(unit_scale_type),   intent(in)  :: US    !< A dimensional unit scaling type
   type(param_file_type),   intent(in)  :: param_file !< Parameter file handle
-  real,                    intent(out) :: S_ref      !< Reference salinity (ppt)
-  real,                    intent(out) :: dSdz       !< Salinity stratification (ppt/Z)
-  real,                    intent(out) :: delta_S    !< Salinity difference across baroclinic zone (ppt)
-  real,                    intent(out) :: dSdx       !< Linear salinity gradient (ppt/m)
-  real,                    intent(out) :: T_ref      !< Reference temperature (ppt)
-  real,                    intent(out) :: dTdz       !< Temperature stratification (ppt/Z)
-  real,                    intent(out) :: delta_T    !< Temperature difference across baroclinic zone (ppt)
-  real,                    intent(out) :: dTdx       !< Linear temperature gradient (ppt/m)
-  real,                    intent(out) :: L_zone     !< Width of baroclinic zone (m)
+  real,                    intent(out) :: S_ref      !< Reference salinity [ppt]
+  real,                    intent(out) :: dSdz       !< Salinity stratification [ppt Z-1 ~> ppt m-1]
+  real,                    intent(out) :: delta_S    !< Salinity difference across baroclinic zone [ppt]
+  real,                    intent(out) :: dSdx       !< Linear salinity gradient [ppt m-1]
+  real,                    intent(out) :: T_ref      !< Reference temperature [degC]
+  real,                    intent(out) :: dTdz       !< Temperature stratification [degC Z-1 ~> degC m-1]
+  real,                    intent(out) :: delta_T    !< Temperature difference across baroclinic zone [degC]
+  real,                    intent(out) :: dTdx       !< Linear temperature gradient [degC m-1]
+  real,                    intent(out) :: L_zone     !< Width of baroclinic zone [m]
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
 
diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90
index ff5e7a1619..095884d232 100644
--- a/src/user/dumbbell_surface_forcing.F90
+++ b/src/user/dumbbell_surface_forcing.F90
@@ -26,11 +26,10 @@ module dumbbell_surface_forcing
   logical :: use_temperature !< If true, temperature and salinity are used as
                              !! state variables.
   logical :: restorebuoy     !< If true, use restoring surface buoyancy forcing.
-  real :: Rho0               !<   The density used in the Boussinesq
-                             !! approximation [kg m-3].
-  real :: G_Earth            !<   The gravitational acceleration [m s-2]
-  real :: Flux_const         !<   The restoring rate at the surface [m s-1].
-  real :: gust_const         !<   A constant unresolved background gustiness
+  real :: Rho0               !< The density used in the Boussinesq approximation [kg m-3].
+  real :: G_Earth            !< The gravitational acceleration [m s-2]
+  real :: Flux_const         !< The restoring rate at the surface [m s-1].
+  real :: gust_const         !< A constant unresolved background gustiness
                              !! that contributes to ustar [Pa].
   real :: slp_amplitude      !< The amplitude of pressure loading [Pa] applied
                              !! to the reservoirs
@@ -98,7 +97,7 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, CS)
     ! Set whichever fluxes are to be used here.  Any fluxes that
     ! are always zero do not need to be changed here.
     do j=js,je ; do i=is,ie
-      ! Fluxes of fresh water through the surface are in units of kg m-2 s-1
+      ! Fluxes of fresh water through the surface are in units of [kg m-2 s-1]
       ! and are positive downward - i.e. evaporation should be negative.
       fluxes%evap(i,j) = -0.0 * G%mask2dT(i,j)
       fluxes%lprec(i,j) = 0.0 * G%mask2dT(i,j)
@@ -106,7 +105,7 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, CS)
       ! vprec will be set later, if it is needed for salinity restoring.
       fluxes%vprec(i,j) = 0.0
 
-      !   Heat fluxes are in units of W m-2 and are positive into the ocean.
+      ! Heat fluxes are in units of [W m-2] and are positive into the ocean.
       fluxes%lw(i,j) = 0.0 * G%mask2dT(i,j)
       fluxes%latent(i,j) = 0.0 * G%mask2dT(i,j)
       fluxes%sens(i,j) = 0.0 * G%mask2dT(i,j)
diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90
index 2c7dc2f530..6180ff2e00 100644
--- a/src/user/seamount_initialization.F90
+++ b/src/user/seamount_initialization.F90
@@ -195,8 +195,8 @@ subroutine seamount_initialize_temperature_salinity ( T, S, h, G, GV, param_file
                                                   eqn_of_state, just_read_params)
   type(ocean_grid_type),                     intent(in)  :: G !< Ocean grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< Vertical grid structure
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC)
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (ppt)
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature [degC]
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity [ppt]
   real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in)  :: h !< Layer thickness [H ~> m or kg m-2]
   type(param_file_type),                     intent(in)  :: param_file !< Parameter file structure
   type(EOS_type),                            pointer     :: eqn_of_state !< Equation of state structure
diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90
index 37d7c1b8c9..990d43fda4 100644
--- a/src/user/sloshing_initialization.F90
+++ b/src/user/sloshing_initialization.F90
@@ -180,8 +180,8 @@ subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, GV, param_file
                                                       eqn_of_state, just_read_params)
   type(ocean_grid_type),                     intent(in)  :: G !< Ocean grid structure.
   type(verticalGrid_type),                   intent(in)  :: GV !< The ocean's vertical grid structure.
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC).
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (ppt).
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature [degC].
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity [ppt].
   real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in)  :: h !< Layer thickness [H ~> m or kg m-2].
   type(param_file_type),                     intent(in)  :: param_file !< A structure indicating the
                                                             !! open file to parse for model
diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90
index 4d9fb72a4e..5a29614506 100644
--- a/src/user/user_change_diffusivity.F90
+++ b/src/user/user_change_diffusivity.F90
@@ -66,8 +66,8 @@ subroutine user_change_diff(h, tv, G, GV, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_a
   ! Local variables
   real :: Rcv(SZI_(G),SZK_(G)) ! The coordinate density in layers [kg m-3].
   real :: p_ref(SZI_(G))       ! An array of tv%P_Ref pressures.
-  real :: rho_fn      ! The density dependence of the input function, 0-1, ND.
-  real :: lat_fn      ! The latitude dependence of the input function, 0-1, ND.
+  real :: rho_fn      ! The density dependence of the input function, 0-1 [nondim].
+  real :: lat_fn      ! The latitude dependence of the input function, 0-1 [nondim].
   logical :: use_EOS  ! If true, density is calculated from T & S using an
                       ! equation of state.
   logical :: store_Kd_add  ! Save the added diffusivity as a diagnostic if true.
diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90
index 7e6fc1f07e..1de9c3664a 100644
--- a/src/user/user_initialization.F90
+++ b/src/user/user_initialization.F90
@@ -137,8 +137,8 @@ end subroutine USER_initialize_velocity
 !! into T(:,:,:) and S(:,:,:).
 subroutine USER_init_temperature_salinity(T, S, G, param_file, eqn_of_state, just_read_params)
   type(ocean_grid_type),                     intent(in)  :: G !< Ocean grid structure.
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC).
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (ppt).
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature [degC].
+  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity [ppt].
   type(param_file_type),                     intent(in)  :: param_file !< A structure indicating the
                                                             !! open file to parse for model
                                                             !! parameter values.