Corrected the unit documentation for 31 variables

  Corrected the documentation of the units for 31 variables in various modules.
All answers and output are bitwise identical.

src/core/MOM_PressureForce_FV.F90

@@ -153,7 +153,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
                         ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1].
   real :: H_to_RL2_T2   ! A factor to convert from thickness units (H) to pressure
                         ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].
-!  real :: oneatm = 101325.0  ! 1 atm in [Pa] = [kg m-1 s-2]
+!  real :: oneatm       ! 1 standard atmosphere of pressure in [R L2 T-2 ~> Pa]
   real, parameter :: C1_6 = 1.0/6.0
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
@@ -187,6 +187,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
       p(i,j,1) = p_atm(i,j)
     enddo ; enddo
   else
+    ! oneatm = 101325.0 * US%kg_m3_to_R * US%m_s_to_L_T**2 ! 1 atm scaled to [R L2 T-2 ~> Pa]
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       p(i,j,1) = 0.0 ! or oneatm

src/core/MOM_continuity_PPM.F90

@@ -924,8 +924,8 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0,
     FAmt_0, &     ! test velocities [H L ~> m2 or kg m-1].
     uhtot_L, &    ! The summed transport with the westerly (uhtot_L) and
     uhtot_R       ! and easterly (uhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
-  real :: FA_0    ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m].
-  real :: FA_avg  ! The average effective face area [L H ~> m2 or kg m], nominally given by
+  real :: FA_0    ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m-1].
+  real :: FA_avg  ! The average effective face area [L H ~> m2 or kg m-1], nominally given by
                   ! the realized transport divided by the barotropic velocity.
   real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]. This
                        ! limiting is necessary to keep the inverse of visc_rem

src/core/MOM_dynamics_split_RK2.F90

@@ -337,8 +337,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
 
   real, pointer, dimension(:,:,:) :: &
     ! These pointers are used to alter which fields are passed to btstep with various options:
-    u_ptr => NULL(), &   ! A pointer to a zonal velocity [L T-1]
-    v_ptr => NULL(), &   ! A pointer to a meridional velocity [L T-1]
+    u_ptr => NULL(), &   ! A pointer to a zonal velocity [L T-1 ~> m s-1]
+    v_ptr => NULL(), &   ! A pointer to a meridional velocity [L T-1 ~> m s-1]
     uh_ptr => NULL(), &  ! A pointer to a zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
     vh_ptr => NULL(), &  ! A pointer to a meridional volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
     ! These pointers are just used as shorthand for CS%u_av, CS%v_av, and CS%h_av.

src/diagnostics/MOM_wave_structure.F90

@@ -171,20 +171,20 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
   real, dimension(SZK_(GV))   :: dz           !< thicknesses of merged layers (same as Hc I hope) [Z ~> m]
   ! real, dimension(SZK_(GV)+1) :: dWdz_profile !< profile of dW/dz
   real                        :: w2avg   !< average of squared vertical velocity structure funtion [Z ~> m]
-  real                        :: int_dwdz2
-  real                        :: int_w2
-  real                        :: int_N2w2
-  real                        :: KE_term !< terms in vertically averaged energy equation
-  real                        :: PE_term !< terms in vertically averaged energy equation
+  real                        :: int_dwdz2 !< Vertical integral of the square of u_strct [Z ~> m]
+  real                        :: int_w2   !< Vertical integral of the square of w_strct [Z ~> m]
+  real                        :: int_N2w2 !< Vertical integral of N2 [Z T-2 ~> m s-2]
+  real                        :: KE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2]
+  real                        :: PE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2]
   real                        :: W0      !< A vertical velocity magnitude [Z T-1 ~> m s-1]
   real                        :: gp_unscaled !< A version of gprime rescaled to [L T-2 ~> m s-2].
   real, dimension(SZK_(GV)-1) :: lam_z   !< product of eigen value and gprime(k); one value for each
                                          !< interface (excluding surface and bottom)
   real, dimension(SZK_(GV)-1) :: a_diag, b_diag, c_diag
                                          !< diagonals of tridiagonal matrix; one value for each
                                          !< interface (excluding surface and bottom)
-  real, dimension(SZK_(GV)-1) :: e_guess !< guess at eigen vector with unit amplitde (for TDMA)
-  real, dimension(SZK_(GV)-1) :: e_itt   !< improved guess at eigen vector (from TDMA)
+  real, dimension(SZK_(GV)-1) :: e_guess !< guess at eigen vector with unit amplitude (for TDMA) [nondim]
+  real, dimension(SZK_(GV)-1) :: e_itt   !< improved guess at eigen vector (from TDMA) [nondim]
   real    :: Pi
   integer :: kc
   integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop
@@ -523,7 +523,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
 
             ! Back-calculate amplitude from energy equation
             if (present(En) .and. (freq**2*Kmag2 > 0.0)) then
-              ! Units here are [R
+              ! Units here are [R Z ~> kg m-2]
               KE_term = 0.25*GV%Rho0*( ((freq**2 + f2) / (freq**2*Kmag2))*int_dwdz2 + int_w2 )
               PE_term = 0.25*GV%Rho0*( int_N2w2 / freq**2 )
               if (En(i,j) >= 0.0) then

src/initialization/MOM_shared_initialization.F90

@@ -814,7 +814,7 @@ subroutine reset_face_lengths_list(G, param_file, US)
   real, allocatable, dimension(:) :: &
     Dmin_u, Dmax_u, Davg_u   ! Porous barrier monomial fit params [m]
   real, allocatable, dimension(:) :: &
-    Dmin_v, Dmax_v, Davg_v
+    Dmin_v, Dmax_v, Davg_v   ! Porous barrier monomial fit params [m]
   real    :: lat, lon     ! The latitude and longitude of a point.
   real    :: len_lon      ! The periodic range of longitudes, usually 360 degrees.
   real    :: len_lat      ! The range of latitudes, usually 180 degrees.

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -489,7 +489,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
     ! netMassInOut = water [H ~> m or kg m-2] added/removed via surface fluxes
     ! netMassOut   = water [H ~> m or kg m-2] removed via evaporating surface fluxes
     ! net_heat     = heat via surface fluxes [degC H ~> degC m or degC kg m-2]
-    ! net_salt     = salt via surface fluxes [ppt H ~> dppt m or gSalt m-2]
+    ! net_salt     = salt via surface fluxes [ppt H ~> ppt m or gSalt m-2]
     ! Pen_SW_bnd   = components to penetrative shortwave radiation
     call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
                   CS%H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
@@ -1527,7 +1527,7 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, &
   real :: TKE_full_ent  ! The TKE remaining if a layer is fully entrained
                         ! [Z L2 T-2 ~> m3 s-2].
   real :: dRL       ! Work required to mix water from the next layer
-                    ! across the mixed layer [L2 T-2 ~> L2 s-2].
+                    ! across the mixed layer [L2 T-2 ~> m2 s-2].
   real :: Pen_En_Contrib  ! Penetrating SW contributions to the changes in
                           ! TKE, divided by layer thickness in m [L2 T-2 ~> m2 s-2].
   real :: Cpen1     ! A temporary variable [L2 T-2 ~> m2 s-2].

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -1171,7 +1171,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
     !                netMassOut < 0 means mass leaves ocean.
     ! netHeat      = heat via surface fluxes [degC H ~> degC m or degC kg m-2], excluding the part
     !                contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
-    ! netSalt      = surface salt fluxes [ppt H ~> dppt m or gSalt m-2]
+    ! netSalt      = surface salt fluxes [ppt H ~> ppt m or gSalt m-2]
     ! Pen_SW_bnd   = components to penetrative shortwave radiation split according to bands.
     !                This field provides that portion of SW from atmosphere that in fact
     !                enters to the ocean and participates in pentrative SW heating.

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -200,7 +200,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
   real :: hwtot            ! Sum of the thicknesses used to calculate
                            ! the near-bottom velocity magnitude [H ~> m or kg m-2].
   real :: hutot            ! Running sum of thicknesses times the
-                           ! velocity magnitudes [H T T-1 ~> m2 s-1 or kg m-1 s-1].
+                           ! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
   real :: Thtot            ! Running sum of thickness times temperature [degC H ~> degC m or degC kg m-2].
   real :: Shtot            ! Running sum of thickness times salinity [ppt H ~> ppt m or ppt kg m-2].
   real :: hweight          ! The thickness of a layer that is within Hbbl
@@ -597,7 +597,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
     ! When stratification dominates h_N<<h_f, and vice versa.
     do i=is,ie ; if (do_i(i)) then
       ! The 400.0 in this expression is the square of a Ci introduced in KW99, eq. 2.22.
-      ustarsq = Rho0x400_G * ustar(i)**2 ! Note not in units of u*^2 but [H R ~> kg m-2 or kg m-5]
+      ustarsq = Rho0x400_G * ustar(i)**2 ! Note not in units of u*^2 but [H R ~> kg m-2 or kg2 m-5]
       htot = 0.0
 
       ! Calculate the thickness of a stratification limited BBL ignoring rotation:

src/tracer/MOM_lateral_boundary_diffusion.F90

@@ -596,7 +596,7 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_
   real, allocatable :: F_layer_z(:)  !< Diffusive flux at U/V-point in the ztop grid    [H L2 conc ~> m3 conc]
   real              :: h_vel(ke)     !< Thicknesses at u- and v-points in the native grid
                                      !! The harmonic mean is used to avoid zero values      [H ~> m or kg m-2]
-  real    :: khtr_avg                !< Thickness-weighted diffusivity at the velocity-point [L2 T-1 ~> m s-1]
+  real    :: khtr_avg                !< Thickness-weighted diffusivity at the velocity-point [L2 T-1 ~> m2 s-1]
                                      !! This is just to remind developers that khtr_avg should be
                                      !! computed once khtr is 3D.
   real    :: htot                    !< Total column thickness                              [H ~> m or kg m-2]

src/tracer/MOM_tracer_registry.F90

@@ -84,7 +84,8 @@ module MOM_tracer_registry
 !  real, dimension(:,:,:), pointer :: diff_cont_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
 !                                                              !! [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1]
 !  real, dimension(:,:,:), pointer :: diff_conc_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
-!                                                              !! expressed as a change in concentration [conc T-1]
+!                                                              !! expressed as a change in concentration
+!                                                              !! [conc T-1 ~> conc s-1]
   real, dimension(:,:,:), pointer :: t_prev         => NULL() !< tracer concentration array at a previous
                                                               !! timestep used for diagnostics [conc]
   real, dimension(:,:,:), pointer :: Trxh_prev      => NULL() !< layer integrated tracer concentration array

src/user/Idealized_Hurricane.F90

@@ -595,7 +595,7 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C
   U_TS = CS%hurr_translation_spd*0.5*cos(transdir)
   V_TS = CS%hurr_translation_spd*0.5*sin(transdir)
 
-  ! Set the surface wind stresses, in [Pa]. A positive taux
+  ! Set the surface wind stresses, in [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
   ! calculation of ustar - otherwise the lower bound would be Isq.

src/user/MOM_controlled_forcing.F90

@@ -46,7 +46,7 @@ module MOM_controlled_forcing
   real    :: lam_prec       !< A constant of proportionality between SSS anomalies
                             !! (normalised by mean SSS) and precipitation [R Z T-1 ~> kg m-2 s-1]
   real    :: lam_cyc_heat   !< A constant of proportionality between cyclical SST
-                            !! anomalies and corrective heat fluxes [W m-2 degC-1]
+                            !! anomalies and corrective heat fluxes [Q R Z T-1 degC-1 ~> W m-2 degC-1]
   real    :: lam_cyc_prec   !< A constant of proportionality between cyclical SSS
                             !! anomalies (normalised by mean SSS) and corrective
                             !! precipitation [R Z T-1 ~> kg m-2 s-1]
@@ -270,7 +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 s] or [ppt s].
+    ! These loops temporarily change the units of the CS%avg_ variables to [degC T ~> degC s]
+    ! or [ppt 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)

src/user/SCM_CVMix_tests.F90

@@ -36,8 +36,8 @@ module SCM_CVMix_tests
   logical :: UseHeatFlux    !< True to use heat flux
   logical :: UseEvaporation !< True to use evaporation
   logical :: UseDiurnalSW   !< True to use diurnal sw radiation
-  real :: tau_x !< (Constant) Wind stress, X [Pa]
-  real :: tau_y !< (Constant) Wind stress, Y [Pa]
+  real :: tau_x !< (Constant) Wind stress, X [R L Z T-2 ~> Pa]
+  real :: tau_y !< (Constant) Wind stress, Y [R L Z T-2 ~> Pa]
   real :: surf_HF !< (Constant) Heat flux [degC Z T-1 ~> m degC s-1]
   real :: surf_evap !< (Constant) Evaporation rate [Z T-1 ~> m s-1]
   real :: Max_sw !< maximum of diurnal sw radiation [degC Z T-1 ~> degC m s-1]
@@ -56,7 +56,7 @@ subroutine SCM_CVMix_tests_TS_init(T, S, h, G, GV, US, param_file, just_read)
   type(ocean_grid_type),                     intent(in)  :: G  !< Grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< Vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T  !< Potential temperature [degC]
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [psu]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [ppt]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h  !< Layer thickness [H ~> m or kg m-2]
   type(unit_scale_type),                     intent(in)  :: US !< A dimensional unit scaling type
   type(param_file_type),                     intent(in)  :: param_file !< Input parameter structure

src/user/baroclinic_zone_initialization.F90

@@ -36,7 +36,7 @@ subroutine bcz_params(G, GV, US, param_file, S_ref, dSdz, delta_S, dSdx, T_ref,
   real,                    intent(out) :: S_ref      !< Reference salinity [ppt]
   real,                    intent(out) :: dSdz       !< Salinity stratification [ppt Z-1 ~> ppt m-1]
   real,                    intent(out) :: delta_S    !< Salinity difference across baroclinic zone [ppt]
-  real,                    intent(out) :: dSdx       !< Linear salinity gradient [ppt m-1]
+  real,                    intent(out) :: dSdx       !< Linear salinity gradient [ppt G%xaxis_units-1]
   real,                    intent(out) :: T_ref      !< Reference temperature [degC]
   real,                    intent(out) :: dTdz       !< Temperature stratification [degC Z-1 ~> degC m-1]
   real,                    intent(out) :: delta_T    !< Temperature difference across baroclinic zone [degC]