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 f8c46c8926..8e04c5f116 100644
--- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
+++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
@@ -85,14 +85,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-2 ~> Pa]
+  real :: gust_const            !< Constant unresolved background gustiness for ustar [R Z2 T-2 ~> Pa]
   logical :: read_gust_2d       !< If true, use a 2-dimensional gustiness supplied from an input file.
   real, pointer, dimension(:,:) :: &
     BBL_tidal_dis => NULL()     !< Tidal energy dissipation in the bottom boundary layer that can act as a
                                 !! source of energy for bottom boundary layer mixing [R Z L2 T-3 ~> W m-2]
   real, pointer, dimension(:,:) :: &
     gust => NULL()              !< A spatially varying unresolved background gustiness that
-                                !! contributes to ustar [R L Z T-2 ~> Pa].  gust is used when read_gust_2d is true.
+                                !! contributes to ustar [R Z2 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]
@@ -698,7 +698,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
     rigidity_at_h, &  ! Ice rigidity at tracer points [L4 Z-1 T-1 ~> m3 s-1]
     net_mass_src, &   ! A temporary of net mass sources [R Z T-1 ~> kg m-2 s-1].
     ustar_tmp, &      ! A temporary array of ustar values [Z T-1 ~> m s-1].
-    tau_mag_tmp       ! A temporary array of surface stress magnitudes [R Z L T-2 ~> Pa]
+    tau_mag_tmp       ! A temporary array of surface stress magnitudes [R Z2 T-2 ~> Pa]
 
   real :: I_GEarth      ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-1]
   real :: Kv_rho_ice    ! (CS%Kv_sea_ice / CS%density_sea_ice) [L4 Z-2 T-1 R-1 ~> m5 s-1 kg-1]
@@ -939,10 +939,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
                                                  !! any contributions from gustiness [Z T-1 ~> m s-1].
   real, dimension(SZI_(G),SZJ_(G)), &
                  optional, intent(inout) :: mag_tau !< The magintude of the wind stress at tracer points
-                                                 !! including subgridscale variability and gustiness [R Z L T-2 ~> Pa]
+                                                 !! including subgridscale variability and gustiness [R Z2 T-2 ~> Pa]
   real, dimension(SZI_(G),SZJ_(G)), &
                  optional, intent(out) :: gustless_mag_tau !< The magintude of the wind stress at tracer points
-                                                 !! without any contributions from gustiness [R Z L T-2 ~> Pa]
+                                                 !! without any contributions from gustiness [R Z2 T-2 ~> Pa]
   integer,       optional, intent(in)    :: tau_halo !< The halo size of wind stresses to set, 0 by default.
 
   ! Local variables
@@ -953,11 +953,12 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   real, dimension(SZIB_(G),SZJB_(G)) :: taux_in_B ! Zonal wind stresses [R Z L T-2 ~> Pa] at q points
   real, dimension(SZIB_(G),SZJB_(G)) :: tauy_in_B ! Meridional wind stresses [R Z L T-2 ~> Pa] at q points
 
-  real :: gustiness     ! unresolved gustiness that contributes to ustar [R Z L T-2 ~> Pa]
-  real :: Irho0         ! Inverse of the mean density rescaled to [Z L-1 R-1 ~> m3 kg-1]
+  real :: gustiness     ! unresolved gustiness that contributes to ustar [R Z2 T-2 ~> Pa]
+  real :: Irho0         ! Inverse of the Boussinesq mean density [R-1 ~> m3 kg-1]
   real :: taux2, tauy2  ! squared wind stresses [R2 Z2 L2 T-4 ~> Pa2]
-  real :: tau_mag       ! magnitude of the wind stress [R Z L T-2 ~> Pa]
+  real :: tau_mag       ! magnitude of the wind stress [R Z2 T-2 ~> Pa]
   real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1]
+  real :: Pa_to_RZ2_T2  ! The combination of unit conversion factors used for mag_tau [R Z2 T-2 Pa-1 ~> 1]
 
   logical :: do_ustar, do_gustless, do_tau_mag, do_gustless_tau_mag
   integer :: wind_stagger  ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains)
@@ -969,7 +970,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   Isqh = G%IscB-halo ; Ieqh  = G%IecB+halo ; Jsqh = G%JscB-halo ; Jeqh = G%JecB+halo
   i0 = is - index_bounds(1) ; j0 = js - index_bounds(3)
 
-  IRho0 = US%L_to_Z / CS%Rho0
+  IRho0 = 1.0 / CS%Rho0
   stress_conversion = US%Pa_to_RLZ_T2 * CS%wind_stress_multiplier
 
   do_ustar = present(ustar) ; do_gustless = present(gustless_ustar)
@@ -1073,6 +1074,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
     ! parametizations. The background gustiness (for example with a relatively small value
     ! of 0.02 Pa) is intended to give reasonable behavior in regions of very weak winds.
     if (associated(IOB%stress_mag)) then
+      Pa_to_RZ2_T2 = US%Pa_to_RLZ_T2 * US%L_to_Z
+
       if (do_ustar .or. do_tau_mag) then ; do j=js,je ; do i=is,ie
         gustiness = CS%gust_const
         if (CS%read_gust_2d) then
@@ -1084,19 +1087,19 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
             gustiness = CS%gust(i,j)
         endif
         if (do_tau_mag) &
-          mag_tau(i,j) = gustiness + US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0)
+          mag_tau(i,j) = gustiness + Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0)
         if (do_gustless_tau_mag) &
-          gustless_mag_tau(i,j) = US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0)
+          gustless_mag_tau(i,j) = Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0)
         if (do_ustar) &
-          ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0))
+          ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0))
       enddo ; enddo ; endif
       if (CS%answer_date < 20190101) then
         if (do_gustless) then ; do j=js,je ; do i=is,ie
-          gustless_ustar(i,j) = sqrt(US%Pa_to_RLZ_T2*US%L_to_Z*IOB%stress_mag(i-i0,j-j0) / CS%Rho0)
+          gustless_ustar(i,j) = sqrt(Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0) / CS%Rho0)
         enddo ; enddo ; endif
       else
         if (do_gustless) then ; do j=js,je ; do i=is,ie
-          gustless_ustar(i,j) = sqrt(IRho0 * US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0))
+          gustless_ustar(i,j) = sqrt(IRho0 * Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0))
         enddo ; enddo ; endif
       endif
     elseif (wind_stagger == BGRID_NE) then
@@ -1104,7 +1107,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         tau_mag = 0.0 ; gustiness = CS%gust_const
         if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + &
              (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0) then
-          tau_mag = sqrt(((G%mask2dBu(I,J)*((taux_in_B(I,J)**2) + (tauy_in_B(I,J)**2)) + &
+          tau_mag = US%L_to_Z * sqrt(((G%mask2dBu(I,J)*((taux_in_B(I,J)**2) + (tauy_in_B(I,J)**2)) + &
               G%mask2dBu(I-1,J-1)*((taux_in_B(I-1,J-1)**2) + (tauy_in_B(I-1,J-1)**2))) + &
              (G%mask2dBu(I,J-1)*((taux_in_B(I,J-1)**2) + (tauy_in_B(I,J-1)**2)) + &
               G%mask2dBu(I-1,J)*((taux_in_B(I-1,J)**2) + (tauy_in_B(I-1,J)**2))) ) / &
@@ -1115,21 +1118,21 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
         if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
         if (CS%answer_date < 20190101) then
-          if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
+          if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0)
         else
           if (do_gustless) gustless_ustar(i,j) = sqrt(IRho0 * tau_mag)
         endif
       enddo ; enddo
     elseif (wind_stagger == AGRID) then
       do j=js,je ; do i=is,ie
-        tau_mag = G%mask2dT(i,j) * sqrt((taux_in_A(i,j)**2) + (tauy_in_A(i,j)**2))
+        tau_mag = G%mask2dT(i,j) * US%L_to_Z * sqrt((taux_in_A(i,j)**2) + (tauy_in_A(i,j)**2))
         gustiness = CS%gust_const
         if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j)
         if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag)
         if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
         if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
         if (CS%answer_date < 20190101) then
-          if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
+          if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0)
         else
           if (do_gustless) gustless_ustar(i,j) = sqrt(IRho0 * tau_mag)
         endif
@@ -1143,7 +1146,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0.0) &
           tauy2 = (G%mask2dCv(i,J-1)*(tauy_in_C(i,J-1)**2) + G%mask2dCv(i,J)*(tauy_in_C(i,J)**2)) / &
                   (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
-        tau_mag = sqrt(taux2 + tauy2)
+        tau_mag = US%L_to_Z * sqrt(taux2 + tauy2)
 
         gustiness = CS%gust_const
         if (CS%read_gust_2d) gustiness = CS%gust(i,j)
@@ -1152,7 +1155,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
         if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
         if (CS%answer_date < 20190101) then
-          if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
+          if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0)
         else
           if (do_gustless) gustless_ustar(i,j) = sqrt(IRho0 * tau_mag)
         endif
@@ -1616,7 +1619,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
                  "an input file", default=.false.)
   call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
                  "The background gustiness in the winds.", &
-                 units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2)
+                 units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z)
   if (CS%read_gust_2d) then
     call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
                  "The file in which the wind gustiness is found in "//&
@@ -1627,7 +1630,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
     ! NOTE: There are certain cases where FMS is unable to read this file, so
     ! we use read_netCDF_data in place of MOM_read_data.
     call read_netCDF_data(gust_file, 'gustiness', CS%gust, G%Domain, &
-                          rescale=US%Pa_to_RLZ_T2) ! units in file should be [Pa]
+                          rescale=US%Pa_to_RLZ_T2*US%L_to_Z) ! units in file should be [Pa]
   endif
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90
index 9c4359bf60..e3b7b0cec7 100644
--- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90
+++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90
@@ -32,7 +32,6 @@ module ocean_model_mod
 use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type
 use MOM_forcing_type, only : forcing, mech_forcing, allocate_forcing_type
 use MOM_forcing_type, only : fluxes_accumulate, get_net_mass_forcing
-use MOM_forcing_type, only : copy_back_forcing_fields
 use MOM_forcing_type, only : forcing_diagnostics, mech_forcing_diags
 use MOM_get_input, only : Get_MOM_Input, directories
 use MOM_grid, only : ocean_grid_type
diff --git a/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90
index b0562851ea..e5c5943d4f 100644
--- a/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90
+++ b/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90
@@ -774,7 +774,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
              ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) )
         if (CS%read_gust_2d) gustiness = CS%gust(i,j)
       endif
-      forces%tau_mag(i,j) = gustiness + tau_mag
+      forces%tau_mag(i,j) = US%L_to_Z*(gustiness + tau_mag)
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag)
     enddo ; enddo
 
@@ -800,7 +800,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
     do j=js,je ; do i=is,ie
       gustiness = CS%gust_const
       if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j)
-      forces%tau_mag(i,j) = gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2))
+      forces%tau_mag(i,j) = US%L_to_Z*(gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2)))
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * &
                                sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2)))
     enddo ; enddo
@@ -822,10 +822,10 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
                  G%mask2dCv(i,J)*(forces%tauy(i,J)**2)) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
 
       if (CS%read_gust_2d) then
-        forces%tau_mag(i,j) = CS%gust(i,j) + sqrt(taux2 + tauy2)
+        forces%tau_mag(i,j) = US%L_to_Z*(CS%gust(i,j) + sqrt(taux2 + tauy2))
         forces%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2))
       else
-        forces%tau_mag(i,j) = CS%gust_const + sqrt(taux2 + tauy2)
+        forces%tau_mag(i,j) = US%L_to_Z*(CS%gust_const + sqrt(taux2 + tauy2))
         forces%ustar(i,j) = sqrt(CS%gust_const*Irho0 + Irho0*sqrt(taux2 + tauy2))
       endif
     enddo ; enddo
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 960c15f39b..4287dfeb43 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
@@ -835,7 +835,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
           ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) )
         if (CS%read_gust_2d) gustiness = CS%gust(i,j)
       endif
-      forces%tau_mag(i,j) = gustiness + tau_mag
+      forces%tau_mag(i,j) = US%L_to_Z*(gustiness + tau_mag)
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag)
     enddo ; enddo
     call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1)
@@ -861,7 +861,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
     do j=js,je ; do i=is,ie
       gustiness = CS%gust_const
       if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j)
-      forces%tau_mag(i,j) = gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2))
+      forces%tau_mag(i,j) = US%L_to_Z*(gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2)))
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * &
                                sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2)))
       !forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j))
@@ -884,10 +884,10 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
                  G%mask2dCv(i,J)*(forces%tauy(i,J)**2)) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
 
       if (CS%read_gust_2d) then
-        forces%tau_mag(i,j) = CS%gust(i,j) + sqrt(taux2 + tauy2)
+        forces%tau_mag(i,j) = US%L_to_Z*(CS%gust(i,j) + sqrt(taux2 + tauy2))
         forces%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2))
       else
-        forces%tau_mag(i,j) = CS%gust_const + sqrt(taux2 + tauy2)
+        forces%tau_mag(i,j) = US%L_to_Z*(CS%gust_const + sqrt(taux2 + tauy2))
         forces%ustar(i,j) = sqrt(CS%gust_const*Irho0 + Irho0*sqrt(taux2 + tauy2))
       endif
     enddo ; enddo
diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90
index 448b16b5e4..81edcdb7c3 100644
--- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90
@@ -92,7 +92,7 @@ module MOM_surface_forcing
   real :: taux_mag              !< Peak magnitude of the zonal wind stress for several analytic
                                 !! profiles [R L Z T-2 ~> Pa]
 
-  real    :: gust_const                 !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa]
+  real    :: gust_const                 !< constant unresolved background gustiness for ustar [R Z2 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-2 ~> Pa]
                                         !! gust is used when read_gust_2d is true.
@@ -419,14 +419,14 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS)
   type(surface_forcing_CS), pointer       :: CS   !< pointer to control structure returned by
                                                   !! a previous surface_forcing_init call
   ! Local variables
-  real :: mag_tau  ! Magnitude of the wind stress [R Z L T-2 ~> Pa]
+  real :: mag_tau  ! Magnitude of the wind stress [R Z2 T-2 ~> Pa]
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
 
   call callTree_enter("wind_forcing_const, MOM_surface_forcing.F90")
   is   = G%isc  ; ie   = G%iec  ; js   = G%jsc  ; je   = G%jec
   Isq  = G%IscB ; Ieq  = G%IecB ; Jsq  = G%JscB ; Jeq  = G%JecB
 
-  mag_tau = sqrt( tau_x0**2 + tau_y0**2)
+  mag_tau = US%L_to_Z * sqrt( tau_x0**2 + tau_y0**2)
 
   ! Set the steady surface wind stresses, in units of [R L Z T-2 ~> Pa].
   do j=js,je ; do I=is-1,Ieq
@@ -439,14 +439,14 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS)
 
   if (CS%read_gust_2d) then
     if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = sqrt( US%L_to_Z * ( mag_tau + CS%gust(i,j) ) / CS%Rho0 )
+      forces%ustar(i,j) = sqrt( ( mag_tau + CS%gust(i,j) ) / CS%Rho0 )
     enddo ; enddo ; endif
     if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
       forces%tau_mag(i,j) = mag_tau + CS%gust(i,j)
     enddo ; enddo ; endif
   else
     if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = sqrt( US%L_to_Z * ( mag_tau + CS%gust_const ) / CS%Rho0 )
+      forces%ustar(i,j) = sqrt( ( mag_tau + CS%gust_const ) / CS%Rho0 )
     enddo ; enddo ; endif
     if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
       forces%tau_mag(i,j) = mag_tau + CS%gust_const
@@ -562,13 +562,13 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)
   ! set the friction velocity
   if (CS%answer_date < 20190101) then
     if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
-      forces%tau_mag(i,j) = CS%gust_const + sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                                                      ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+      forces%tau_mag(i,j) = CS%gust_const + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                                                  ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
     enddo ; enddo ; endif
     if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = sqrt(US%L_to_Z * ((CS%gust_const/CS%Rho0) + &
-              sqrt(0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2) + &
-                        (forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))/CS%Rho0) )
+      forces%ustar(i,j) = sqrt( (CS%gust_const/CS%Rho0) + &
+              US%L_to_Z * sqrt(0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2) + &
+                                    (forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))/CS%Rho0 )
     enddo ; enddo ; endif
   else
     call stresses_to_ustar(forces, G, US, CS)
@@ -710,7 +710,7 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
   real    :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R L Z T-2 ~> Pa]
   real    :: ustar_loc(SZI_(G),SZJ_(G)) ! The local value of ustar [Z T-1 ~> m s-1]
   real    :: tau_mag    ! The magnitude of the wind stress including any contributions from
-                        ! sub-gridscale variability or gustiness [R L Z T-2 ~> Pa]
+                        ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa]
   integer :: time_lev                ! The time level that is used for a field.
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
   logical :: read_Ustar
@@ -745,19 +745,19 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
       if (.not.read_Ustar) then
         if (CS%read_gust_2d) then
           if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
-            forces%tau_mag(i,j) = CS%gust(i,j) + sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
+            forces%tau_mag(i,j) = CS%gust(i,j) + US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
           enddo ; enddo ; endif
           if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
-            tau_mag = CS%gust(i,j) + sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
-            forces%ustar(i,j) = sqrt(tau_mag * US%L_to_Z / CS%Rho0)
+            tau_mag = CS%gust(i,j) + US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
+            forces%ustar(i,j) = sqrt(tau_mag / CS%Rho0)
           enddo ; enddo ; endif
         else
           if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
-            forces%tau_mag(i,j) = CS%gust_const + sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
+            forces%tau_mag(i,j) = CS%gust_const + US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
           enddo ; enddo ; endif
           if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
-            forces%ustar(i,j) = sqrt(US%L_to_Z * (CS%gust_const/CS%Rho0 + &
-                    sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) / CS%Rho0) )
+            forces%ustar(i,j) = sqrt( CS%gust_const/CS%Rho0 + &
+                    US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) / CS%Rho0 )
           enddo ; enddo ; endif
         endif
       endif
@@ -799,25 +799,25 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
         if (CS%read_gust_2d) then
           if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
             forces%tau_mag(i,j) = CS%gust(i,j) + &
-                    sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                              ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+                    US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                          ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
           enddo ; enddo ; endif
           if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
             tau_mag = CS%gust(i,j) + &
-                    sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                              ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
-            forces%ustar(i,j) = sqrt( tau_mag * US%L_to_Z / CS%Rho0 )
+                   US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                         ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+            forces%ustar(i,j) = sqrt( tau_mag / CS%Rho0 )
           enddo ; enddo ; endif
         else
           if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
             forces%tau_mag(i,j) = CS%gust_const + &
-                  sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                            ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+                  US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                        ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
           enddo ; enddo ; endif
           if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
-             forces%ustar(i,j) = sqrt(US%L_to_Z * ( (CS%gust_const/CS%Rho0) + &
-                  sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                            ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))/CS%Rho0))
+             forces%ustar(i,j) = sqrt( CS%gust_const/CS%Rho0 + &
+                  US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                        ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))/CS%Rho0 )
           enddo ; enddo ; endif
         endif
       endif
@@ -830,7 +830,7 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
       call MOM_read_data(filename, CS%Ustar_var, ustar_loc(:,:), &
                          G%Domain, timelevel=time_lev, scale=US%m_to_Z*US%T_to_s)
       if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
-        forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * ustar_loc(i,j)**2
+        forces%tau_mag(i,j) = CS%Rho0 * ustar_loc(i,j)**2
       enddo ; enddo ; endif
       if (associated(forces%ustar)) then ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
         forces%ustar(i,j) = ustar_loc(i,j)
@@ -861,7 +861,7 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS)
   real :: ustar_prev(SZI_(G),SZJ_(G)) ! The pre-override value of ustar [Z T-1 ~> m s-1]
   real :: ustar_loc(SZI_(G),SZJ_(G)) ! The value of ustar, perhaps altered by data override [Z T-1 ~> m s-1]
   real :: tau_mag       ! The magnitude of the wind stress including any contributions from
-                        ! sub-gridscale variability or gustiness [R L Z T-2 ~> Pa]
+                        ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa]
   integer :: i, j
 
   call callTree_enter("wind_forcing_by_data_override, MOM_surface_forcing.F90")
@@ -885,24 +885,24 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS)
   enddo ; enddo
 
   if (CS%read_gust_2d) then
-    call data_override(G%Domain, 'gust', CS%gust, day, scale=US%Pa_to_RLZ_T2)
+    call data_override(G%Domain, 'gust', CS%gust, day, scale=US%Pa_to_RLZ_T2*US%L_to_Z)
     if (associated(forces%tau_mag)) then ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      forces%tau_mag(i,j) = sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j)
+      forces%tau_mag(i,j) = US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j)
     enddo ; enddo ; endif
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      tau_mag = sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j)
-      ustar_loc(i,j) = sqrt( tau_mag * US%L_to_Z / CS%Rho0 )
+      tau_mag = US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j)
+      ustar_loc(i,j) = sqrt( tau_mag / CS%Rho0 )
     enddo ; enddo
   else
     if (associated(forces%tau_mag)) then
       do j=G%jsc,G%jec ; do i=G%isc,G%iec
-        forces%tau_mag(i,j) = sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust_const
-      ! ustar_loc(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 )
+        forces%tau_mag(i,j) = US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust_const
+      ! ustar_loc(i,j) = sqrt( forces%tau_mag(i,j) / CS%Rho0 )
       enddo ; enddo
     endif
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      ustar_loc(i,j) = sqrt(US%L_to_Z * (sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))/CS%Rho0 + &
-                            CS%gust_const/CS%Rho0))
+      ustar_loc(i,j) = sqrt(US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))/CS%Rho0 + &
+                            CS%gust_const/CS%Rho0)
     enddo ; enddo
   endif
 
@@ -913,7 +913,7 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS)
   ! Only reset values where data override of ustar has occurred
   if (associated(forces%tau_mag)) then
     do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (ustar_prev(i,j) /= ustar_loc(i,j)) then
-      forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * ustar_loc(i,j)**2
+      forces%tau_mag(i,j) = CS%Rho0 * ustar_loc(i,j)**2
     endif ; enddo ; enddo
   endif
 
@@ -934,38 +934,37 @@ subroutine stresses_to_ustar(forces, G, US, CS)
   type(surface_forcing_CS), pointer       :: CS     !< pointer to control structure returned by
                                                     !! a previous surface_forcing_init call
   ! Local variables
-  real :: I_rho         ! The inverse of the reference density times a ratio of scaling
-                        ! factors [Z L-1 R-1 ~> m3 kg-1]
+  real :: I_rho         ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]
   real :: tau_mag       ! The magnitude of the wind stress including any contributions from
-                        ! sub-gridscale variability or gustiness [R L Z T-2 ~> Pa]
+                        ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa]
   integer :: i, j, is, ie, js, je
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
-  I_rho = US%L_to_Z / CS%Rho0
+  I_rho = 1.0 / CS%Rho0
 
   if (CS%read_gust_2d) then
     if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
       forces%tau_mag(i,j) = CS%gust(i,j) + &
-              sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                        ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+              US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                    ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
     enddo ; enddo ; endif
     if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
       tau_mag = CS%gust(i,j) + &
-              sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                        ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+              US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                    ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
       forces%ustar(i,j) = sqrt( tau_mag * I_rho )
     enddo ; enddo ; endif
   else
     if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
       forces%tau_mag(i,j) = CS%gust_const + &
-              sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                        ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+              US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                    ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
     enddo ; enddo ; endif
     if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
       tau_mag = CS%gust_const + &
-              sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
-                        ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
+              US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + &
+                                    ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))
       forces%ustar(i,j) = sqrt( tau_mag * I_rho )
     enddo ; enddo ; endif
   endif
@@ -2007,7 +2006,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
 
   call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
                  "The background gustiness in the winds.", &
-                 units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2)
+                 units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z)
 
   call get_param(param_file, mdl, "USTAR_GUSTLESS_BUG", CS%ustar_gustless_bug, &
                  "If true include a bug in the time-averaging of the gustless wind friction velocity", &
@@ -2047,7 +2046,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
     ! NOTE: There are certain cases where FMS is unable to read this file, so
     ! we use read_netCDF_data in place of MOM_read_data.
     call read_netCDF_data(filename, 'gustiness', CS%gust, G%Domain, &
-                          rescale=US%Pa_to_RLZ_T2) ! units in file should be [Pa]
+                          rescale=US%Pa_to_RLZ_T2*US%L_to_Z) ! units in file should be [Pa]
   endif
 
 !  All parameter settings are now known.
diff --git a/config_src/drivers/solo_driver/user_surface_forcing.F90 b/config_src/drivers/solo_driver/user_surface_forcing.F90
index 559291b225..55b1be1172 100644
--- a/config_src/drivers/solo_driver/user_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/user_surface_forcing.F90
@@ -38,7 +38,7 @@ module user_surface_forcing
   real :: rho_restore        !< The density that is used to convert piston velocities into salt
                              !! or heat fluxes with salinity or temperature restoring [R ~> kg m-3]
   real :: gust_const         !< A constant unresolved background gustiness
-                             !! that contributes to ustar [R L Z T-2 ~> Pa].
+                             !! that contributes to ustar [R Z2 T-2 ~> Pa].
 
   type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                              !! timing of diagnostic output.
@@ -91,10 +91,10 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS)
   if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
     !  This expression can be changed if desired, but need not be.
     forces%tau_mag(i,j) = G%mask2dT(i,j) * (CS%gust_const + &
-            sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + &
-                 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2))))
+            US%L_to_Z*sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + &
+                           0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2))))
     if (associated(forces%ustar)) &
-      forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(forces%tau_mag(i,j) * (US%L_to_Z/CS%Rho0))
+      forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(forces%tau_mag(i,j) * (1.0/CS%Rho0))
   enddo ; enddo ; endif
 
 end subroutine USER_wind_forcing
@@ -275,7 +275,7 @@ subroutine USER_surface_forcing_init(Time, G, US, param_file, diag, CS)
                  units="kg m-3", default=1035.0, scale=US%kg_m3_to_R)
   call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
                  "The background gustiness in the winds.", &
-                 units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2)
+                 units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z)
 
   call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, &
                  "If true, the buoyancy fluxes drive the model back "//&
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index abdf7828ee..f37e54f621 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -79,12 +79,14 @@ module MOM_forcing_type
     ustar         => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1].
     tau_mag       => NULL(), & !< Magnitude of the wind stress averaged over tracer cells,
                                !! including any contributions from sub-gridscale variability
-                               !! or gustiness [R L Z T-2 ~> Pa]
+                               !! or gustiness, rescaled to units that are more convenient for
+                               !! calculating turbulent fluxes and friction velocities [R Z2 T-2 ~> Pa]
     ustar_gustless => NULL(), & !< surface friction velocity scale without any
                                !! any augmentation for gustiness [Z T-1 ~> m s-1].
     tau_mag_gustless => NULL() !< Magnitude of the wind stress averaged over tracer cells,
                                !! without any augmentation for sub-gridscale variability
-                               !! or gustiness [R L Z T-2 ~> Pa]
+                               !! or gustiness, rescaled to units that are more convenient for
+                               !! calculating turbulent fluxes and friction velocities [R Z2 T-2 ~> Pa]
 
   ! surface buoyancy force, used when temperature is not a state variable
   real, pointer, dimension(:,:) :: &
@@ -1132,9 +1134,9 @@ subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units)
                                                !! of [H T-1 ~> m s-1 or kg m-2 s-1]
 
   ! Local variables
-  real :: I_rho        ! The inverse of the reference density times a ratio of scaling
-                       ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases
-                       ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3]
+  real :: I_rho        ! The inverse of the reference density [R-1 ~> m3 kg-1]
+                       ! or in some semi-Boussinesq cases the reference
+                       ! density [H2 R-1 ~> m3 kg-1 or kg m-3]
   logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is
                        ! returned in [H T-1 ~> m s-1 or kg m-2 s-1]
   integer :: i, j, k, is, ie, js, je, hs
@@ -1164,16 +1166,16 @@ subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units)
         "find_ustar_fluxes called in non-Boussinesq mode with insufficient valid values of SpV_avg.")
     if (Z_T_units) then
       do j=js,je ; do i=is,ie
-        U_star(i,j) = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1))
+        U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1))
       enddo ; enddo
     else
       do j=js,je ; do i=is,ie
-        U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1))
+        U_star(i,j) = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1))
       enddo ; enddo
     endif
   else
-    I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H
-    if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0
+    I_rho = GV%Z_to_H * GV%RZ_to_H
+    if (Z_T_units) I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0
     do j=js,je ; do i=is,ie
       U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * I_rho)
     enddo ; enddo
@@ -1198,9 +1200,8 @@ subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_unit
                                                !! of [H T-1 ~> m s-1 or kg m-2 s-1]
 
   ! Local variables
-  real :: I_rho        ! The inverse of the reference density times a ratio of scaling
-                       ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases
-                       ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3]
+  real :: I_rho        ! The inverse of the reference density [R-1 ~> m3 kg-1] or in some semi-Boussinesq cases
+                       ! the rescaled reference density [H2 R-1 ~> m3 kg-1 or kg m-3]
   logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is
                        ! returned in [H T-1 ~> m s-1 or kg m-2 s-1]
   integer :: i, j, k, is, ie, js, je, hs
@@ -1230,16 +1231,16 @@ subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_unit
         "find_ustar_mech called in non-Boussinesq mode with insufficient valid values of SpV_avg.")
     if (Z_T_units) then
       do j=js,je ; do i=is,ie
-        U_star(i,j) = sqrt(US%L_to_Z*forces%tau_mag(i,j) * tv%SpV_avg(i,j,1))
+        U_star(i,j) = sqrt(forces%tau_mag(i,j) * tv%SpV_avg(i,j,1))
       enddo ; enddo
     else
       do j=js,je ; do i=is,ie
-        U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*forces%tau_mag(i,j) / tv%SpV_avg(i,j,1))
+        U_star(i,j) = GV%RZ_to_H * sqrt(forces%tau_mag(i,j) / tv%SpV_avg(i,j,1))
       enddo ; enddo
     endif
   else
-    I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H
-    if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0
+    I_rho = GV%Z_to_H * GV%RZ_to_H
+    if (Z_T_units) I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0
     do j=js,je ; do i=is,ie
       U_star(i,j) = sqrt(forces%tau_mag(i,j) * I_rho)
     enddo ; enddo
@@ -1266,7 +1267,7 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
   if (associated(fluxes%ustar)) &
     call hchksum(fluxes%ustar, mesg//" fluxes%ustar", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T)
   if (associated(fluxes%tau_mag)) &
-    call hchksum(fluxes%tau_mag, mesg//" fluxes%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa)
+    call hchksum(fluxes%tau_mag, mesg//" fluxes%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa*US%Z_to_L)
   if (associated(fluxes%buoy)) &
     call hchksum(fluxes%buoy, mesg//" fluxes%buoy ", G%HI, haloshift=hshift, unscale=US%L_to_m**2*US%s_to_T**3)
   if (associated(fluxes%sw)) &
@@ -1373,7 +1374,7 @@ subroutine MOM_mech_forcing_chksum(mesg, forces, G, US, haloshift)
   if (associated(forces%ustar)) &
     call hchksum(forces%ustar, mesg//" forces%ustar", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T)
   if (associated(forces%tau_mag)) &
-    call hchksum(forces%tau_mag, mesg//" forces%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa)
+    call hchksum(forces%tau_mag, mesg//" forces%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa*US%Z_to_L)
   if (associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) &
     call uvchksum(mesg//" forces%rigidity_ice_[uv]", forces%rigidity_ice_u, &
         forces%rigidity_ice_v, G%HI, haloshift=hshift, symmetric=.true., &
@@ -1503,7 +1504,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
 
   handles%id_tau_mag = register_diag_field('ocean_model', 'tau_mag', diag%axesT1, Time, &
         'Average magnitude of the wind stress including contributions from gustiness', &
-        'Pa', conversion=US%RLZ_T2_to_Pa)
+        'Pa', conversion=US%RLZ_T2_to_Pa*US%Z_to_L)
 
   handles%id_ustar = register_diag_field('ocean_model', 'ustar', diag%axesT1, Time, &
       'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', &
@@ -2460,7 +2461,7 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, US, Rho0)
         endif
       endif
       if (associated(fluxes%tau_mag_gustless)) then
-        fluxes%tau_mag_gustless(i,j) = sqrt(taux2 + tauy2)
+        fluxes%tau_mag_gustless(i,j) = US%L_to_Z*sqrt(taux2 + tauy2)
       endif
     enddo ; enddo
   endif
@@ -3830,14 +3831,14 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar)
                                                  !! or updated from mean tau.
 
   real :: tx_mean, ty_mean ! Mean wind stresses [R L Z T-2 ~> Pa]
-  real :: tau_mag      ! The magnitude of the wind stresses [R L Z T-2 ~> Pa]
-  real :: Irho0        ! Inverse of the mean density rescaled to [Z L-1 R-1 ~> m3 kg-1]
+  real :: tau_mag      ! The magnitude of the wind stresses [R Z2 T-2 ~> Pa]
+  real :: Irho0        ! Inverse of the mean density [R-1 ~> m3 kg-1]
   logical :: do_stress, do_ustar, do_taumag, do_shelf, do_press, do_iceberg, tau2ustar
   integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB
 
-  Irho0 = US%L_to_Z / Rho0
+  Irho0 = 1.0 / Rho0
 
   tau2ustar = .false.
   if (present(UpdateUstar)) tau2ustar = UpdateUstar
@@ -3855,7 +3856,7 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar)
       if (G%mask2dCv(i,J) > 0.0) forces%tauy(i,J) = ty_mean
     enddo ; enddo
     if (tau2ustar) then
-      tau_mag = sqrt((tx_mean**2) + (ty_mean**2))
+      tau_mag = US%L_to_Z*sqrt((tx_mean**2) + (ty_mean**2))
       if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         forces%tau_mag(i,j) = tau_mag
       endif ; enddo ; enddo ; endif
@@ -3866,13 +3867,13 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar)
       if (associated(forces%ustar)) &
         call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T)
       if (associated(forces%tau_mag)) &
-        call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa)
+        call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L)
     endif
   else
     if (associated(forces%ustar)) &
       call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T)
     if (associated(forces%tau_mag)) &
-      call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa)
+      call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L)
   endif
 
   if (do_shelf) then
@@ -3915,9 +3916,9 @@ subroutine homogenize_forcing(fluxes, G, GV, US)
     call homogenize_field_t(fluxes%ustar_gustless, G, tmp_scale=US%Z_to_m*US%s_to_T)
 
   if (associated(fluxes%tau_mag)) &
-    call homogenize_field_t(fluxes%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa)
+    call homogenize_field_t(fluxes%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L)
   if (associated(fluxes%tau_mag_gustless)) &
-    call homogenize_field_t(fluxes%tau_mag_gustless, G, tmp_scale=US%RLZ_T2_to_Pa)
+    call homogenize_field_t(fluxes%tau_mag_gustless, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L)
 
   if (do_water) then
     call homogenize_field_t(fluxes%evap, G, tmp_scale=US%RZ_T_to_kg_m2s)
diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
index d9d9cf6b14..d1fcd43b3c 100644
--- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
+++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
@@ -1609,8 +1609,8 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_F
       TKE(i) = (dt*CS%mstar)*((GV%Z_to_H*(U_star*U_Star*U_Star))*exp_kh) + &
                (exp_kh * dKE_conv + nstar_FC*Conv_En(i) + nstar_CA * TKE_CA)
     else
-      ! Note that GV%Z_to_H*U_star**3 = GV%RZ_to_H * US%L_to_Z*fluxes%tau_mag(i,j) * U_star
-      TKE(i) = (dt*CS%mstar) * ((GV%RZ_to_H*US%L_to_Z * fluxes%tau_mag(i,j) * U_star)*exp_kh) + &
+      ! Note that GV%Z_to_H*U_star**3 = GV%RZ_to_H * fluxes%tau_mag(i,j) * U_star
+      TKE(i) = (dt*CS%mstar) * ((GV%RZ_to_H * fluxes%tau_mag(i,j) * U_star)*exp_kh) + &
                (exp_kh * dKE_conv + nstar_FC*Conv_En(i) + nstar_CA * TKE_CA)
     endif
 
@@ -1622,7 +1622,7 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_F
       if (GV%Boussinesq .or. GV%semi_Boussinesq .or. .not.(associated(fluxes%tau_mag))) then
         wind_TKE_src = CS%mstar*(GV%Z_to_H*U_star*U_Star*U_Star) * diag_wt
       else
-        wind_TKE_src = CS%mstar*(GV%RZ_to_H * US%L_to_Z*fluxes%tau_mag(i,j) * U_star) * diag_wt
+        wind_TKE_src = CS%mstar*(GV%RZ_to_H * fluxes%tau_mag(i,j) * U_star) * diag_wt
       endif
       CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + &
           ( wind_TKE_src + TKE_river(i) * diag_wt )
diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90
index 693485292e..c94e1032fe 100644
--- a/src/parameterizations/vertical/MOM_energetic_PBL.F90
+++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90
@@ -410,8 +410,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV,
   real :: u_star_BBL ! The bottom boundary layer friction velocity [H T-1 ~> m s-1 or kg m-2 s-1].
   real :: BBL_TKE   ! The mechanically generated turbulent kinetic energy available for bottom
                     ! boundary layer mixing within a timestep [R Z3 T-2 ~> J m-2]
-  real :: I_rho     ! The inverse of the Boussinesq reference density times a ratio of scaling
-                    ! factors [Z L-1 R-1 ~> m3 kg-1]
+  real :: I_rho     ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]
   real :: I_dt      ! The Adcroft reciprocal of the timestep [T-1 ~> s-1]
   real :: I_rho0dt  ! The inverse of the Boussinesq reference density times the time
                     ! step [R-1 T-1 ~> m3 kg-1 s-1]
@@ -493,7 +492,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV,
 
 
   h_neglect = GV%H_subroundoff
-  I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 ! This is not used when fully non-Boussinesq.
+  I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0 ! This is not used when fully non-Boussinesq.
   I_dt = 0.0 ; if (dt > 0.0) I_dt = 1.0 / dt
   I_rho0dt = 1.0 / (GV%Rho0 * dt)  ! This is not used when fully non-Boussinesq.
   BBL_mixing = ((CS%ePBL_BBL_effic > 0.0) .or. (CS%ePBL_tidal_effic > 0.0))
@@ -604,14 +603,14 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV,
         u_star_Mean = fluxes%ustar_gustless(i,j)
         mech_TKE = dt * GV%Rho0 * u_star**3
       elseif (allocated(tv%SpV_avg)) then
-        u_star = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1))
-        u_star_Mean = sqrt(US%L_to_Z*fluxes%tau_mag_gustless(i,j) * tv%SpV_avg(i,j,1))
-        mech_TKE = dt * u_star * US%L_to_Z*fluxes%tau_mag(i,j)
+        u_star = sqrt(fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1))
+        u_star_Mean = sqrt(fluxes%tau_mag_gustless(i,j) * tv%SpV_avg(i,j,1))
+        mech_TKE = dt * u_star * fluxes%tau_mag(i,j)
       else
         u_star = sqrt(fluxes%tau_mag(i,j) * I_rho)
-        u_star_Mean = sqrt(US%L_to_Z*fluxes%tau_mag_gustless(i,j) * I_rho)
+        u_star_Mean = sqrt(fluxes%tau_mag_gustless(i,j) * I_rho)
         mech_TKE = dt * GV%Rho0 * u_star**3
-        ! The line above is equivalent to: mech_TKE = dt * u_star * US%L_to_Z*fluxes%tau_mag(i,j)
+        ! The line above is equivalent to: mech_TKE = dt * u_star * fluxes%tau_mag(i,j)
       endif
 
       if (allocated(tv%SpV_avg) .and. .not.GV%Boussinesq) then
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index 0c0891930b..138a87754d 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -1777,8 +1777,7 @@ subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_t
   real :: ustar_sq          ! ustar squared [Z2 T-2 ~> m2 s-2]
   real :: Kd_mlr            ! A diffusivity associated with mixed layer turbulence radiation
                             ! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
-  real :: I_rho             ! The inverse of the reference density times a ratio of scaling
-                            ! factors [Z L-1 R-1 ~> m3 kg-1]
+  real :: I_rho             ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]
   real :: C1_6              ! 1/6 [nondim]
   real :: Omega2            ! rotation rate squared [T-2 ~> s-2].
   real :: z1                ! layer thickness times I_decay [nondim]
@@ -1794,7 +1793,7 @@ subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_t
   C1_6      = 1.0 / 6.0
   kml       = GV%nkml
   dz_neglect = GV%dz_subroundoff
-  I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 ! This is not used when fully non-Boussinesq.
+  I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0 ! This is not used when fully non-Boussinesq.
 
   if (.not.CS%ML_radiation) return
 
@@ -1816,12 +1815,12 @@ subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_t
       ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2
       u_star_H = GV%Z_to_H * fluxes%ustar(i,j)
     elseif (allocated(tv%SpV_avg)) then
-      ustar_sq = max(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1), CS%ustar_min**2)
-      u_star_H = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1))
+      ustar_sq = max(fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1), CS%ustar_min**2)
+      u_star_H = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1))
     else  ! This semi-Boussinesq form is mathematically equivalent to the Boussinesq version above.
      ! Differs at roundoff:  ustar_sq = max(fluxes%tau_mag(i,j) * I_rho, CS%ustar_min**2)
       ustar_sq = max((sqrt(fluxes%tau_mag(i,j) * I_rho))**2, CS%ustar_min**2)
-      u_star_H = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * GV%Rho0)
+      u_star_H = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) * GV%Rho0)
     endif
     TKE_ml_flux(i) = (CS%mstar * CS%ML_rad_coeff) * (ustar_sq * u_star_H)
     I_decay_len2_TKE = CS%TKE_decay**2 * (f_sq / ustar_sq)
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 9fa32d6e90..8a9b93faab 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -1980,7 +1980,6 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i,
   real :: h_shear ! The distance over which shears occur [Z ~> m].
   real :: dhc     ! The distance between the center of adjacent layers [Z ~> m].
   real :: visc_ml ! The mixed layer viscosity [H Z T-1 ~> m2 s-1 or Pa s].
-  real :: tau_scale  ! A scaling factor for the interpolated wind stress magnitude [H R-1 L-1 ~> m3 kg-1 or nondim]
   real :: I_Hmix  ! The inverse of the mixed layer thickness [Z-1 ~> m-1].
   real :: a_ml    ! The layer coupling coefficient across an interface in
                   ! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].
@@ -2008,8 +2007,6 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i,
   nz = GV%ke
   h_neglect = GV%dZ_subroundoff
 
-  tau_scale = US%L_to_Z * GV%RZ_to_H
-
   if (CS%answer_date < 20190101) then
     !   The maximum coupling coefficient was originally introduced to avoid
     ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90
index 99cc7b88c5..d9d46f7d6e 100644
--- a/src/user/Idealized_Hurricane.F90
+++ b/src/user/Idealized_Hurricane.F90
@@ -51,7 +51,7 @@ module Idealized_hurricane
   real    :: max_windspeed        !< Maximum wind speeds [L T-1 ~> m s-1]
   real    :: hurr_translation_spd !< Hurricane translation speed [L T-1 ~> m s-1]
   real    :: hurr_translation_dir !< Hurricane translation direction [radians]
-  real    :: gustiness            !< Gustiness (optional, used in u*) [R L Z T-2 ~> Pa]
+  real    :: gustiness            !< Gustiness (used in u*) [R Z2 T-2 ~> Pa]
   real    :: Rho0                 !< A reference ocean density [R ~> kg m-3]
   real    :: Hurr_cen_Y0          !< The initial y position of the hurricane
                                   !!  This experiment is conducted in a Cartesian
@@ -313,7 +313,7 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS)
                  units="kg m-3", default=1035.0, scale=US%kg_m3_to_R, do_not_log=.true.)
   call get_param(param_file, mdl, "GUST_CONST", CS%gustiness, &
                  "The background gustiness in the winds.", &
-                 units="Pa", default=0.0, scale=US%kg_m2s_to_RZ_T*US%m_s_to_L_T, do_not_log=.true.)
+                 units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z, do_not_log=.true.)
 
   if (CS%rad_edge >= CS%rad_ambient) call MOM_error(FATAL, &
     "idealized_hurricane_wind_init: IDL_HURR_RAD_AMBIENT must be larger than IDL_HURR_RAD_EDGE.")
@@ -437,16 +437,16 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS)
   !> Get Ustar
   if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
     !  This expression can be changed if desired, but need not be.
-    forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(US%L_to_Z * (CS%gustiness/CS%Rho0 + &
-            sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + &
-                 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))/CS%Rho0))
+    forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(CS%gustiness/CS%Rho0 + &
+            US%L_to_Z * sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + &
+                             0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))/CS%Rho0)
   enddo ; enddo ; endif
 
-  !> Get tau_mag [R L Z T-2 ~> Pa]
+  !> Get tau_mag [R Z2 T-2 ~> Pa]
   if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
     forces%tau_mag(i,j) = G%mask2dT(i,j) * (CS%gustiness + &
-            sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + &
-                 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2))))
+            US%L_to_Z * sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + &
+                             0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2))))
   enddo ; enddo ; endif
 
 end subroutine idealized_hurricane_wind_forcing
diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90
index 46cf6423d4..def4c59568 100644
--- a/src/user/SCM_CVMix_tests.F90
+++ b/src/user/SCM_CVMix_tests.F90
@@ -202,7 +202,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  ! The magnitude of the wind stress [R L Z T-2 ~> Pa]
+  real    :: mag_tau  ! The magnitude of the wind stress [R Z2 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
@@ -217,9 +217,9 @@ subroutine SCM_CVMix_tests_wind_forcing(sfc_state, forces, day, G, US, CS)
   enddo ; enddo
   call pass_vector(forces%taux, forces%tauy, G%Domain, To_All)
 
-  mag_tau = sqrt((CS%tau_x*CS%tau_x) + (CS%tau_y*CS%tau_y))
+  mag_tau = US%L_to_Z * sqrt((CS%tau_x*CS%tau_x) + (CS%tau_y*CS%tau_y))
   if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
-    forces%ustar(i,j) = sqrt( US%L_to_Z * mag_tau / CS%Rho0 )
+    forces%ustar(i,j) = sqrt( mag_tau / CS%Rho0 )
   enddo ; enddo ; endif
 
   if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie