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 f70cd34012..251f37290d 100644
--- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
+++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
@@ -624,13 +624,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   endif
 
   ! Set the wind stresses and ustar.
-  if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless)) then
+  if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless) .and. associated(fluxes%tau_mag)) then
     call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar, &
-                              gustless_ustar=fluxes%ustar_gustless)
-  elseif (associated(fluxes%ustar)) then
-    call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar)
-  elseif (associated(fluxes%ustar_gustless)) then
-    call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, gustless_ustar=fluxes%ustar_gustless)
+                              mag_tau=fluxes%tau_mag, gustless_ustar=fluxes%ustar_gustless)
+  else
+    if (associated(fluxes%ustar)) &
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar)
+    if (associated(fluxes%ustar_gustless)) &
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, gustless_ustar=fluxes%ustar_gustless)
+    if (associated(fluxes%tau_mag)) &
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, mag_tau=fluxes%tau_mag)
   endif
 
   if (coupler_type_initialized(fluxes%tr_fluxes) .and. &
@@ -674,7 +677,8 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   real, dimension(SZI_(G),SZJ_(G)) :: &
     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].
+    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]
 
   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]
@@ -778,12 +782,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   ! Set the wind stresses and ustar.
   if (wt1 <= 0.0) then
     call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux=forces%taux, tauy=forces%tauy, &
-                              ustar=forces%ustar, tau_halo=1)
+                              ustar=forces%ustar, mag_tau=forces%tau_mag, tau_halo=1)
   else
     call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux=forces%taux, tauy=forces%tauy, &
-                              ustar=ustar_tmp, tau_halo=1)
+                              ustar=ustar_tmp, mag_tau=tau_mag_tmp, tau_halo=1)
     do j=js,je ; do i=is,ie
       forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j)
+      forces%tau_mag(i,j) = wt1*forces%tau_mag(i,j) + wt2*tau_mag_tmp(i,j)
     enddo ; enddo
   endif
 
@@ -880,7 +885,7 @@ end subroutine convert_IOB_to_forces
 !! Ice_ocean_boundary_type into optional argument arrays, including changes of units, sign
 !! conventions, and putting the fields into arrays with MOM-standard sized halos.
 subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, ustar, &
-                                gustless_ustar, tau_halo)
+                                gustless_ustar, mag_tau, tau_halo)
   type(ice_ocean_boundary_type), &
                    target, intent(in)    :: IOB  !< An ice-ocean boundary type with fluxes to drive
                                                  !! the ocean in a coupled model
@@ -900,6 +905,9 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   real, dimension(SZI_(G),SZJ_(G)), &
                  optional, intent(out)   :: gustless_ustar !< The surface friction velocity without
                                                  !! 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]
   integer,       optional, intent(in)    :: tau_halo !< The halo size of wind stresses to set, 0 by default.
 
   ! Local variables
@@ -916,7 +924,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   real :: tau_mag       ! magnitude of the wind stress [R Z L T-2 ~> Pa]
   real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1]
 
-  logical :: do_ustar, do_gustless
+  logical :: do_ustar, do_gustless, do_tau_mag
   integer :: wind_stagger  ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains)
   integer :: i, j, is, ie, js, je, ish, ieh, jsh, jeh, Isqh, Ieqh, Jsqh, Jeqh, i0, j0, halo
 
@@ -929,7 +937,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   IRho0 = US%L_to_Z / CS%Rho0
   stress_conversion = US%Pa_to_RLZ_T2 * CS%wind_stress_multiplier
 
-  do_ustar = present(ustar) ; do_gustless = present(gustless_ustar)
+  do_ustar = present(ustar) ; do_gustless = present(gustless_ustar) ; do_tau_mag = present(mag_tau)
 
   wind_stagger = CS%wind_stagger
   if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. &
@@ -1022,13 +1030,13 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
     endif   ! endif for extracting wind stress fields with various staggerings
   endif
 
-  if (do_ustar .or. do_gustless) then
+  if (do_ustar .or. do_tau_mag .or. do_gustless) then
     ! Set surface friction velocity directly or as a function of staggering.
     ! ustar is required for the bulk mixed layer formulation and other turbulent mixing
     ! 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
-      if (do_ustar) then ; do j=js,je ; do i=is,ie
+      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
           if ((wind_stagger == CGRID_NE) .or. &
@@ -1038,7 +1046,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
                 (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0)) ) &
             gustiness = CS%gust(i,j)
         endif
-        ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0))
+        if (do_tau_mag) &
+          mag_tau(i,j) = gustiness + US%Pa_to_RLZ_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))
       enddo ; enddo ; endif
       if (CS%answer_date < 20190101) then
         if (do_gustless) then ; do j=js,je ; do i=is,ie
@@ -1062,6 +1073,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
           if (CS%read_gust_2d) gustiness = CS%gust(i,j)
         endif
         if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag)
+        if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
         if (CS%answer_date < 20190101) then
           if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
         else
@@ -1074,6 +1086,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         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 (CS%answer_date < 20190101) then
           if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
         else
@@ -1095,6 +1108,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         if (CS%read_gust_2d) 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 (CS%answer_date < 20190101) then
           if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
         else
diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
index 9b858af94e..ec5dab57a7 100644
--- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
+++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
@@ -774,6 +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%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag)
     enddo ; enddo
 
@@ -799,6 +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%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
@@ -820,8 +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%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%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 b8162b4b59..6d65ae4d28 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
@@ -851,6 +851,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%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag)
     enddo ; enddo
     call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1)
@@ -876,6 +877,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%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
@@ -897,8 +899,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%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2))
       else
+        forces%tau_mag(i,j) = CS%gust(i,j) + 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 0e8aedb8d0..859bfd81c8 100644
--- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90
@@ -407,10 +407,16 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS)
     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 )
     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 )
     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
+    enddo ; enddo ; endif
   endif
 
   call callTree_leave("wind_forcing_const")
@@ -522,9 +528,11 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)
   ! set the friction velocity
   if (CS%answer_date < 20190101) 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%ustar(i,j) = sqrt(US%L_to_Z * ((CS%gust_const/CS%Rho0) + &
-              sqrt(0.5*(forces%tauy(i,j-1)*forces%tauy(i,j-1) + forces%tauy(i,j)*forces%tauy(i,j) + &
-                        forces%taux(i-1,j)*forces%taux(i-1,j) + forces%taux(i,j)*forces%taux(i,j)))/CS%Rho0) )
+              sqrt(0.5*(forces%tauy(i,J-1)*forces%tauy(i,J-1) + forces%tauy(i,J)*forces%tauy(i,J) + &
+                        forces%taux(I-1,j)*forces%taux(I-1,j) + forces%taux(I,j)*forces%taux(I,j)))/CS%Rho0) )
     enddo ; enddo
   else
     call stresses_to_ustar(forces, G, US, CS)
@@ -725,11 +733,12 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
       if (.not.read_Ustar) then
         if (CS%read_gust_2d) then
           do j=js,je ; do i=is,ie
-            forces%ustar(i,j) = sqrt((CS%gust(i,j) + &
-                    sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j))) * US%L_to_Z / CS%Rho0)
+            forces%tau_mag(i,j) = CS%gust(i,j) + sqrt(temp_x(i,j)**2 + temp_y(i,j)**2)
+            forces%ustar(i,j) = sqrt(forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0)
           enddo ; enddo
         else
           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%ustar(i,j) = sqrt(US%L_to_Z * (CS%gust_const/CS%Rho0 + &
                     sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j)) / CS%Rho0) )
           enddo ; enddo
@@ -772,15 +781,19 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
       if (.not.read_Ustar) then
         if (CS%read_gust_2d) then
           do j=js,je ; do i=is,ie
-            forces%ustar(i,j) = sqrt((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 / CS%Rho0 )
+            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)))
+            forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 )
           enddo ; enddo
         else
           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%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))
+                    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
@@ -792,6 +805,9 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
     if (read_Ustar) then
       call MOM_read_data(filename, CS%Ustar_var, forces%ustar(:,:), &
                          G%Domain, timelevel=time_lev, scale=US%m_to_Z*US%T_to_s)
+      do j=js,je ; do i=is,ie
+        forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * forces%ustar(i,j)**2
+      enddo ; enddo
     endif
 
     CS%wind_last_lev = time_lev
@@ -815,6 +831,7 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS)
   ! Local variables
   real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal wind stresses at h-points [R Z L T-2 ~> Pa].
   real :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R Z L T-2 ~> Pa].
+  real :: ustar_tmp(SZI_(G),SZJ_(G)) ! The pre-override value of ustar [Z T-1 ~> m s-1]
   integer :: i, j
 
   call callTree_enter("wind_forcing_by_data_override, MOM_surface_forcing.F90")
@@ -840,17 +857,25 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS)
   if (CS%read_gust_2d) then
     call data_override(G%Domain, 'gust', CS%gust, day, scale=US%Pa_to_RLZ_T2)
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      forces%ustar(i,j) = sqrt((sqrt(temp_x(i,j)**2 + temp_y(i,j)**2) + &
-                               CS%gust(i,j)) * US%L_to_Z / CS%Rho0)
+      forces%tau_mag(i,j) = sqrt(temp_x(i,j)**2 + temp_y(i,j)**2) + CS%gust(i,j)
+      forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 )
     enddo ; enddo
   else
     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
+      ! forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 )
       forces%ustar(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
+
   ! Give the data override the option to modify the newly calculated forces%ustar.
+  ustar_tmp(:,:) = forces%ustar(:,:)
   call data_override(G%Domain, 'ustar', forces%ustar, day, scale=US%m_to_Z*US%T_to_s)
+  ! Only reset values where data override of ustar has occurred
+  do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (ustar_tmp(i,j) /= forces%ustar(i,j)) then
+    forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * forces%ustar(i,j)**2
+  endif ; enddo ; enddo
 
   call pass_vector(forces%taux, forces%tauy, G%Domain, To_All)
 
@@ -875,15 +900,17 @@ subroutine stresses_to_ustar(forces, G, US, CS)
 
   if (CS%read_gust_2d) then
     do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = sqrt( (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))) ) * I_rho )
+      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)))
+      forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * I_rho )
     enddo ; enddo
   else
     do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = sqrt( (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))) ) * I_rho )
+      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%ustar(i,j) = sqrt( forces%tau_mag(i,j) * I_rho )
     enddo ; enddo
   endif
 
diff --git a/config_src/drivers/solo_driver/user_surface_forcing.F90 b/config_src/drivers/solo_driver/user_surface_forcing.F90
index 42e732bb73..d7d3b89a8a 100644
--- a/config_src/drivers/solo_driver/user_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/user_surface_forcing.F90
@@ -88,9 +88,10 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS)
   !  is always positive.
   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((CS%gust_const + &
+    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/CS%Rho0))
+                 0.5*(forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2)))
+    forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(forces%tau_mag(i,j) * (US%L_to_Z/CS%Rho0))
   enddo ; enddo ; endif
 
 end subroutine USER_wind_forcing
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index af8481fd1c..7b9f2f9d3f 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -651,6 +651,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
     call create_group_pass(pass_tau_ustar_psurf, forces%taux, forces%tauy, G%Domain)
     if (associated(forces%ustar)) &
       call create_group_pass(pass_tau_ustar_psurf, forces%ustar, G%Domain)
+    if (associated(forces%tau_mag)) &
+      call create_group_pass(pass_tau_ustar_psurf, forces%tau_mag, G%Domain)
     if (associated(forces%p_surf)) &
       call create_group_pass(pass_tau_ustar_psurf, forces%p_surf, G%Domain)
     if (G%nonblocking_updates) then
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index 9a5e1f48f5..a59c33d525 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -68,6 +68,9 @@ module MOM_forcing_type
   ! surface stress components and turbulent velocity scale
   real, pointer, dimension(:,:) :: &
     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]
     ustar_gustless => NULL()   !< surface friction velocity scale without any
                                !! any augmentation for gustiness [Z T-1 ~> m s-1].
 
@@ -222,6 +225,8 @@ module MOM_forcing_type
   real, pointer, dimension(:,:) :: &
     taux  => NULL(), & !< zonal wind stress [R L Z T-2 ~> Pa]
     tauy  => NULL(), & !< meridional wind stress [R L Z T-2 ~> Pa]
+    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]
     ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1].
     net_mass_src => NULL() !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1]
 
@@ -359,6 +364,7 @@ module MOM_forcing_type
   integer :: id_taux  = -1
   integer :: id_tauy  = -1
   integer :: id_ustar = -1
+  integer :: id_tau_mag = -1
 
   integer :: id_psurf     = -1
   integer :: id_TKE_tidal = -1
@@ -1083,6 +1089,8 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
   ! and js...je as their extent.
   if (associated(fluxes%ustar)) &
     call hchksum(fluxes%ustar, mesg//" fluxes%ustar", G%HI, haloshift=hshift, scale=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, scale=US%RLZ_T2_to_Pa)
   if (associated(fluxes%buoy)) &
     call hchksum(fluxes%buoy, mesg//" fluxes%buoy ", G%HI, haloshift=hshift, scale=US%L_to_m**2*US%s_to_T**3)
   if (associated(fluxes%sw)) &
@@ -1186,11 +1194,13 @@ subroutine MOM_mech_forcing_chksum(mesg, forces, G, US, haloshift)
   ! and js...je as their extent.
   if (associated(forces%taux) .and. associated(forces%tauy)) &
     call uvchksum(mesg//" forces%tau[xy]", forces%taux, forces%tauy, G%HI, &
-                  haloshift=hshift, symmetric=.true., scale=US%RZ_T_to_kg_m2s*US%L_T_to_m_s)
+                  haloshift=hshift, symmetric=.true., scale=US%RLZ_T2_to_Pa)
   if (associated(forces%p_surf)) &
     call hchksum(forces%p_surf, mesg//" forces%p_surf", G%HI, haloshift=hshift, scale=US%RL2_T2_to_Pa)
   if (associated(forces%ustar)) &
     call hchksum(forces%ustar, mesg//" forces%ustar", G%HI, haloshift=hshift, scale=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, scale=US%RLZ_T2_to_Pa)
   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., &
@@ -1237,6 +1247,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg)
   write(0,'(2a)') 'MOM_forcing_type, forcing_SinglePointPrint: Called from ',mesg
   write(0,'(a,2es15.3)') 'MOM_forcing_type, forcing_SinglePointPrint: lon,lat = ',G%geoLonT(i,j),G%geoLatT(i,j)
   call locMsg(fluxes%ustar,'ustar')
+  call locMsg(fluxes%tau_mag,'tau_mag')
   call locMsg(fluxes%buoy,'buoy')
   call locMsg(fluxes%sw,'sw')
   call locMsg(fluxes%sw_vis_dir,'sw_vis_dir')
@@ -1305,18 +1316,22 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
 
   handles%id_taux = register_diag_field('ocean_model', 'taux', diag%axesCu1, Time,  &
         'Zonal surface stress from ocean interactions with atmos and ice', &
-        'Pa', conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s, &
+        'Pa', conversion=US%RLZ_T2_to_Pa, &
         standard_name='surface_downward_x_stress', cmor_field_name='tauuo',         &
         cmor_units='N m-2', cmor_long_name='Surface Downward X Stress',             &
         cmor_standard_name='surface_downward_x_stress')
 
   handles%id_tauy = register_diag_field('ocean_model', 'tauy', diag%axesCv1, Time,  &
         'Meridional surface stress ocean interactions with atmos and ice', &
-        'Pa', conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s, &
+        'Pa', conversion=US%RLZ_T2_to_Pa, &
         standard_name='surface_downward_y_stress', cmor_field_name='tauvo',        &
         cmor_units='N m-2', cmor_long_name='Surface Downward Y Stress',            &
         cmor_standard_name='surface_downward_y_stress')
 
+  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)
+
   handles%id_ustar = register_diag_field('ocean_model', 'ustar', diag%axesT1, Time, &
       'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', &
       'm s-1', conversion=US%Z_to_m*US%s_to_T)
@@ -2046,6 +2061,7 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces)
       fluxes%p_surf_full(i,j) = forces%p_surf_full(i,j)
 
       fluxes%ustar(i,j) = wt1*fluxes%ustar(i,j) + wt2*forces%ustar(i,j)
+      fluxes%tau_mag(i,j) = wt1*fluxes%tau_mag(i,j) + wt2*forces%tau_mag(i,j)
     enddo ; enddo
   else
     do j=js,je ; do i=is,ie
@@ -2053,6 +2069,7 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces)
       fluxes%p_surf_full(i,j) = flux_tmp%p_surf_full(i,j)
 
       fluxes%ustar(i,j) = wt1*fluxes%ustar(i,j) + wt2*flux_tmp%ustar(i,j)
+      fluxes%tau_mag(i,j) = wt1*fluxes%tau_mag(i,j) + wt2*flux_tmp%tau_mag(i,j)
     enddo ; enddo
   endif
 
@@ -2173,6 +2190,12 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres)
     enddo ; enddo
   endif
 
+  if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then
+    do j=js,je ; do i=is,ie
+      fluxes%tau_mag(i,j) = forces%tau_mag(i,j)
+    enddo ; enddo
+  endif
+
   if (do_pres) then
     if (associated(forces%p_surf) .and. associated(fluxes%p_surf)) then
       do j=js,je ; do i=is,ie
@@ -2304,6 +2327,12 @@ subroutine copy_back_forcing_fields(fluxes, forces, G)
     enddo ; enddo
   endif
 
+  if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then
+    do j=js,je ; do i=is,ie
+      forces%tau_mag(i,j) = fluxes%tau_mag(i,j)
+    enddo ; enddo
+  endif
+
 end subroutine copy_back_forcing_fields
 
 !> Offer mechanical forcing fields for diagnostics for those
@@ -2946,6 +2975,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
     if ((handles%id_buoy > 0) .and. associated(fluxes%buoy))                         &
       call post_data(handles%id_buoy, fluxes%buoy, diag)
 
+    if ((handles%id_tau_mag > 0) .and. associated(fluxes%tau_mag)) &
+      call post_data(handles%id_tau_mag, fluxes%tau_mag, diag)
+
     if ((handles%id_ustar > 0) .and. associated(fluxes%ustar)) &
       call post_data(handles%id_ustar, fluxes%ustar, diag)
 
@@ -3015,6 +3047,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
 
   call myAlloc(fluxes%ustar,isd,ied,jsd,jed, ustar)
   call myAlloc(fluxes%ustar_gustless,isd,ied,jsd,jed, ustar)
+  call myAlloc(fluxes%tau_mag,isd,ied,jsd,jed, ustar)
 
   call myAlloc(fluxes%evap,isd,ied,jsd,jed, water)
   call myAlloc(fluxes%lprec,isd,ied,jsd,jed, water)
@@ -3150,6 +3183,7 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, &
   call myAlloc(forces%tauy,isd,ied,JsdB,JedB, stress)
 
   call myAlloc(forces%ustar,isd,ied,jsd,jed, ustar)
+  call myAlloc(forces%tau_mag,isd,ied,jsd,jed, ustar)
 
   call myAlloc(forces%p_surf,isd,ied,jsd,jed, press)
   call myAlloc(forces%p_surf_full,isd,ied,jsd,jed, press)
@@ -3218,8 +3252,7 @@ subroutine get_forcing_groups(fluxes, water, heat, ustar, press, shelf, &
   !   to some degree.  But since this would be enforced at the driver level,
   !   we handle them here as independent flags.
 
-  ustar = associated(fluxes%ustar) &
-      .and. associated(fluxes%ustar_gustless)
+  ustar = associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless)
   ! TODO: Check for all associated fields, but for now just check one as a marker
   water = associated(fluxes%evap)
   heat = associated(fluxes%seaice_melt_heat)
@@ -3276,6 +3309,7 @@ subroutine deallocate_forcing_type(fluxes)
 
   if (associated(fluxes%ustar))                deallocate(fluxes%ustar)
   if (associated(fluxes%ustar_gustless))       deallocate(fluxes%ustar_gustless)
+  if (associated(fluxes%tau_mag))              deallocate(fluxes%tau_mag)
   if (associated(fluxes%buoy))                 deallocate(fluxes%buoy)
   if (associated(fluxes%sw))                   deallocate(fluxes%sw)
   if (associated(fluxes%seaice_melt_heat))     deallocate(fluxes%seaice_melt_heat)
@@ -3334,9 +3368,10 @@ end subroutine deallocate_forcing_type
 subroutine deallocate_mech_forcing(forces)
   type(mech_forcing), intent(inout) :: forces  !< Forcing fields structure
 
-  if (associated(forces%taux))  deallocate(forces%taux)
-  if (associated(forces%tauy))  deallocate(forces%tauy)
-  if (associated(forces%ustar)) deallocate(forces%ustar)
+  if (associated(forces%taux))           deallocate(forces%taux)
+  if (associated(forces%tauy))           deallocate(forces%tauy)
+  if (associated(forces%ustar))          deallocate(forces%ustar)
+  if (associated(forces%tau_mag))        deallocate(forces%tau_mag)
   if (associated(forces%p_surf))         deallocate(forces%p_surf)
   if (associated(forces%p_surf_full))    deallocate(forces%p_surf_full)
   if (associated(forces%net_mass_src))   deallocate(forces%net_mass_src)
@@ -3365,6 +3400,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns)
   if (do_ustar) then
     call rotate_array(fluxes_in%ustar, turns, fluxes%ustar)
     call rotate_array(fluxes_in%ustar_gustless, turns, fluxes%ustar_gustless)
+    call rotate_array(fluxes_in%tau_mag, turns, fluxes%tau_mag)
   endif
 
   if (do_water) then
@@ -3495,8 +3531,10 @@ subroutine rotate_mech_forcing(forces_in, turns, forces)
     call rotate_vector(forces_in%taux, forces_in%tauy, turns, &
         forces%taux, forces%tauy)
 
-  if (do_ustar) &
+  if (do_ustar) then
     call rotate_array(forces_in%ustar, turns, forces%ustar)
+    call rotate_array(forces_in%tau_mag, turns, forces%tau_mag)
+  endif
 
   if (do_shelf) then
     call rotate_array_pair( &
@@ -3555,24 +3593,27 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar)
                               do_press, do_iceberg)
 
   if (do_stress) then
-    tx_mean = global_area_mean_u(forces%taux, G, tmp_scale=US%Z_to_L*US%RL2_T2_to_Pa)
+    tx_mean = global_area_mean_u(forces%taux, G, tmp_scale=US%RLZ_T2_to_Pa)
     do j=js,je ; do i=isB,ieB
       if (G%mask2dCu(I,j) > 0.0) forces%taux(I,j) = tx_mean
     enddo ; enddo
-    ty_mean = global_area_mean_v(forces%tauy, G, tmp_scale=US%Z_to_L*US%RL2_T2_to_Pa)
+    ty_mean = global_area_mean_v(forces%tauy, G, tmp_scale=US%RLZ_T2_to_Pa)
     do j=jsB,jeB ; do i=is,ie
       if (G%mask2dCv(i,J) > 0.0) forces%tauy(i,J) = ty_mean
     enddo ; enddo
     if (tau2ustar) then
-      do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0.0) forces%ustar(i,j) = sqrt(sqrt(tx_mean**2 + ty_mean**2)*Irho0)
-      enddo ; enddo
+      do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
+        forces%tau_mag(i,j) = sqrt(tx_mean**2 + ty_mean**2)
+        forces%ustar(i,j) = sqrt(forces%tau_mag(i,j) * Irho0)
+      endif ; enddo ; enddo
     else
       call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T)
+      call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa)
     endif
   else
     if (do_ustar) then
       call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T)
+      call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa)
     endif
   endif
 
@@ -3613,6 +3654,7 @@ subroutine homogenize_forcing(fluxes, G, GV, US)
   if (do_ustar) then
     call homogenize_field_t(fluxes%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T)
     call homogenize_field_t(fluxes%ustar_gustless, G, tmp_scale=US%Z_to_m*US%s_to_T)
+    call homogenize_field_t(fluxes%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa)
   endif
 
   if (do_water) then
diff --git a/src/framework/MOM_unit_scaling.F90 b/src/framework/MOM_unit_scaling.F90
index 482c2eec7a..868352102e 100644
--- a/src/framework/MOM_unit_scaling.F90
+++ b/src/framework/MOM_unit_scaling.F90
@@ -52,6 +52,7 @@ module MOM_unit_scaling
   real :: RZ3_T3_to_W_m2  !< Convert turbulent kinetic energy fluxes from R Z3 T-3 to W m-2    [W T3 R-1 Z-3 m-2 ~> 1]
   real :: W_m2_to_RZ3_T3  !< Convert turbulent kinetic energy fluxes from W m-2 to R Z3 T-3     [R Z3 m2 T-3 W-1 ~> 1]
   real :: RL2_T2_to_Pa    !< Convert pressures from R L2 T-2 to Pa                                [Pa T2 R-1 L-2 ~> 1]
+  real :: RLZ_T2_to_Pa    !< Convert wind stresses from R L Z T-2 to Pa                       [Pa T2 R-1 L-1 Z-1 ~> 1]
   real :: Pa_to_RL2_T2    !< Convert pressures from Pa to R L2 T-2                                [R L2 T-2 Pa-1 ~> 1]
   real :: Pa_to_RLZ_T2    !< Convert wind stresses from Pa to R L Z T-2                          [R L Z T-2 Pa-1 ~> 1]
 
@@ -221,6 +222,7 @@ subroutine set_unit_scaling_combos(US)
   US%RL2_T2_to_Pa = US%R_to_kg_m3 * US%L_T_to_m_s**2
   US%Pa_to_RL2_T2 = US%kg_m3_to_R * US%m_s_to_L_T**2
   ! Wind stresses:
+  US%RLZ_T2_to_Pa = US%R_to_kg_m3 * US%L_T_to_m_s**2 * US%Z_to_L
   US%Pa_to_RLZ_T2 = US%kg_m3_to_R * US%m_s_to_L_T**2 * US%L_to_Z
 
 end subroutine set_unit_scaling_combos