Document units of 62 real variables

  Added standard-format unit descriptions for 62 real variables in comments
scattered across 17 modules in the core, tracer, and both parameterizations
directories.  Another 6 comments were added noting dimensionally inconsistent
expressions in the recently added subroutine vertFPmix.  Only comments are
changed and all answers are bitwise identical.

src/core/MOM_barotropic.F90

@@ -409,7 +409,7 @@ module MOM_barotropic
 !>@}
 
 !> A negligible parameter which avoids division by zero, but is too small to
-!! modify physical values.
+!! modify physical values [nondim].
 real, parameter :: subroundoff = 1e-30
 
 contains
@@ -662,7 +662,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     time_bt_start, &  ! The starting time of the barotropic steps.
     time_step_end, &  ! The end time of a barotropic step.
     time_end_in       ! The end time for diagnostics when this routine started.
-  real :: time_int_in ! The diagnostics' time interval when this routine started.
+  real :: time_int_in ! The diagnostics' time interval when this routine started [s]
   real :: Htot_avg    ! The average total thickness of the tracer columns adjacent to a
                       ! velocity point [H ~> m or kg m-2]
   logical :: do_hifreq_output  ! If true, output occurs every barotropic step.
@@ -3972,7 +3972,7 @@ subroutine set_local_BT_cont_types(BT_cont, BTCL_u, BTCL_v, G, US, MS, BT_Domain
     vBT_NN, vBT_SS, &  ! Meridional velocities at which the form of the fit changes [L T-1 ~> m s-1]
     FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS ! Meridional face areas [H L ~> m2 or kg m-1]
   real :: dt ! The baroclinic timestep [T ~> s] or 1.0 [nondim]
-  real, parameter :: C1_3 = 1.0/3.0
+  real, parameter :: C1_3 = 1.0/3.0  ! [nondim]
   integer :: i, j, is, ie, js, je, hs
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
@@ -4107,7 +4107,7 @@ subroutine adjust_local_BT_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, &
 
   ! Local variables
   real :: dt ! The baroclinic timestep [T ~> s] or 1.0 [nondim]
-  real, parameter :: C1_3 = 1.0/3.0
+  real, parameter :: C1_3 = 1.0/3.0  ! [nondim]
   integer :: i, j, is, ie, js, je, hs
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec

src/core/MOM_dynamics_split_RK2.F90

@@ -387,7 +387,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
     h_av    ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
 
-  real, dimension(SZI_(G),SZJ_(G)) :: hbl           ! Boundary layer depth from Cvmix
+  real, dimension(SZI_(G),SZJ_(G)) :: hbl       ! Boundary layer depth from Cvmix [H ~> m or kg m-2]
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
   real :: Idt_bc    ! Inverse of the baroclinic timestep [T-1 ~> s-1]
   logical :: dyn_p_surf
@@ -1055,8 +1055,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   endif
 
   if (CS%fpmix) then
-    if (CS%id_uold > 0) call post_data(CS%id_uold , uold, CS%diag)
-    if (CS%id_vold > 0) call post_data(CS%id_vold , vold, CS%diag)
+    if (CS%id_uold > 0) call post_data(CS%id_uold, uold, CS%diag)
+    if (CS%id_vold > 0) call post_data(CS%id_vold, vold, CS%diag)
   endif
 
   ! The time-averaged free surface height has already been set by the last call to btstep.
@@ -1072,8 +1072,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   if (CS%id_CAv > 0) call post_data(CS%id_CAv, CS%CAv, CS%diag)
 
   ! Here the thickness fluxes are offered for time averaging.
-  if (CS%id_uh         > 0) call post_data(CS%id_uh , uh,                   CS%diag)
-  if (CS%id_vh         > 0) call post_data(CS%id_vh , vh,                   CS%diag)
+  if (CS%id_uh         > 0) call post_data(CS%id_uh,  uh,                   CS%diag)
+  if (CS%id_vh         > 0) call post_data(CS%id_vh,  vh,                   CS%diag)
   if (CS%id_uav        > 0) call post_data(CS%id_uav, u_av,                 CS%diag)
   if (CS%id_vav        > 0) call post_data(CS%id_vav, v_av,                 CS%diag)
   if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag)
@@ -1301,7 +1301,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
                                     intent(inout) :: u          !< zonal velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                                     intent(inout) :: v          !< merid velocity [L T-1 ~> m s-1]
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) , &
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  &
                                     intent(inout) :: h          !< layer thickness [H ~> m or kg m-2]
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                             target, intent(inout) :: uh    !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

src/core/MOM_forcing_type.F90

@@ -2169,7 +2169,7 @@ subroutine forcing_accumulate(flux_tmp, forces, fluxes, G, wt2)
   type(forcing),         intent(inout) :: fluxes !< A structure containing time-averaged
                                                  !! thermodynamic forcing fields
   type(ocean_grid_type), intent(inout) :: G      !< The ocean's grid structure
-  real,                  intent(out)   :: wt2    !< The relative weight of the new fluxes
+  real,                  intent(out)   :: wt2    !< The relative weight of the new fluxes [nondim]
 
   ! This subroutine copies mechancal forcing from flux_tmp to fluxes and
   ! stores the time-weighted averages of the various buoyancy fluxes in fluxes,
@@ -2187,7 +2187,7 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces)
   type(forcing),             intent(inout) :: fluxes !< A structure containing time-averaged
                                                      !! thermodynamic forcing fields
   type(ocean_grid_type),     intent(inout) :: G      !< The ocean's grid structure
-  real,                      intent(out)   :: wt2    !< The relative weight of the new fluxes
+  real,                      intent(out)   :: wt2    !< The relative weight of the new fluxes [nondim]
   type(mech_forcing), optional, intent(in) :: forces !< A structure with the driving mechanical forces
 
   ! This subroutine copies mechanical forcing from flux_tmp to fluxes and

src/core/MOM_grid.F90

@@ -136,14 +136,18 @@ module MOM_grid
     IareaBu      !< IareaBu = 1/areaBu [L-2 ~> m-2].
 
   real, pointer, dimension(:) :: &
-    gridLatT => NULL(), & !< The latitude of T points for the purpose of labeling the output axes.
+    gridLatT => NULL(), & !< The latitude of T points for the purpose of labeling the output axes,
+                          !! often in units of [degrees_N] or [km] or [m] or [gridpoints].
                           !! On many grids this is the same as geoLatT.
-    gridLatB => NULL()    !< The latitude of B points for the purpose of labeling the output axes.
+    gridLatB => NULL()    !< The latitude of B points for the purpose of labeling the output axes,
+                          !! often in units of [degrees_N] or [km] or [m] or [gridpoints].
                           !! On many grids this is the same as geoLatBu.
   real, pointer, dimension(:) :: &
-    gridLonT => NULL(), & !< The longitude of T points for the purpose of labeling the output axes.
+    gridLonT => NULL(), & !< The longitude of T points for the purpose of labeling the output axes,
+                          !! often in units of [degrees_E] or [km] or [m] or [gridpoints].
                           !! On many grids this is the same as geoLonT.
-    gridLonB => NULL()    !< The longitude of B points for the purpose of labeling the output axes.
+    gridLonB => NULL()    !< The longitude of B points for the purpose of labeling the output axes,
+                          !! often in units of [degrees_E] or [km] or [m] or [gridpoints].
                           !! On many grids this is the same as geoLonBu.
   character(len=40) :: &
     ! Except on a Cartesian grid, these are usually some variant of "degrees".
@@ -187,8 +191,8 @@ module MOM_grid
   ! initialization routines (but not all)
   real :: south_lat     !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m]
   real :: west_lon      !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m]
-  real :: len_lat       !< The latitudinal (or y-coord) extent of physical domain
-  real :: len_lon       !< The longitudinal (or x-coord) extent of physical domain
+  real :: len_lat       !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m]
+  real :: len_lon       !< The longitudinal (or x-coord) extent of physical domain [degrees_E] or [km] or [m]
   real :: Rad_Earth     !< The radius of the planet [m]
   real :: Rad_Earth_L   !< The radius of the planet in rescaled units [L ~> m]
   real :: max_depth     !< The maximum depth of the ocean in depth units [Z ~> m]

src/core/MOM_interface_heights.F90

@@ -66,7 +66,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)
   real :: p(SZI_(G),SZJ_(G),SZK_(GV)+1)   ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
   real :: dz_geo(SZI_(G),SZJ_(G),SZK_(GV)) ! The change in geopotential height
                                            ! across a layer [L2 T-2 ~> m2 s-2].
-  real :: dilate(SZI_(G))                 ! non-dimensional dilation factor
+  real :: dilate(SZI_(G))                 ! A non-dimensional dilation factor [nondim]
   real :: htot(SZI_(G))                   ! total thickness [H ~> m or kg m-2]
   real :: I_gEarth          ! The inverse of the gravitational acceleration times the
                             ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]

src/parameterizations/lateral/MOM_mixed_layer_restrat.F90

@@ -837,7 +837,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
   real :: muza            ! mu(z) at top of the layer [nondim]
   real :: dh              ! Portion of the layer thickness that is in the mixed layer [H ~> m or kg m-2]
   real :: res_scaling_fac ! The resolution-dependent scaling factor [nondim]
-  real, parameter :: two_thirds = 2./3.
+  real, parameter :: two_thirds = 2./3.  ! [nondim]
   logical :: line_is_empty, keep_going
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
@@ -1156,7 +1156,7 @@ real elemental function rmean2ts(signal, filtered, tau_growing, tau_decaying, dt
   real, intent(in) :: tau_decaying ! Time scale for decaying signal [T ~> s]
   real, intent(in) :: dt           ! Time step [T ~> s]
   ! Local variables
-  real :: afac, bfac ! Non-dimensional weights
+  real :: afac, bfac ! Non-dimensional fractional weights [nondim]
   real :: rt ! Reciprocal time scale [T-1 ~> s-1]
 
   if (signal>=filtered) then

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -1671,11 +1671,11 @@ subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: int_slope_u !< Ratio that determine how much of
                                                                       !! the isopycnal slopes are taken directly from
                                                                       !! the interface slopes without consideration
-                                                                      !! of density gradients.
+                                                                      !! of density gradients [nondim].
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: int_slope_v !< Ratio that determine how much of
                                                                       !! the isopycnal slopes are taken directly from
                                                                       !! the interface slopes without consideration
-                                                                      !! of density gradients.
+                                                                      !! of density gradients [nondim].
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
     de_top     ! The distances between the top of a layer and the top of the

src/parameterizations/lateral/MOM_tidal_forcing.F90

@@ -256,7 +256,7 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS)
   character(len=40)  :: mdl = "MOM_tidal_forcing" ! This module's name.
   character(len=128) :: mesg
   character(len=200) :: tidal_input_files(4*MAX_CONSTITUENTS)
-  real :: tide_sal_scalar_value
+  real :: tide_sal_scalar_value ! The constant of proportionality with the scalar approximation to SAL [nondim]
   integer :: i, j, c, is, ie, js, je, isd, ied, jsd, jed, nc
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -108,7 +108,7 @@ module MOM_energetic_PBL
 
   !/ mstar_scheme == 0
   real    :: fixed_mstar     !< Mstar is the ratio of the friction velocity cubed to the TKE available to
-                             !! drive entrainment, nondimensional. This quantity is the vertically
+                             !! drive entrainment [nondim]. This quantity is the vertically
                              !! integrated shear production minus the vertically integrated
                              !! dissipation of TKE produced by shear.  This value is used if the option
                              !! for using a fixed mstar is used.

src/parameterizations/vertical/MOM_opacity.F90

@@ -422,7 +422,7 @@ function opacity_morel(chl_data)
 !> This sets the penetrating shortwave fraction according to the scheme proposed by
 !! Morel and Antoine (1994).
 function SW_pen_frac_morel(chl_data)
-  real, intent(in)  :: chl_data !< The chlorophyll-A concentration in mg m-3.
+  real, intent(in)  :: chl_data !< The chlorophyll-A concentration [mg m-3]
   real :: SW_pen_frac_morel     !< The returned penetrating shortwave fraction [nondim]
 
   !   The following are coefficients for the optical model taken from Morel and
@@ -608,7 +608,7 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l
   real :: SW_trans          ! fraction of shortwave radiation that is not
                             ! absorbed in a layer [nondim]
   real :: unabsorbed        ! fraction of the shortwave radiation that
-                            ! is not absorbed because the layers are too thin
+                            ! is not absorbed because the layers are too thin [nondim]
   real :: Ih_limit          ! inverse of the total depth at which the
                             ! surface fluxes start to be limited [H-1 ~> m-1 or m2 kg-1]
   real :: h_min_heat        ! minimum thickness layer that should get heated [H ~> m or kg m-2]

src/parameterizations/vertical/MOM_vert_friction.F90

@@ -234,11 +234,15 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2w_u !< angle between mtm flux and wind at u-pts [rad]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2w_v !< angle between mtm flux and wind at v-pts [rad]
 
-  real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp, omega_tmp !< constants and dummy variables
-  real :: du, dv, depth, sigma, Wind_x, Wind_y          !< intermediate variables
-  real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables
-  real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables
-  real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles
+  real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp  !< constants and dummy variables [nondim]
+  real :: omega_tmp        !< A dummy angle [radians]
+  real :: du, dv           !< Velocity increments [L T-1 ~> m s-1]
+  real :: depth            !< Cumulative layer thicknesses [H ~> m or kg m=2]
+  real :: sigma            !< Fractional depth in the mixed layer [nondim]
+  real :: Wind_x, Wind_y   !< intermediate wind stress componenents [L2 T-2 ~> m2 s-2]
+  real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables [L2 T-2 ~> m2 s-2]
+  real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables [L2 T-2 ~> m2 s-2]
+  real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles [radians]
   integer :: kblmin, kbld, kp1, k, nz !< vertical indices
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq ! horizontal indices
 
@@ -321,6 +325,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB
   enddo
 
   if (CS%debug) then
+    !### These checksum calls are missing necessary dimensionally scaling factors.
     call uvchksum("surface tau[xy]_[uv] ", taux_u, tauy_v, G%HI, haloshift=1, scalar_pair=.true.)
     call uvchksum("ustar2", ustar2_u, ustar2_v, G%HI, haloshift=0, scalar_pair=.true.)
     call uvchksum(" hbl", hbl_u ,   hbl_v , G%HI, haloshift=0, scalar_pair=.true.)
@@ -427,6 +432,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB
         kbld  = min( (kbl_u(I,j)) , (nz-2) )
         if ( tau_u(I,j,kbld+2) > tau_u(I,j,kbld+1) ) kbld = kbld + 1
 
+        !### This expression is dimensionally inconsistent.
         tauh  =  tau_u(I,j,kbld+1) + GV%H_subroundoff
         ! surface boundary conditions
         depth   = 0.
@@ -437,6 +443,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB
 
           ! linear stress mag
           tau_MAG   = (ustar2_u(I,j) * (1.-sigma) )  + (tauh * sigma )
+          !### The following expressions are dimensionally inconsistent.
           cos_tmp   = tauxDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff)
           sin_tmp   = tauyDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff)
 
@@ -457,6 +464,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB
           tauNLdn  = tauNL_X
 
           ! nonlocal increment and update to uold
+          !### The following expression is dimensionally inconsistent and missing parentheses.
           du = (tauNLup - tauNLdn) * (dt/CS%h_u(I,j,k) + GV%H_subroundoff)
           ui(I,j,k)    = uold(I,j,k)  + du
           uold(I,j,k)  = du
@@ -496,6 +504,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB
 
           ! linear stress
           tau_MAG   = (ustar2_v(i,J) * (1.-sigma))  + (tauh * sigma)
+          !### The following expressions are dimensionally inconsistent.
           cos_tmp   = tauxDG_v(i,J,k+1) / (tau_v(i,J,k+1)  + GV%H_subroundoff)
           sin_tmp   = tauyDG_v(i,J,k+1) / (tau_v(i,J,k+1)  + GV%H_subroundoff)
 
@@ -514,6 +523,8 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB
           tauNL_X  = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp)
           tauNL_Y  = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp)
           tauNLdn  = tauNL_Y
+          !### The following expression is dimensionally inconsistent, [L T-1] vs. [L2 H-1 T-1] on the right,
+          !    and it is inconsistent with the counterpart expression for du.
           dv            = (tauNLup - tauNLdn) * (dt/(CS%h_v(i,J,k)) )
           vi(i,J,k)    = vold(i,J,k) + dv
           vold(i,J,k)  = dv
@@ -553,10 +564,10 @@ end subroutine vertFPmix
 
 !> Returns the empirical shape-function given sigma.
 real function G_sig(sigma)
-  real , intent(in) :: sigma   !< non-dimensional normalized boundary layer depth [m]
+  real , intent(in) :: sigma   !< non-dimensional normalized boundary layer depth [nondim]
 
   ! local variables
-  real :: p1, c2, c3  !< parameters used to fit and match empirycal shape-functions.
+  real :: p1, c2, c3  !< parameters used to fit and match empirical shape-functions [nondim]
 
   ! parabola
   p1 = 0.287
@@ -2634,7 +2645,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
 # include "version_variable.h"
   character(len=40)  :: mdl = "MOM_vert_friction" ! This module's name.
   character(len=40)  :: thickness_units
-  real :: Kv_mks ! KVML in MKS
+  real :: Kv_mks ! KVML in MKS [m2 s-1]
 
   if (associated(CS)) then
     call MOM_error(WARNING, "vertvisc_init called with an associated "// &

src/tracer/MOM_CFC_cap.F90

@@ -361,7 +361,7 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C
 
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
-  real :: flux_scale
+  real :: flux_scale ! A dimensional rescaling factor for fluxes [H R-1 Z-1 ~> m3 kg-1 or nondim]
   integer :: i, j, k, is, ie, js, je, nz, m
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke