diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
index f732cb9a56..88d2cb3f42 100644
--- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
+++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
@@ -83,14 +83,14 @@ module MOM_surface_forcing_gfdl
                                 !! type without any further adjustments to drive the ocean dynamics.
                                 !! The actual net mass source may differ due to corrections.
 
-  real :: gust_const            !< Constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
+  real :: gust_const            !< Constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa]
   logical :: read_gust_2d       !< If true, use a 2-dimensional gustiness supplied from an input file.
   real, pointer, dimension(:,:) :: &
     TKE_tidal => NULL()         !< Turbulent kinetic energy introduced to the bottom boundary layer
                                 !! by drag on the tidal flows [R Z3 T-3 ~> W m-2].
   real, pointer, dimension(:,:) :: &
     gust => NULL()              !< A spatially varying unresolved background gustiness that
-                                !! contributes to ustar [R L Z T-1 ~> Pa].  gust is used when read_gust_2d is true.
+                                !! contributes to ustar [R L Z T-2 ~> Pa].  gust is used when read_gust_2d is true.
   real, pointer, dimension(:,:) :: &
     ustar_tidal => NULL()       !< Tidal contribution to the bottom friction velocity [Z T-1 ~> m s-1]
   real :: cd_tides              !< Drag coefficient that applies to the tides [nondim]
diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
index b34f1a3b35..0364d46ddc 100644
--- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
+++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
@@ -76,14 +76,14 @@ module MOM_surface_forcing_mct
                                 !! the correction for the atmospheric (and sea-ice)
                                 !! pressure limited by max_p_surf instead of the
                                 !! full atmospheric pressure.  The default is true.
-  real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
+  real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa]
   logical :: read_gust_2d       !< If true, use a 2-dimensional gustiness supplied
                                 !! from an input file.
   real, pointer, dimension(:,:) :: &
     TKE_tidal => NULL(), &      !< turbulent kinetic energy introduced to the
                                 !! bottom boundary layer by drag on the tidal flows [R Z3 T-3 ~> W m-2]
     gust => NULL(), &           !< spatially varying unresolved background
-                                !! gustiness that contributes to ustar [R L Z T-1 ~> Pa].
+                                !! gustiness that contributes to ustar [R L Z T-2 ~> Pa].
                                 !! gust is used when read_gust_2d is true.
     ustar_tidal => NULL()       !< tidal contribution to the bottom friction velocity [Z T-1 ~> m s-1]
   real :: cd_tides              !< drag coefficient that applies to the tides (nondimensional)
diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
index e2710bd18c..921d43d9b8 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
@@ -81,14 +81,14 @@ module MOM_surface_forcing_nuopc
   logical :: use_CFC            !< enables the MOM_CFC_cap tracer package.
   logical :: enthalpy_cpl       !< Controls if enthalpy terms are provided by the coupler or computed
                                 !! internally.
-  real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
+  real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa]
   logical :: read_gust_2d       !< If true, use a 2-dimensional gustiness supplied
                                 !! from an input file.
   real, pointer, dimension(:,:) :: &
     TKE_tidal => NULL(), &      !< turbulent kinetic energy introduced to the
                                 !! bottom boundary layer by drag on the tidal flows [R Z3 T-3 ~> W m-2]
     gust => NULL(), &           !< spatially varying unresolved background
-                                !! gustiness that contributes to ustar [R L Z T-1 ~> Pa].
+                                !! gustiness that contributes to ustar [R L Z T-2 ~> Pa].
                                 !! gust is used when read_gust_2d is true.
     ustar_tidal => NULL()       !< tidal contribution to the bottom friction velocity [Z T-1 ~> m s-1]
   real :: cd_tides              !< drag coefficient that applies to the tides (nondimensional)
diff --git a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 b/config_src/drivers/solo_driver/MESO_surface_forcing.F90
index 18c3c33fdb..12f1b6b78d 100644
--- a/config_src/drivers/solo_driver/MESO_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/MESO_surface_forcing.F90
@@ -31,7 +31,7 @@ module MESO_surface_forcing
   real :: G_Earth            !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
   real :: Flux_const         !< The restoring rate at the surface [Z T-1 ~> m s-1].
   real :: gust_const         !< A constant unresolved background gustiness
-                             !! that contributes to ustar [R L Z T-1 ~> Pa]
+                             !! that contributes to ustar [R L Z T-2 ~> Pa]
   real, dimension(:,:), pointer :: &
     T_Restore(:,:) => NULL(), & !< The temperature to restore the SST toward [C ~> degC].
     S_Restore(:,:) => NULL(), & !< The salinity to restore the sea surface salnity toward [S ~> ppt]
diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90
index 74ed5c3023..092bc9e513 100644
--- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90
@@ -85,13 +85,13 @@ module MOM_surface_forcing
   real :: latent_heat_fusion    !< latent heat of fusion times [Q ~> J kg-1]
   real :: latent_heat_vapor     !< latent heat of vaporization [Q ~> J kg-1]
   real :: tau_x0                !< Constant zonal wind stress used in the WIND_CONFIG="const"
-                                !! forcing [R L Z T-1 ~> Pa]
+                                !! forcing [R L Z T-2 ~> Pa]
   real :: tau_y0                !< Constant meridional wind stress used in the WIND_CONFIG="const"
-                                !! forcing [R L Z T-1 ~> Pa]
+                                !! forcing [R L Z T-2 ~> Pa]
 
-  real    :: gust_const                 !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
+  real    :: gust_const                 !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa]
   logical :: read_gust_2d               !< if true, use 2-dimensional gustiness supplied from a file
-  real, pointer :: gust(:,:) => NULL()  !< spatially varying unresolved background gustiness [R L Z T-1 ~> Pa]
+  real, pointer :: gust(:,:) => NULL()  !< spatially varying unresolved background gustiness [R L Z T-2 ~> Pa]
                                         !! gust is used when read_gust_2d is true.
 
   real, pointer :: T_Restore(:,:)    => NULL()  !< temperature to damp (restore) the SST to [C ~> degC]
@@ -102,9 +102,9 @@ module MOM_surface_forcing
 
   ! if WIND_CONFIG=='gyres' then use the following as  = A, B, C and n respectively for
   ! taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L)
-  real :: gyres_taux_const   !< A constant wind stress [R L Z T-1 ~> Pa].
-  real :: gyres_taux_sin_amp !< The amplitude of cosine wind stress gyres [R L Z T-1 ~> Pa], if WIND_CONFIG=='gyres'
-  real :: gyres_taux_cos_amp !< The amplitude of cosine wind stress gyres [R L Z T-1 ~> Pa], if WIND_CONFIG=='gyres'
+  real :: gyres_taux_const   !< A constant wind stress [R L Z T-2 ~> Pa].
+  real :: gyres_taux_sin_amp !< The amplitude of cosine wind stress gyres [R L Z T-2 ~> Pa], if WIND_CONFIG=='gyres'
+  real :: gyres_taux_cos_amp !< The amplitude of cosine wind stress gyres [R L Z T-2 ~> Pa], if WIND_CONFIG=='gyres'
   real :: gyres_taux_n_pis   !< The number of sine lobes in the basin if WIND_CONFIG=='gyres' [nondim]
   integer :: answer_date     !< This 8-digit integer gives the approximate date with which the order
                              !! of arithmetic and expressions were added to the code.
@@ -115,7 +115,7 @@ module MOM_surface_forcing
                                             !! gustless wind friction velocity.
   ! if WIND_CONFIG=='scurves' then use the following to define a piecewise scurve profile
   real :: scurves_ydata(20) = 90. !< Latitudes of scurve nodes [degreesN]
-  real :: scurves_taux(20) = 0.   !< Zonal wind stress values at scurve nodes [R L Z T-1 ~> Pa]
+  real :: scurves_taux(20) = 0.   !< Zonal wind stress values at scurve nodes [R L Z T-2 ~> Pa]
 
   real :: T_north   !< Target temperatures at north used in buoyancy_forcing_linear [C ~> degC]
   real :: T_south   !< Target temperatures at south used in buoyancy_forcing_linear [C ~> degC]
@@ -392,7 +392,7 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS)
 
   mag_tau = sqrt( tau_x0**2 + tau_y0**2)
 
-  ! Set the steady surface wind stresses, in units of [R L Z T-1 ~> Pa].
+  ! Set the steady surface wind stresses, in units of [R L Z T-2 ~> Pa].
   do j=js,je ; do I=is-1,Ieq
     forces%taux(I,j) = tau_x0
   enddo ; enddo
@@ -438,7 +438,7 @@ subroutine wind_forcing_2gyre(sfc_state, forces, day, G, US, CS)
   Pa_to_RLZ_T2 = US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z
   PI = 4.0*atan(1.0)
 
-  ! Set the steady surface wind stresses, in units of [R L Z T-1 ~> Pa].
+  ! Set the steady surface wind stresses, in units of [R L Z T-2 ~> Pa].
   do j=js,je ; do I=is-1,Ieq
     forces%taux(I,j) = 0.1 * Pa_to_RLZ_T2 * &
                       (1.0 - cos(2.0*PI*(G%geoLatCu(I,j)-CS%South_lat) / CS%len_lat))
@@ -513,7 +513,7 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)
 
   PI = 4.0*atan(1.0)
 
-  ! steady surface wind stresses [R L Z T-1 ~> Pa]
+  ! steady surface wind stresses [R L Z T-2 ~> Pa]
   do j=js-1,je+1 ; do I=is-1,Ieq
     y = (G%geoLatCu(I,j)-CS%South_lat) / CS%len_lat
     forces%taux(I,j) = CS%gyres_taux_const + &
@@ -670,8 +670,8 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
                                                   !! a previous surface_forcing_init call
   ! Local variables
   character(len=200) :: filename  ! The name of the input file.
-  real    :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal wind stresses at h-points [R L Z T-1 ~> Pa]
-  real    :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R L Z T-1 ~> Pa]
+  real    :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal wind stresses at h-points [R L Z T-2 ~> Pa]
+  real    :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R L Z T-2 ~> Pa]
   real    :: Pa_to_RLZ_T2            ! A unit conversion factor from Pa to the internal units
                                      ! for wind stresses [R Z L T-2 Pa-1 ~> 1]
   integer :: time_lev_daily          ! The time levels to read for fields with
diff --git a/config_src/drivers/solo_driver/user_surface_forcing.F90 b/config_src/drivers/solo_driver/user_surface_forcing.F90
index ae3f854335..fc803c27e6 100644
--- a/config_src/drivers/solo_driver/user_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/user_surface_forcing.F90
@@ -33,10 +33,10 @@ module user_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 [R ~> kg m-3].
-  real :: G_Earth            !< The gravitational acceleration [L2 Z-1 s-2 ~> m s-2].
+  real :: G_Earth            !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
   real :: Flux_const         !< The restoring rate at the surface [Z T-1 ~> m s-1].
   real :: gust_const         !< A constant unresolved background gustiness
-                             !! that contributes to ustar [R L Z T-1 ~> Pa].
+                             !! that contributes to ustar [R L Z T-2 ~> Pa].
 
   type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                              !! timing of diagnostic output.
@@ -71,7 +71,7 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS)
   ! Allocate the forcing arrays, if necessary.
   call allocate_mech_forcing(G, forces, stress=.true., ustar=.true.)
 
-  !  Set the surface wind stresses, in units of [R L Z T-1 ~> Pa].  A positive taux
+  !  Set the surface wind stresses, in units of [R L Z T-2 ~> Pa].  A positive taux
   !  accelerates the ocean to the (pseudo-)east.
 
   !  The i-loop extends to is-1 so that taux can be used later in the
diff --git a/config_src/drivers/unit_drivers/MOM_sum_driver.F90 b/config_src/drivers/unit_drivers/MOM_sum_driver.F90
index 0a6191b286..7a1ba82843 100644
--- a/config_src/drivers/unit_drivers/MOM_sum_driver.F90
+++ b/config_src/drivers/unit_drivers/MOM_sum_driver.F90
@@ -183,7 +183,7 @@ subroutine benchmark_init_topog_local(D, G, param_file, max_depth)
   real :: PI                   ! 3.1415926... calculated as 4*atan(1) [nondim]
   real :: D0                   ! A constant to make the maximum
                                ! basin depth MAXIMUM_DEPTH [m]
-  real :: m_to_Z  ! A dimensional rescaling factor [m ~> Z]
+  real :: m_to_Z  ! A dimensional rescaling factor [Z m-1 ~> 1]
   real :: x ! A fractional position in the x-direction [nondim]
   real :: y ! A fractional position in the y-direction [nondim]
   ! This include declares and sets the variable "version".
diff --git a/src/ALE/MOM_hybgen_remap.F90 b/src/ALE/MOM_hybgen_remap.F90
index 213c6c677e..5ab3e162db 100644
--- a/src/ALE/MOM_hybgen_remap.F90
+++ b/src/ALE/MOM_hybgen_remap.F90
@@ -263,7 +263,7 @@ subroutine hybgen_weno_coefs(s, h_src, edges, nk, ns, thin, PCM_lay)
 !  real, parameter :: dsmll=1.0e-8  ! This has units of [A2], and hence can not be a parameter.
 !
   real :: curv_cell   ! An estimate of the tracer curvature centered on a cell times the grid
-                      ! spacing [A H-1 ~> A m-1 or A kg m-2]
+                      ! spacing [A H-1 ~> A m-1 or A m2 kg-1]
   real :: seh1, seh2  ! Tracer slopes at the cell edges times the cell grid spacing [A]
   real :: q01, q02    ! Various tracer differences between a cell average and the edge values [A]
   real :: q001, q002  ! Tracer slopes at the cell edges times the cell grid spacing [A]
@@ -277,7 +277,7 @@ subroutine hybgen_weno_coefs(s, h_src, edges, nk, ns, thin, PCM_lay)
                       ! concentrations and the left and right edges [A2]
   real :: min_ratio   ! The minimum ratio of the values of zw used to interpolate the edge values [nondim]
   real :: wt1         ! The weight of the upper layer in the interpolated shared edge value [nondim]
-  real :: slope_edge(nk+1)  ! Tracer slopes at the edges [A H-1 ~> A m-1 or A kg m-2]
+  real :: slope_edge(nk+1)  ! Tracer slopes at the edges [A H-1 ~> A m-1 or A m2 kg-1]
   real :: val_edge(nk+1)    ! A weighted average edge concentration [A]
   integer :: i, k
 
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index 00e748bdc1..84eb5fc90a 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -1976,7 +1976,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   real :: salin_underflow      ! A tiny value of salinity below which the it is set to 0 [S ~> ppt]
   real :: temp_underflow       ! A tiny magnitude of temperatures below which they are set to 0 [C ~> degC]
   real :: conv2watt            ! A conversion factor from temperature fluxes to heat
-                               ! fluxes [J m-2 H-1 degC-1 ~> J m-3 degC-1 or J kg-1 degC-1]
+                               ! fluxes [J m-2 H-1 C-1 ~> J m-3 degC-1 or J kg-1 degC-1]
   real :: conv2salt            ! A conversion factor for salt fluxes [m H-1 ~> 1] or [kg m-2 H-1 ~> 1]
   real :: RL2_T2_rescale, Z_rescale, QRZ_rescale ! Unit conversion factors
   character(len=48) :: S_flux_units
diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90
index 54eecd20c3..090d1ee0fb 100644
--- a/src/core/MOM_continuity_PPM.F90
+++ b/src/core/MOM_continuity_PPM.F90
@@ -1044,7 +1044,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(in)    :: h_in !< Layer thickness used to
                                                                     !! calculate fluxes [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out)   :: vh   !< Volume flux through meridional
-                                                                    !! faces = v*h*dx [H L2 s-1 ~> m3 s-1 or kg s-1]
+                                                                    !! faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]
   real,                                       intent(in)    :: dt   !< Time increment [T ~> s].
   type(unit_scale_type),                      intent(in)    :: US   !< A dimensional unit scaling type
   type(continuity_PPM_CS),                    intent(in)    :: CS   !< This module's control structure.G
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index f005ac1a0e..9a5e1f48f5 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -969,7 +969,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
                                                       ! [H T-1 ~> m s-1 or kg m-2 s-1]
   real, dimension(SZI_(G))              :: netHeat    ! net temp flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
   real, dimension(max(nsw,1), SZI_(G))  :: penSWbnd   ! penetrating SW radiation by band
-                                                      ! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
+                                                      ! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
   real, dimension(SZI_(G))              :: pressure   ! pressure at the surface [R L2 T-2 ~> Pa]
   real, dimension(SZI_(G))              :: dRhodT     ! density partial derivative wrt temp [R C-1 ~> kg m-3 degC-1]
   real, dimension(SZI_(G))              :: dRhodS     ! density partial derivative wrt saln [R S-1 ~> kg m-3 ppt-1]
@@ -996,7 +996,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
   ! The surface forcing is contained in the fluxes type.
   ! We aggregate the thermodynamic forcing for a time step into the following:
   ! netH       = water added/removed via surface fluxes [H T-1 ~> m s-1 or kg m-2 s-1]
-  ! netHeat    = heat via surface fluxes [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
+  ! netHeat    = heat via surface fluxes [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
   ! netSalt    = salt via surface fluxes [S H T-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 (in arbitrary time units)
@@ -1015,12 +1015,12 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
                                 tv%eqn_of_state, EOS_domain(G%HI))
 
   ! Adjust netSalt to reflect dilution effect of FW flux
-  ! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
+  ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
   netSalt(G%isc:G%iec) = netSalt(G%isc:G%iec) - Salt(G%isc:G%iec,j,1) * netH(G%isc:G%iec)
 
   ! Add in the SW heating for purposes of calculating the net
   ! surface buoyancy flux affecting the top layer.
-  ! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
+  ! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
   !netHeat(:) = netHeatMinusSW(:) + sum( penSWbnd, dim=1 )
   netHeat(G%isc:G%iec) = netHeatMinusSW(G%isc:G%iec) + netPen(G%isc:G%iec,1)
 
diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90
index 4e812b65d7..ebe3907469 100644
--- a/src/core/MOM_porous_barriers.F90
+++ b/src/core/MOM_porous_barriers.F90
@@ -285,7 +285,7 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt)
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(out) :: eta_v !< Layer interface heights at v points [Z ~> m]
 
   ! local variables
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Layer interface heights [Z ~> m or 1/eta_to_m].
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Layer interface heights [Z ~> m].
   real :: h_neglect ! Negligible thicknesses [Z ~> m]
   integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq
 
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 35cdf3038a..6aa94f584f 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -101,7 +101,7 @@ module MOM_variables
   real, dimension(:,:), pointer :: salt_deficit => NULL()
                          !<   The salt needed to maintain the ocean column
                          !! at a minimum salinity of MIN_SALINITY since the last time
-                         !! that calculate_surface_state was called, [ppt R Z ~> gSalt m-2].
+                         !! that calculate_surface_state was called, [S R Z ~> gSalt m-2].
   real, dimension(:,:), pointer :: TempxPmE => NULL()
                          !<   The net inflow of water into the ocean times the
                          !! temperature at which this inflow occurs since the
diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90
index 12b87ebc64..4ddedf85a8 100644
--- a/src/equation_of_state/MOM_EOS.F90
+++ b/src/equation_of_state/MOM_EOS.F90
@@ -632,11 +632,11 @@ end subroutine calc_spec_vol_1d
 
 !> Calls the appropriate subroutine to calculate the freezing point for scalar inputs.
 subroutine calculate_TFreeze_scalar(S, pressure, T_fr, EOS, pres_scale, scale_from_EOS)
-  real,           intent(in)  :: S    !< Salinity, [ppt] or [Z ~> ppt] depending on scale_from_EOS
+  real,           intent(in)  :: S    !< Salinity, [ppt] or [S ~> ppt] depending on scale_from_EOS
   real,           intent(in)  :: pressure !< Pressure, in [Pa] or [R L2 T-2 ~> Pa] depending on
                                       !! pres_scale or scale_from_EOS
   real,           intent(out) :: T_fr !< Freezing point potential temperature referenced to the
-                                      !! surface [degC] or [degC ~> C] depending on scale_from_EOS
+                                      !! surface [degC] or [C ~> degC] depending on scale_from_EOS
   type(EOS_type), intent(in)  :: EOS  !< Equation of state structure
   real, optional, intent(in)  :: pres_scale  !< A multiplicative factor to convert pressure
                                       !! into Pa [Pa T2 R-1 L-2 ~> 1].
diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90
index 07464861cb..dd45e6cd81 100644
--- a/src/equation_of_state/MOM_EOS_linear.F90
+++ b/src/equation_of_state/MOM_EOS_linear.F90
@@ -331,7 +331,7 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, &
                         intent(in)  :: T         !< Potential temperature relative to the surface
                                                  !! [C ~> degC].
   real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
-                        intent(in)  :: S         !< Salinity [S ?~> PSU].
+                        intent(in)  :: S         !< Salinity [S ~> PSU].
   real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                         intent(in)  :: z_t       !< Height at the top of the layer in depth units [Z ~> m].
   real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
diff --git a/src/framework/MOM_unit_scaling.F90 b/src/framework/MOM_unit_scaling.F90
index caf73ae8e1..bfc2189188 100644
--- a/src/framework/MOM_unit_scaling.F90
+++ b/src/framework/MOM_unit_scaling.F90
@@ -41,7 +41,7 @@ module MOM_unit_scaling
   real :: L_T_to_m_s      !< Convert lateral velocities from L T-1 to m s-1                         [T m L-1 s-1 ~> 1]
   real :: m_s_to_L_T      !< Convert lateral velocities from m s-1 to L T-1                         [L s T-1 m-1 ~> 1]
   real :: L_T2_to_m_s2    !< Convert lateral accelerations from L T-2 to m s-2                     [L s2 T-2 m-1 ~> 1]
-  real :: Z2_T_to_m2_s    !< Convert vertical diffusivities from Z2 T-1 to m2 s-1                 [T1 m2 Z-2 s-1 ~> 1]
+  real :: Z2_T_to_m2_s    !< Convert vertical diffusivities from Z2 T-1 to m2 s-1                  [T m2 Z-2 s-1 ~> 1]
   real :: m2_s_to_Z2_T    !< Convert vertical diffusivities from m2 s-1 to Z2 T-1                  [Z2 s T-1 m-2 ~> 1]
   real :: W_m2_to_QRZ_T   !< Convert heat fluxes from W m-2 to Q R Z T-1                       [Q R Z m2 T-1 W-1 ~> 1]
   real :: QRZ_T_to_W_m2   !< Convert heat fluxes from Q R Z T-1 to W m-2                    [W T Q-1 R-1 Z-1 m-2 ~> 1]
diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
index 6691095b08..3049cae00c 100644
--- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90
+++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
@@ -46,9 +46,9 @@ module MOM_ice_shelf_dynamics
   real, pointer, dimension(:,:) :: v_shelf => NULL() !< the meridional velocity of the ice shelf/sheet
                                        !! on q-points (B grid) [L T-1 ~> m s-1]
   real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the driving stress of the ice shelf/sheet
-                                       !! on q-points (C grid) [Pa ~> Pa]
+                                       !! on q-points (C grid) [R L2 T-2 ~> Pa]
   real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional stress of the ice shelf/sheet
-                                       !! on q-points (C grid) [Pa ~> Pa]
+                                       !! on q-points (C grid) [R L2 T-2 ~> Pa]
   real, pointer, dimension(:,:) :: u_face_mask => NULL() !< mask for velocity boundary conditions on the C-grid
                                        !! u-face - this is because the FEM cares about FACES THAT GET INTEGRATED OVER,
                                        !! not vertices. Will represent boundary conditions on computational boundary
@@ -147,7 +147,7 @@ module MOM_ice_shelf_dynamics
   logical :: moving_shelf_front  !< Specify whether to advance shelf front (and calve).
   logical :: calve_to_mask       !< If true, calve off the ice shelf when it passes the edge of a mask.
   real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m].
-  real :: T_shelf_missing   !< An ice shelf temperature to use where there is no ice shelf [degC ~> C]
+  real :: T_shelf_missing   !< An ice shelf temperature to use where there is no ice shelf [C ~> degC]
   real :: cg_tolerance !< The tolerance in the CG solver, relative to initial residual, that
                        !! determines when to stop the conjugate gradient iterations [nondim].
   real :: nonlinear_tolerance !< The fractional nonlinear tolerance, relative to the initial error,
@@ -235,7 +235,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS)
   type(MOM_restart_CS),   intent(inout) :: restart_CS !< MOM restart control struct
 
   ! Local variables
-  real :: T_shelf_missing ! An ice shelf temperature to use where there is no ice shelf [degC ~> C]
+  real :: T_shelf_missing ! An ice shelf temperature to use where there is no ice shelf [C ~> degC]
   logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics
   character(len=40)  :: mdl = "MOM_ice_shelf_dyn"  ! This module's name.
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
@@ -696,10 +696,10 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled
   logical,      optional, intent(in)    :: coupled_grounding !< If true, the grounding line is
                                               !! determined by coupled ice-ocean dynamics
   logical,      optional, intent(in)    :: must_update_vel !< Always update the ice velocities if true.
-  real, dimension(SZDIB_(G),SZDJB_(G))  ::taud_x,taud_y   !<area-averaged driving stress [R L2T-2 ~> Pa]
-  real, dimension(SZDI_(G),SZDJ_(G))  :: ice_visc ! <area-averaged vertically integrated ice viscosity
-                                              !![R2 L2T-3 ~> Pa s-1 m]
-  real, dimension(SZDI_(G),SZDJ_(G))  :: basal_tr ! <area-averaged basal traction [R L2T-2 ~> Pa]
+  real, dimension(SZDIB_(G),SZDJB_(G))  :: taud_x, taud_y  !<area-averaged driving stress [R L2 T-2 ~> Pa]
+  real, dimension(SZDI_(G),SZDJ_(G))  :: ice_visc !< area-averaged vertically integrated ice viscosity
+                                              !! [R L2 Z T-1 ~> Pa s m]
+  real, dimension(SZDI_(G),SZDJ_(G))  :: basal_tr !< area-averaged basal traction [R L2 T-2 ~> Pa]
   integer :: iters
   logical :: update_ice_vel, coupled_GL
 
@@ -1806,10 +1806,10 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
   real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: OD  !< ocean floor depth at tracer points [Z ~> m].
   real, dimension(SZDIB_(G),SZDJB_(G)), &
-                         intent(inout) :: taudx  !< X-direction driving stress at q-points [kg L s-2 ~> kg m s-2]
+                         intent(inout) :: taudx  !< X-direction driving stress at q-points [R L3 Z T-2 ~> kg m s-2]
   real, dimension(SZDIB_(G),SZDJB_(G)), &
-                         intent(inout) :: taudy  !< Y-direction driving stress at q-points [kg L s-2 ~> kg m s-2]
-                                                  ! This will become  [R L3 Z T-2 ~> kg m s-2]
+                         intent(inout) :: taudy  !< Y-direction driving stress at q-points [R L3 Z T-2 ~> kg m s-2]
+
 
 ! driving stress!
 
@@ -1827,7 +1827,6 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
   real, pointer, dimension(:,:,:,:) :: Phi => NULL() ! The gradients of bilinear basis elements at Gaussian
                                                 ! quadrature points surrounding the cell vertices [L-1 ~> m-1].
 
-
   real    :: rho, rhow, rhoi_rhow ! Ice and ocean densities [R ~> kg m-3]
   real    :: sx, sy    ! Ice shelf top slopes [Z L-1 ~> nondim]
   real    :: neumann_val ! [R Z L2 T-2 ~> kg s-2]
@@ -3200,8 +3199,7 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time)
   real, dimension(SZDI_(G),SZDJ_(G))   :: th_after_uflux, th_after_vflux, TH ! Integrated temperatures [C Z ~> degC m]
   integer                           :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
   real :: Tsurf ! Surface air temperature [C ~> degC].  This is hard coded but should be an input argument.
-  real :: adot  ! A surface heat exchange coefficient divided by the heat capacity of
-                ! ice [R Z T-1 degC-1 ~> kg m-2 s-1 degC-1].
+  real :: adot  ! A surface heat exchange coefficient [R Z T-1 ~> kg m-2 s-1].
 
 
   ! For now adot and Tsurf are defined here adot=surf acc 0.1m/yr, Tsurf=-20oC, vary them later
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 0e50ebb67f..bd0931c694 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -1192,7 +1192,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: S_t, S_b ! Top and bottom edge values for reconstructions
                                                         ! of salinity within each layer [S ~> ppt]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: T_t, T_b ! Top and bottom edge values for reconstructions
-                                                        ! of temperature within each layer [T ~> degC]
+                                                        ! of temperature within each layer [C ~> degC]
   character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
   real :: scale_factor   ! A file-dependent scaling factor for the input pressure [various].
   real :: min_thickness  ! The minimum layer thickness, recast into Z units [Z ~> m].
diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
index 94f4468433..ffdf236152 100644
--- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
+++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
@@ -813,7 +813,7 @@ real function growth_time(u_star, hBL, absf, h_neg, vonKar, Kv_rest, restrat_coe
 
   ! Local variables
   real :: mom_mixrate  ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
-  real :: Kv_eff       ! An effective overall viscosity [Z1 T-1 ~> m2 s-1]
+  real :: Kv_eff       ! An effective overall viscosity [Z2 T-1 ~> m2 s-1]
   real :: pi2          ! A scaling constant that is approximately pi^2 [nondim]
 
   ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) + Kv_water
diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
index 4e12aeeaad..a7ff2f1c0a 100644
--- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90
+++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
@@ -722,7 +722,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
   integer :: nk_linear  ! The number of layers over which the streamfunction goes to 0.
   real :: G_rho0        ! g/Rho0 [L2 R-1 Z-1 T-2 ~> m4 kg-1 s-2].
   real :: N2_floor      ! A floor for N2 to avoid degeneracy in the elliptic solver
-                        ! times unit conversion factors [T-2 L2 Z-2 ~> s-2]
+                        ! times unit conversion factors [L2 Z-2 T-2 ~> s-2]
   real :: Tl(5)         ! copy of T in local stencil [C ~> degC]
   real :: mn_T          ! mean of T in local stencil [C ~> degC]
   real :: mn_T2         ! mean of T**2 in local stencil [C2 ~> degC2]
diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
index 5e530bea3d..66e2dfa6b2 100644
--- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
+++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
@@ -912,7 +912,7 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot,      &
   real, dimension(SZI_(G),SZK0_(GV)), &
                             intent(in)    :: T     !< Layer temperatures [C ~> degC].
   real, dimension(SZI_(G),SZK0_(GV)), &
-                            intent(in)    :: S     !< Layer salinities [C ~> ppt].
+                            intent(in)    :: S     !< Layer salinities [S ~> ppt].
   real, dimension(SZI_(G),SZK0_(GV)), &
                             intent(in)    :: R0    !< Potential density referenced to
                                                    !! surface pressure [R ~> kg m-3].
@@ -1890,7 +1890,7 @@ subroutine resort_ML(h, T, S, R0, Rcv, RcvTgt, eps, d_ea, d_eb, ksort, G, GV, CS
   ! Local variables
   real    :: h_move     ! The thickness of water being moved between layers [H ~> m or kg m-2]
   real    :: h_tgt_old  ! The previous thickness of the recipient layer [H ~> m or kg m-2]
-  real    :: I_hnew     ! The inverse of a new layer thickness [H-1 ~> m-1 or m3 kg-1]
+  real    :: I_hnew     ! The inverse of a new layer thickness [H-1 ~> m-1 or m2 kg-1]
   real    :: dT_dS_wt2  ! The square of the relative weighting of temperature and salinity changes
                         ! when extrapolating to match a target density [C2 S-2 ~> degC2 ppt-2]
   real    :: dT_dR      ! The ratio of temperature changes to density changes when
diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90
index 4e07b1d8ed..78ec0d9391 100644
--- a/src/parameterizations/vertical/MOM_kappa_shear.F90
+++ b/src/parameterizations/vertical/MOM_kappa_shear.F90
@@ -653,7 +653,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, dz, &
                 ! 1/dz_Int, as they have different uses.
     S2, &       ! The squared shear at an interface [T-2 ~> s-2].
     a1, &       ! a1 is the coupling between adjacent interfaces in the TKE,
-                ! velocity, and density equations [Z s-1 ~> m s-1] or [Z ~> m]
+                ! velocity, and density equations [Z ~> m]
     c1, &       ! c1 is used in the tridiagonal (and similar) solvers [nondim].
     k_src, &    ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
     kappa_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90
index 43462131ca..77de5d13cd 100644
--- a/src/parameterizations/vertical/MOM_opacity.F90
+++ b/src/parameterizations/vertical/MOM_opacity.F90
@@ -38,7 +38,7 @@ module MOM_opacity
       !< The maximum wavelength in each band of penetrating shortwave radiation [nm]
 
   real :: PenSW_flux_absorb !< A heat flux that is small enough to be completely absorbed in the next
-                        !! sufficiently thick layer [H degC T-1 ~> degC m s-1 or degC kg m-2 s-1].
+                        !! sufficiently thick layer [C H T-1 ~> degC m s-1 or degC kg m-2 s-1].
   real :: PenSW_absorb_Invlen !< The inverse of the thickness that is used to absorb the remaining
                         !! shortwave heat flux when it drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2].
   integer :: answer_date  !< The vintage of the order of arithmetic and expressions in the optics
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index bb159b2199..0dec7a40c0 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -900,7 +900,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
 
   real, dimension(SZI_(G)) :: &
     pres,      &  ! pressure at each interface [R L2 T-2 ~> Pa]
-    Temp_int,  &  ! temperature at each interface [C ~>degC]
+    Temp_int,  &  ! temperature at each interface [C ~> degC]
     Salin_int, &  ! salinity at each interface [S ~> ppt]
     drho_bot,  &  ! A density difference [R ~> kg m-3]
     h_amp,     &  ! The topographic roughness amplitude [Z ~> m].
diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90
index 2a9e7deeba..1e3bf258d8 100644
--- a/src/parameterizations/vertical/MOM_set_viscosity.F90
+++ b/src/parameterizations/vertical/MOM_set_viscosity.F90
@@ -287,7 +287,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
   real :: tmp              ! A temporary variable, sometimes in [Z ~> m]
   real :: tmp_val_m1_to_p1 ! A temporary variable [nondim]
   real :: curv_tol         ! Numerator of curvature cubed, used to estimate
-                           ! accuracy of a single L(:) Newton iteration [H5 ~> m3 or kg5 m-10]
+                           ! accuracy of a single L(:) Newton iteration [H5 ~> m5 or kg5 m-10]
   logical :: use_L0, do_one_L_iter    ! Control flags for L(:) Newton iteration
   logical :: use_BBL_EOS, do_i(SZIB_(G))
   integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
diff --git a/src/tracer/MOM_offline_aux.F90 b/src/tracer/MOM_offline_aux.F90
index d42355f245..7619cac2bd 100644
--- a/src/tracer/MOM_offline_aux.F90
+++ b/src/tracer/MOM_offline_aux.F90
@@ -777,9 +777,9 @@ subroutine update_offline_from_arrays(G, GV, nk_input, ridx_sum, mean_file, sum_
   real, dimension(:,:,:,:), allocatable,     intent(inout) :: hend_all  !< End of timestep layer thickness
                                                                         !! [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: temp      !< Temperature array [C ~> degC]
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: salt      !< Salinity array [ppt ~> S]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: salt      !< Salinity array [S ~> ppt]
   real, dimension(:,:,:,:), allocatable,     intent(inout) :: temp_all  !< Temperature array [C ~> degC]
-  real, dimension(:,:,:,:), allocatable,     intent(inout) :: salt_all  !< Salinity array [ppt ~> S]
+  real, dimension(:,:,:,:), allocatable,     intent(inout) :: salt_all  !< Salinity array [S ~> ppt]
 
   integer :: i, j, k, is, ie, js, je, nz
   real, parameter :: fill_value = 0.
diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90
index d887f5f3be..c089181c16 100644
--- a/src/tracer/MOM_tracer_Z_init.F90
+++ b/src/tracer/MOM_tracer_Z_init.F90
@@ -279,7 +279,7 @@ subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, n
                               intent(in)  :: tr_in !< The z-space array of tracer concentrations
                                                    !! that is read in [A]
   real, dimension(nk_data+1), intent(in)  :: z_edges !< The depths of the cell edges in the input z* data
-                                                   !! [Z ~> m or m]
+                                                   !! [Z ~> m] or [m]
   integer,                    intent(in)  :: nlay  !< The number of vertical layers in the target grid
   real, dimension(SZI_(G),SZJ_(G),nlay+1), &
                               intent(in)  :: e     !< The depths of the target layer interfaces [Z ~> m] or [m]
diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90
index 2244993447..fbc2b28a95 100644
--- a/src/tracer/dye_example.F90
+++ b/src/tracer/dye_example.F90
@@ -272,7 +272,7 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
                                               !! fluxes can be applied [H ~> m or kg m-2]
 
 ! Local variables
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-3]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
   real    :: z_bot    ! Height of the bottom of the layer relative to the sea surface [Z ~> m]
   real    :: z_center ! Height of the center of the layer relative to the sea surface [Z ~> m]
   integer :: i, j, k, is, ie, js, je, nz, m
diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90
index 36b157547f..20685d9711 100644
--- a/src/tracer/oil_tracer.F90
+++ b/src/tracer/oil_tracer.F90
@@ -363,7 +363,7 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
         CS%tr(i,j,k,m) = G%mask2dT(i,j)*max(1. - dt*CS%oil_decay_rate(m),0.)*CS%tr(i,j,k,m) ! Safest
       elseif (CS%oil_decay_rate(m)<0.) then
         decay_timescale = (12.0 * (3.0**(-(tv%T(i,j,k)-20.0*US%degC_to_C)/10.0*US%degC_to_C))) * &
-                          (86400.0*US%s_to_T) ! Timescale [s ~> T]
+                          (86400.0*US%s_to_T) ! Timescale [T ~> s]
         ldecay = 1. / decay_timescale ! Rate [T-1 ~> s-1]
         CS%tr(i,j,k,m) = G%mask2dT(i,j)*max(1. - dt*ldecay,0.)*CS%tr(i,j,k,m)
       endif
diff --git a/src/user/MOM_controlled_forcing.F90 b/src/user/MOM_controlled_forcing.F90
index 5be01bece4..d218b4ea80 100644
--- a/src/user/MOM_controlled_forcing.F90
+++ b/src/user/MOM_controlled_forcing.F90
@@ -270,8 +270,8 @@ subroutine apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, virt_heat, virt_prec
     ! Accumulate the average anomalies for this period.
     dt_wt = wt_per1 * dt
     CS%avg_time(m_mid) = CS%avg_time(m_mid) + dt_wt
-    ! These loops temporarily change the units of the CS%avg_ variables to [degC T ~> degC s]
-    ! or [ppt T ~> ppt s].
+    ! These loops temporarily change the units of the CS%avg_ variables to [C T ~> degC s]
+    ! or [S T ~> ppt s].
     do j=js,je ; do i=is,ie
       CS%avg_SST_anom(i,j,m_mid) = CS%avg_SST_anom(i,j,m_mid) + &
                                    dt_wt * G%mask2dT(i,j) * SST_anom(i,j)
diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90
index 723c8a0595..a548436329 100644
--- a/src/user/MOM_wave_interface.F90
+++ b/src/user/MOM_wave_interface.F90
@@ -198,7 +198,7 @@ module MOM_wave_interface
   real :: rho_air  !< A typical density of air at sea level, as used in wave calculations [R ~> kg m-3]
   real :: nu_air   !< The viscosity of air, as used in wave calculations [Z2 T-1 ~> m2 s-1]
   real :: SWH_from_u10sq !< A factor for converting the square of the 10 m wind speed to the
-                   !! significant wave height [Z T2 L-2 ~> s m-2]
+                   !! significant wave height [Z T2 L-2 ~> s2 m-1]
   real :: Charnock_min !< The minimum value of the Charnock coefficient, which relates the square of
                    !! the air friction velocity divided by the gravitational acceleration to the
                    !! wave roughness length [nondim]
@@ -598,7 +598,7 @@ subroutine set_LF17_wave_params(param_file, mdl, US, CS)
   call get_param(param_file, mdl, "WAVE_HEIGHT_SCALE_FACTOR", CS%SWH_from_u10sq, &
                  "A factor relating the square of the 10 m wind speed to the significant "//&
                  "wave height, with a default value based on the Pierson-Moskowitz spectrum.", &
-                 units="s m-2", default=0.0246, scale=US%m_to_Z*US%L_T_to_m_s**2)
+                 units="s2 m-1", default=0.0246, scale=US%m_to_Z*US%L_T_to_m_s**2)
   call get_param(param_file, mdl, "CHARNOCK_MIN", CS%Charnock_min, &
                  "The minimum value of the Charnock coefficient, which relates the square of "//&
                  "the air friction velocity divided by the gravitational acceleration to the "//&
diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90
index da18a7341c..9ff99b583f 100644
--- a/src/user/Rossby_front_2d_initialization.F90
+++ b/src/user/Rossby_front_2d_initialization.F90
@@ -122,7 +122,7 @@ subroutine Rossby_front_initialize_temperature_salinity(T, S, h, G, GV, US, &
 
   integer   :: i, j, k, is, ie, js, je, nz
   real      :: T_ref        ! Reference temperature within the surface layer [C ~> degC]
-  real      :: S_ref        ! Reference salinity within the surface layer [S ~> [ppt]
+  real      :: S_ref        ! Reference salinity within the surface layer [S ~> ppt]
   real      :: T_range      ! Range of temperatures over the vertical [C ~> degC]
   real      :: zc           ! Position of the middle of the cell [Z ~> m]
   real      :: zi           ! Bottom interface position relative to the sea surface [H ~> m or kg m-2]
diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90
index 60ce7b4a56..4a20f0e9b3 100644
--- a/src/user/tidal_bay_initialization.F90
+++ b/src/user/tidal_bay_initialization.F90
@@ -26,7 +26,7 @@ module tidal_bay_initialization
 !> Control structure for tidal bay open boundaries.
 type, public :: tidal_bay_OBC_CS ; private
   real :: tide_flow = 3.0e6  !< Maximum tidal flux with the tidal bay configuration [L2 Z T-1 ~> m3 s-1]
-  real :: tide_period        !< The period associated with the tidal bay configuration [T ~> s-1]
+  real :: tide_period        !< The period associated with the tidal bay configuration [T ~> s]
   real :: tide_ssh_amp       !< The magnitude of the sea surface height anomalies at the inflow
                              !! with the tidal bay configuration [Z ~> m]
 end type tidal_bay_OBC_CS