diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 3ee203210c..056b171ba8 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -189,12 +189,12 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, pbv, Wav
   real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: &
     PV, &       ! A diagnostic array of the potential vorticities [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
     RV          ! A diagnostic array of the relative vorticities [T-1 ~> s-1].
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) ::     CAuS !
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) ::     CAvS !
-  real :: fv1, fv2, fv3, fv4   ! (f+rv)*v [L T-2 ~> m s-2].
-  real :: fu1, fu2, fu3, fu4   ! -(f+rv)*u [L T-2 ~> m s-2].
-  real :: max_fv, max_fu       ! The maximum or minimum of the neighboring Coriolis
-  real :: min_fv, min_fu       ! accelerations [L T-2 ~> m s-2], i.e. max(min)_fu(v)q.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: CAuS ! Stokes contribution to CAu [L T-2 ~> m s-2]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: CAvS ! Stokes contribution to CAv [L T-2 ~> m s-2]
+  real :: fv1, fv2, fv3, fv4   ! (f+rv)*v at the 4 points surrounding a u points[L T-2 ~> m s-2]
+  real :: fu1, fu2, fu3, fu4   ! -(f+rv)*u at the 4 points surrounding a v point [L T-2 ~> m s-2]
+  real :: max_fv, max_fu       ! The maximum of the neighboring Coriolis accelerations [L T-2 ~> m s-2]
+  real :: min_fv, min_fu       ! The minimum of the neighboring Coriolis accelerations [L T-2 ~> m s-2]
 
   real, parameter :: C1_12 = 1.0 / 12.0 ! C1_12 = 1/12 [nondim]
   real, parameter :: C1_24 = 1.0 / 24.0 ! C1_24 = 1/24 [nondim]
diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90
index 871de51632..80630084b9 100644
--- a/src/core/MOM_checksum_packages.F90
+++ b/src/core/MOM_checksum_packages.F90
@@ -29,9 +29,9 @@ module MOM_checksum_packages
 
 !> A type for storing statistica about a variable
 type :: stats ; private
-  real :: minimum = 1.E34  !< The minimum value
-  real :: maximum = -1.E34 !< The maximum value
-  real :: average = 0.     !< The average value
+  real :: minimum = 1.E34  !< The minimum value [degC] or [ppt] or other units
+  real :: maximum = -1.E34 !< The maximum value [degC] or [ppt] or other units
+  real :: average = 0.     !< The average value [degC] or [ppt] or other units
 end type stats
 
 contains
diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90
index 6cffea5c75..e1fb3d3278 100644
--- a/src/core/MOM_density_integrals.F90
+++ b/src/core/MOM_density_integrals.F90
@@ -143,7 +143,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
   real :: p5(5)      ! Pressures at five quadrature points [R L2 T-2 ~> Pa]
   real :: r5(5)      ! Densities at five quadrature points [R ~> kg m-3]
   real :: rho_anom   ! The depth averaged density anomaly [R ~> kg m-3]
-  real, parameter :: C1_90 = 1.0/90.0  ! Rational constants.
+  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant [nondim]
   real :: GxRho      ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]
   real :: I_Rho      ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1]
   real :: dz         ! The layer thickness [Z ~> m]
@@ -784,7 +784,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
   real :: w_left, w_right  ! Left and right weights [nondim]
   real :: intz(5) ! The gravitational acceleration times the integrals of density
                   ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]
-  real, parameter :: C1_90 = 1.0/90.0  ! Rational constants.
+  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant [nondim]
   real :: GxRho ! The gravitational acceleration times density [R L2 Z-1 T-2 ~> kg m-2 s-2]
   real :: I_Rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1]
   real :: dz ! Layer thicknesses at tracer points [Z ~> m]
@@ -1175,7 +1175,7 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d
   real :: intp(5)    ! The integrals of specific volume with pressure at the
                      ! 5 sub-column locations [L2 T-2 ~> m2 s-2]
   logical :: do_massWeight ! Indicates whether to do mass weighting.
-  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant.
+  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant [nondim]
   integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, n, halo
 
   Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB
diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90
index c1eb749467..4e812b65d7 100644
--- a/src/core/MOM_porous_barriers.F90
+++ b/src/core/MOM_porous_barriers.F90
@@ -29,7 +29,7 @@ module MOM_porous_barriers
   type(diag_ctrl), pointer :: &
       diag => Null()                !< A structure to regulate diagnostic output timing
   logical :: debug                  !< If true, write verbose checksums for debugging purposes.
-  real    :: mask_depth             !< The depth shallower than which porous barrier is not applied.
+  real    :: mask_depth             !< The depth shallower than which porous barrier is not applied [Z ~> m]
   integer :: eta_interp             !< An integer indicating how the interface heights at the velocity
                                     !! points are calculated. Valid values are given by the parameters
                                     !! defined below: MAX, MIN, ARITHMETIC and HARMONIC.
diff --git a/src/ice_shelf/MOM_marine_ice.F90 b/src/ice_shelf/MOM_marine_ice.F90
index 2a1b6d799b..8635eb71b5 100644
--- a/src/ice_shelf/MOM_marine_ice.F90
+++ b/src/ice_shelf/MOM_marine_ice.F90
@@ -28,7 +28,7 @@ module MOM_marine_ice
 type, public :: marine_ice_CS ; private
   real :: kv_iceberg          !< The viscosity of the icebergs [L4 Z-2 T-1 ~> m2 s-1] (for ice rigidity)
   real :: berg_area_threshold !< Fraction of grid cell which iceberg must occupy
-                              !! so that fluxes below are set to zero. (0.5 is a
+                              !! so that fluxes below are set to zero [nondim]. (0.5 is a
                               !! good value to use.) Not applied for negative values.
   real :: latent_heat_fusion  !< Latent heat of fusion [Q ~> J kg-1]
   real :: density_iceberg     !< A typical density of icebergs [R ~> kg m-3] (for ice rigidity)
diff --git a/src/ice_shelf/user_shelf_init.F90 b/src/ice_shelf/user_shelf_init.F90
index c384ef7cee..4d1f263ca8 100644
--- a/src/ice_shelf/user_shelf_init.F90
+++ b/src/ice_shelf/user_shelf_init.F90
@@ -131,7 +131,9 @@ subroutine USER_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, C
   logical,                 intent(in)    :: new_sim !< If true, this the start of a new run.
 
 
-  real :: c1, edge_pos, slope_pos
+  real :: c1        ! The inverse of the range over which the shelf slopes [km-1]
+  real :: edge_pos  ! The time-evolving position the ice shelf edge [km]
+  real :: slope_pos ! The time-evolving position of the start of the ice shelf slope [km]
   integer :: i, j
 
   edge_pos = CS%pos_shelf_edge_0 + CS%shelf_speed*(time_type_to_real(Time) / 86400.0)
diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90
index e622b11805..8bea8fe6e9 100644
--- a/src/initialization/MOM_grid_initialize.F90
+++ b/src/initialization/MOM_grid_initialize.F90
@@ -365,7 +365,7 @@ subroutine set_grid_metrics_cartesian(G, param_file, US)
   real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) ! Axis labels [degrees_E] or [km] or [m]
   real :: dx_everywhere, dy_everywhere ! Grid spacings [L ~> m].
   real :: I_dx, I_dy                   ! Inverse grid spacings [L-1 ~> m-1].
-  real :: PI
+  real :: PI  ! The ratio of the circumference of a circle to its diameter [nondim]
   character(len=80) :: units_temp
   character(len=48) :: mdl  = "MOM_grid_init set_grid_metrics_cartesian"
 
@@ -922,7 +922,7 @@ end function  dL
 !! function fn takes the value fnval, also returning in ittmax the number of iterations of
 !! Newton's method that were used to polish the root.
 function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax)
-  real :: find_root !< The value of y where fn(y) = fnval that will be returned
+  real :: find_root !< The value of y where fn(y) = fnval that will be returned [radians]
   real,      external    :: fn    !< The external function whose root is being sought [gridpoints]
   real,      external    :: dy_df !< The inverse of the derivative of that function [radian gridpoint-1]
   type(GPS), intent(in)  :: GP    !< A structure of grid parameters
@@ -1128,12 +1128,12 @@ end function Int_dj_dy
 
 !> Extrapolates missing metric data into all the halo regions.
 subroutine extrapolate_metric(var, jh, missing)
-  real, dimension(:,:), intent(inout) :: var     !< The array in which to fill in halos [A]
+  real, dimension(:,:), intent(inout) :: var     !< The array in which to fill in halos [abitrary]
   integer,              intent(in)    :: jh      !< The size of the halos to be filled
-  real,       optional, intent(in)    :: missing !< The missing data fill value, 0 by default [A]
+  real,       optional, intent(in)    :: missing !< The missing data fill value, 0 by default [abitrary]
   ! Local variables
-  real :: badval
-  integer :: i,j
+  real :: badval ! A bad data value [abitrary]
+  integer :: i, j
 
   badval = 0.0 ; if (present(missing)) badval = missing
 
@@ -1162,8 +1162,8 @@ end subroutine extrapolate_metric
 !> This function implements Adcroft's rule for reciprocals, namely that
 !!   Adcroft_Inv(x) = 1/x for |x|>0 or 0 for x=0.
 function Adcroft_reciprocal(val) result(I_val)
-  real, intent(in) :: val  !< The value being inverted.
-  real :: I_val            !< The Adcroft reciprocal of val.
+  real, intent(in) :: val  !< The value being inverted [abitrary]
+  real :: I_val            !< The Adcroft reciprocal of val [abitrary-1]
 
   I_val = 0.0
   if (val /= 0.0) I_val = 1.0/val
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 1eec90f568..0e50ebb67f 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -885,10 +885,10 @@ subroutine initialize_thickness_uniform(h, depth_tot, G, GV, param_file, just_re
                                                       !! parameters without changing h.
   ! Local variables
   character(len=40)  :: mdl = "initialize_thickness_uniform" ! This subroutine's name.
-  real :: e0(SZK_(GV)+1)  ! The resting interface heights, in depth units, usually
+  real :: e0(SZK_(GV)+1)  ! The resting interface heights [Z ~> m], usually
                           ! negative because it is positive upward.
-  real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface
-                          ! positive upward, in depth units.
+  real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface,
+                          ! positive upward [Z ~> m].
   integer :: i, j, k, is, ie, js, je, nz
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
diff --git a/src/user/MOM_controlled_forcing.F90 b/src/user/MOM_controlled_forcing.F90
index 24d370e920..5be01bece4 100644
--- a/src/user/MOM_controlled_forcing.F90
+++ b/src/user/MOM_controlled_forcing.F90
@@ -397,7 +397,7 @@ end subroutine apply_ctrl_forcing
 
 !> This function maps rval into an integer in the range from 1 to num_period.
 function periodic_int(rval, num_period) result (m)
-  real,    intent(in) :: rval       !< Input for mapping.
+  real,    intent(in) :: rval       !< Input for mapping [nondim]
   integer, intent(in) :: num_period !< Maximum output.
   integer             :: m          !< Return value.
 
@@ -412,9 +412,9 @@ function periodic_int(rval, num_period) result (m)
 !> This function shifts rval by an integer multiple of num_period so that
 !! 0 <= val_out < num_period.
 function periodic_real(rval, num_period) result(val_out)
-  real,    intent(in) :: rval       !< Input to be shifted into valid range.
+  real,    intent(in) :: rval       !< Input to be shifted into valid range [nondim]
   integer, intent(in) :: num_period !< Maximum valid value.
-  real                :: val_out    !< Return value.
+  real                :: val_out    !< Return value [nondim]
   integer :: nshft
 
   if (rval < 0) then ; nshft = floor(abs(rval) / num_period) + 1
diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90
index f99ca27994..723c8a0595 100644
--- a/src/user/MOM_wave_interface.F90
+++ b/src/user/MOM_wave_interface.F90
@@ -1676,8 +1676,10 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS )
        pointer       :: CS !< Surface wave related control structure.
 
   ! Local variables
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: P_deltaStokes_L ! The stokes induced Pressure anomaly, layer averaged
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: P_deltaStokes_i ! The stokes induced Pressure anomaly at interfaces
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: P_deltaStokes_L ! The Stokes induced pressure anomaly,
+                                                              ! layer averaged [L2 T-2 ~> m2 s-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: P_deltaStokes_i ! The Stokes induced pressure anomaly
+                                                              ! at interfaces [L2 T-2 ~> m2 s-2]
   real :: P_Stokes_l, P_Stokes_r ! Stokes-induced pressure anomaly over layer (left/right of point) [L2 T-2 ~> m2 s-2]
   real :: P_Stokes_l0, P_Stokes_r0 ! Stokes-induced pressure anomaly at interface
                                    ! (left/right of point) [L2 T-2 ~> m2 s-2]
@@ -1690,11 +1692,12 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS )
   real :: zi_l(SZK_(G)+1), zi_r(SZK_(G)+1)   ! The height of the edges of the cells (left/right of point) [Z ~> m].
   real :: idz_l(SZK_(G)), idz_r(SZK_(G)) ! The inverse thickness of the cells (left/right of point) [Z-1 ~> m-1]
   real :: h_l, h_r   ! The thickness of the cell (left/right of point) [Z ~> m].
-  real :: dexp2kzL,dexp4kzL,dexp2kzR,dexp4kzR ! Analytical evaluation of multi-exponential decay contribution
-                                              ! to Stokes pressure anomalies.
-  real :: TwoK, FourK, iTwoK, iFourK ! Wavenumber multipliers/inverses
+  real :: dexp2kzL, dexp4kzL, dexp2kzR, dexp4kzR ! Analytical evaluation of multi-exponential decay
+                                              ! contribution to Stokes pressure anomalies [nondim].
+  real :: TwoK, FourK   ! Wavenumbers multiplied by a factor [Z-1 ~> m-1]
+  real :: iTwoK, iFourK ! Inverses of wavenumbers [Z ~> m]
 
-  integer :: i,j,k,l
+  integer :: i, j, k, l
 
   !---------------------------------------------------------------
   ! Compute the Stokes contribution to the pressure gradient force
diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90
index 06b3ed43d6..62b55bb0a1 100644
--- a/src/user/Phillips_initialization.F90
+++ b/src/user/Phillips_initialization.F90
@@ -233,16 +233,16 @@ subroutine Phillips_initialize_sponges(G, GV, US, tv, param_file, CSp, h)
   real, intent(in), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< Thickness field [H ~> m or kg m-2].
 
   ! Local variables
-  real :: eta0(SZK_(GV)+1)  ! The 1-d nominal positions of the interfaces.
+  real :: eta0(SZK_(GV)+1)  ! The 1-d nominal positions of the interfaces [Z ~> m]
   real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m].
-  real :: temp(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for other variables.
+  real :: temp(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for other variables [various]
   real :: Idamp(SZI_(G),SZJ_(G))    ! The sponge damping rate [T-1 ~> s-1]
   real :: eta_im(SZJ_(G),SZK_(GV)+1) ! A temporary array for zonal-mean eta [Z ~> m].
   real :: Idamp_im(SZJ_(G))         ! The inverse zonal-mean damping rate [T-1 ~> s-1].
   real :: damp_rate    ! The inverse zonal-mean damping rate [T-1 ~> s-1].
-  real :: jet_width    ! The width of the zonal mean jet, in km.
+  real :: jet_width    ! The width of the zonal mean jet [km].
   real :: jet_height   ! The interface height scale associated with the zonal-mean jet [Z ~> m].
-  real :: y_2          ! The y-position relative to the channel center, in km.
+  real :: y_2          ! The y-position relative to the channel center [km].
   real :: half_strat   ! The fractional depth where the straficiation is centered [nondim].
   real :: half_depth   ! The depth where the stratification is centered [Z ~> m].
   real :: pi              ! The ratio of the circumference of a circle to its diameter [nondim]
@@ -309,8 +309,8 @@ end subroutine Phillips_initialize_sponges
 
 !> sech calculates the hyperbolic secant.
 function sech(x)
-  real, intent(in) :: x    !< Input value.
-  real             :: sech !< Result.
+  real, intent(in) :: x    !< Input value [nondim].
+  real             :: sech !< Result [nondim].
 
   ! This is here to prevent overflows or underflows.
   if (abs(x) > 228.) then
@@ -330,9 +330,14 @@ subroutine Phillips_initialize_topography(D, G, param_file, max_depth, US)
   type(unit_scale_type),           intent(in)  :: US !< A dimensional unit scaling type
 
   ! Local variables
-  real :: PI, Htop, Wtop, Ltop, offset, dist
-  real :: x1, x2, x3, x4, y1, y2
-  integer :: i,j,is,ie,js,je
+  real :: PI       ! The ratio of the circumference of a circle to its diameter [nondim]
+  real :: Htop     ! The maximum height of the topography above max_depth [Z ~> m]
+  real :: Wtop     ! meridional width of topographic features [km]
+  real :: Ltop     ! zonal width of topographic features [km]
+  real :: offset   ! meridional offset from the center of topographic features [km]
+  real :: dist     ! zonal width of topographic features [km]
+  real :: x1, x2, x3, x4, y1, y2 ! Various positions in the domain [km]
+  integer :: i, j, is, ie, js, je
   character(len=40)  :: mdl = "Phillips_initialize_topography" ! This subroutine's name.
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
@@ -349,8 +354,8 @@ subroutine Phillips_initialize_topography(D, G, param_file, max_depth, US)
   dist = 0.333*G%len_lon   ! distance between drake and mount
                            ! should be longer than Ltop/2
 
-  y1=G%south_lat+0.5*G%len_lat+offset-0.5*Wtop; y2=y1+Wtop
-  x1=G%west_lon+0.1*G%len_lon; x2=x1+Ltop; x3=x1+dist; x4=x3+3.0/2.0*Ltop
+  y1 = G%south_lat+0.5*G%len_lat+offset-0.5*Wtop ; y2 = y1+Wtop
+  x1 = G%west_lon+0.1*G%len_lon ; x2 = x1+Ltop ; x3 = x1+dist ; x4 = x3+3.0/2.0*Ltop
 
   do j=js,je ; do i=is,ie
     D(i,j)=0.0
diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90
index f681231694..8df8f90e3d 100644
--- a/src/user/SCM_CVMix_tests.F90
+++ b/src/user/SCM_CVMix_tests.F90
@@ -197,7 +197,7 @@ subroutine SCM_CVMix_tests_wind_forcing(sfc_state, forces, day, G, US, CS)
   ! Local variables
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
-  real    :: mag_tau
+  real    :: mag_tau  ! The magnitude of the wind stress [R L Z T-2 ~> Pa]
   ! Bounds for loops and memory allocation
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
@@ -231,7 +231,7 @@ subroutine SCM_CVMix_tests_buoyancy_forcing(sfc_state, fluxes, day, G, US, CS)
   ! Local variables
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
-  real :: PI
+  real :: PI  ! The ratio of the circumference of a circle to its diameter [nondim]
 
   PI = 4.0*atan(1.0)
 
diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90
index a9344a6a30..3920b52729 100644
--- a/src/user/benchmark_initialization.F90
+++ b/src/user/benchmark_initialization.F90
@@ -41,7 +41,7 @@ subroutine benchmark_initialize_topography(D, G, param_file, max_depth, US)
 
   ! Local variables
   real :: min_depth ! The minimum basin depth [Z ~> m]
-  real :: PI        ! 3.1415926... calculated as 4*atan(1)
+  real :: PI        ! 3.1415926... calculated as 4*atan(1) [nondim]
   real :: D0        ! A constant to make the maximum basin depth MAXIMUM_DEPTH [Z ~> m]
   real :: x         ! Longitude relative to the domain edge, normalized by its extent [nondim]
   real :: y         ! Latitude relative to the domain edge, normalized by its extent [nondim]
@@ -113,7 +113,7 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e
     drho_dT, &      ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1].
     drho_dS         ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
   real :: pres(SZK_(GV))  ! Reference pressure [R L2 T-2 ~> Pa].
-  real :: a_exp     ! The fraction of the overall stratification that is exponential.
+  real :: a_exp     ! The fraction of the overall stratification that is exponential [nondim]
   real :: I_ts, I_md ! Inverse lengthscales [Z-1 ~> m-1].
   real :: T_frac    ! A ratio of the interface temperature to the range
                     ! between SST and the bottom temperature [nondim].
@@ -121,7 +121,7 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e
                     ! interface temperature for a given z [nondim]
   real :: derr_dz   ! The derivative of the normalized error between the profile's
                     ! temperature and the interface temperature with z [Z-1 ~> m-1]
-  real :: pi        ! 3.1415926... calculated as 4*atan(1)
+  real :: pi        ! 3.1415926... calculated as 4*atan(1) [nondim]
   real :: z         ! A work variable for the interface position [Z ~> m]
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
@@ -246,7 +246,7 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, &
   real :: drho_dT(SZK_(GV))  ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
   real :: drho_dS(SZK_(GV))  ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
   real :: rho_guess(SZK_(GV)) ! Potential density at T0 & S0 [R ~> kg m-3]
-  real :: PI                 ! 3.1415926... calculated as 4*atan(1)
+  real :: PI                 ! 3.1415926... calculated as 4*atan(1) [nondim]
   real :: SST                !  The initial sea surface temperature [C ~> degC]
   character(len=40)  :: mdl = "benchmark_init_temperature_salinity" ! This subroutine's name.
   integer :: i, j, k, k1, is, ie, js, je, nz, itt
diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90
index 07fc539979..63c5c8a0d4 100644
--- a/src/user/circle_obcs_initialization.F90
+++ b/src/user/circle_obcs_initialization.F90
@@ -44,7 +44,12 @@ subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, jus
   real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface
                            ! positive upward, in depth units [Z ~> m].
   real :: IC_amp           ! The amplitude of the initial height displacement [H ~> m or kg m-2].
-  real :: diskrad, rad, lonC, latC, xOffset
+  real :: diskrad          ! Radius of the elevated disk [km] or [degrees] or [m]
+  real :: rad              ! Distance from the center of the elevated disk [km] or [degrees] or [m]
+  real :: lonC             ! The x-position of a point [km] or [degrees] or [m]
+  real :: latC             ! The y-position of a point [km] or [degrees] or [m]
+  real :: xOffset          ! The x-offset of the elevated disc center relative to the domain
+                           ! center [km] or [degrees] or [m]
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "circle_obcs_initialization"   ! This module's name.
diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90
index 5e0cb65007..81aa4c2b3b 100644
--- a/src/user/dense_water_initialization.F90
+++ b/src/user/dense_water_initialization.F90
@@ -35,9 +35,9 @@ module dense_water_initialization
 subroutine dense_water_initialize_topography(D, G, param_file, max_depth)
   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 the units of depth_max
+                           intent(out) :: D !< Ocean bottom depth [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file !< Parameter file structure
-  real,                    intent(in)  :: max_depth !< Maximum ocean depth in arbitrary units
+  real,                    intent(in)  :: max_depth !< Maximum ocean depth [Z ~> m]
 
   ! Local variables
   real, dimension(5) :: domain_params ! nondimensional widths of all domain sections [nondim]
diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90
index 0b8f59a6e8..4ac5ab3bf9 100644
--- a/src/user/dumbbell_surface_forcing.F90
+++ b/src/user/dumbbell_surface_forcing.F90
@@ -34,7 +34,7 @@ module dumbbell_surface_forcing
                              !! to the reservoirs
   real :: slp_period         !< Period of sinusoidal pressure wave [days]
   real, dimension(:,:), allocatable :: &
-    forcing_mask             !< A mask regulating where forcing occurs
+    forcing_mask             !< A mask regulating where forcing occurs [nondim]
   real, dimension(:,:), allocatable :: &
     S_restore                !< The surface salinity field toward which to restore [S ~> ppt].
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate the
diff --git a/src/user/dyed_channel_initialization.F90 b/src/user/dyed_channel_initialization.F90
index 411ab6ef98..aed7142fad 100644
--- a/src/user/dyed_channel_initialization.F90
+++ b/src/user/dyed_channel_initialization.F90
@@ -94,7 +94,7 @@ subroutine dyed_channel_set_OBC_tracer_data(OBC, G, GV, param_file, tr_Reg)
   character(len=40)  :: mdl = "dyed_channel_set_OBC_tracer_data" ! This subroutine's name.
   character(len=80)  :: name, longname
   integer :: m, n
-  real :: dye
+  real :: dye  ! Inflow dye concentrations [arbitrary]
   type(tracer_type), pointer      :: tr_ptr => NULL()
 
   if (.not.associated(OBC)) call MOM_error(FATAL, 'dyed_channel_initialization.F90: '// &
@@ -142,7 +142,7 @@ subroutine dyed_channel_update_flow(OBC, CS, G, GV, US, Time)
   type(time_type),            intent(in) :: Time !< model time.
   ! Local variables
   real :: flow      ! The OBC velocity [L T-1 ~> m s-1]
-  real :: PI        ! 3.1415926535...
+  real :: PI        ! 3.1415926535... [nondim]
   real :: time_sec  ! The elapsed time since the start of the calendar [T ~> s]
   integer :: i, j, k, l, isd, ied, jsd, jed
   integer :: IsdB, IedB, JsdB, JedB
diff --git a/src/user/dyed_obcs_initialization.F90 b/src/user/dyed_obcs_initialization.F90
index c5efef4905..6248efab2f 100644
--- a/src/user/dyed_obcs_initialization.F90
+++ b/src/user/dyed_obcs_initialization.F90
@@ -41,7 +41,7 @@ subroutine dyed_obcs_set_OBC_data(OBC, G, GV, param_file, tr_Reg)
   character(len=80)  :: name, longname
   integer :: is, ie, js, je, isd, ied, jsd, jed, m, n, nz
   integer :: IsdB, IedB, JsdB, JedB
-  real :: dye
+  real :: dye ! Inflow dye concentration [arbitrary]
   type(tracer_type), pointer      :: tr_ptr => NULL()
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
diff --git a/src/user/external_gwave_initialization.F90 b/src/user/external_gwave_initialization.F90
index 554440dbcb..63cc89342a 100644
--- a/src/user/external_gwave_initialization.F90
+++ b/src/user/external_gwave_initialization.F90
@@ -44,7 +44,8 @@ subroutine external_gwave_initialize_thickness(h, G, GV, US, param_file, just_re
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   integer :: i, j, k, is, ie, js, je, nz
-  real :: PI, Xnondim
+  real :: PI       ! The ratio of the circumference of a circle to its diameter [nondim]
+  real :: Xnondim  ! A normalized x position [nondim]
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
diff --git a/src/user/lock_exchange_initialization.F90 b/src/user/lock_exchange_initialization.F90
index a3418e6482..3b41237c36 100644
--- a/src/user/lock_exchange_initialization.F90
+++ b/src/user/lock_exchange_initialization.F90
@@ -36,8 +36,8 @@ subroutine lock_exchange_initialize_thickness(h, G, GV, US, param_file, just_rea
 
   real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface
                            ! positive upward [Z ~> m].
-  real :: front_displacement ! Vertical displacement acrodd front
-  real :: thermocline_thickness ! Thickness of stratified region
+  real :: front_displacement ! Vertical displacement across front [Z ~> m]
+  real :: thermocline_thickness ! Thickness of stratified region [Z ~> m]
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "lock_exchange_initialize_thickness" ! This subroutine's name.
diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90
index 5d7f1b7e97..df46a142f1 100644
--- a/src/user/shelfwave_initialization.F90
+++ b/src/user/shelfwave_initialization.F90
@@ -138,7 +138,7 @@ subroutine shelfwave_set_OBC_data(OBC, CS, G, GV, US, h, Time)
   type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
   type(verticalGrid_type), intent(in) :: GV   !< The ocean's vertical grid structure
   type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
   type(time_type),         intent(in) :: Time !< model time.
 
   ! The following variables are used to set up the transport in the shelfwave example.
diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90
index f125da0a25..c12d34a721 100644
--- a/src/user/user_change_diffusivity.F90
+++ b/src/user/user_change_diffusivity.F90
@@ -151,7 +151,7 @@ end subroutine user_change_diff
 
 !> This subroutine checks whether the 4 values of range are in ascending order.
 function range_OK(range) result(OK)
-  real, dimension(4), intent(in) :: range  !< Four values to check.
+  real, dimension(4), intent(in) :: range  !< Four values to check [arbitrary]
   logical                        :: OK     !< Return value.
 
   OK = ((range(1) <= range(2)) .and. (range(2) <= range(3)) .and. &
@@ -169,7 +169,7 @@ function val_weights(val, range) result(ans)
   real, dimension(4), intent(in) :: range  !< Range over which the answer is non-zero [arbitrary units].
   real                           :: ans    !< Return value [nondim].
   ! Local variables
-  real :: x   ! A nondimensional number between 0 and 1.
+  real :: x   ! A nondimensional number between 0 and 1 [nondim].
 
   ans = 0.0
   if ((val > range(1)) .and. (val < range(4))) then
diff --git a/src/user/user_revise_forcing.F90 b/src/user/user_revise_forcing.F90
index eb9694a091..ce767d7479 100644
--- a/src/user/user_revise_forcing.F90
+++ b/src/user/user_revise_forcing.F90
@@ -21,7 +21,7 @@ module user_revise_forcing
 
 !> Control structure for user_revise_forcing
 type, public :: user_revise_forcing_CS ; private
-  real :: cdrag  !< The quadratic bottom drag coefficient.
+  real :: cdrag  !< The quadratic bottom drag coefficient [nondim]
 end type user_revise_forcing_CS
 
 contains