diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90
index 7b1b4b3946..104a2b0312 100644
--- a/src/user/SCM_CVMix_tests.F90
+++ b/src/user/SCM_CVMix_tests.F90
@@ -42,6 +42,8 @@ module SCM_CVMix_tests
   real :: surf_evap !< (Constant) Evaporation rate [Z T-1 ~> m s-1]
   real :: Max_sw !< maximum of diurnal sw radiation [C Z T-1 ~> degC m s-1]
   real :: Rho0 !< reference density [R ~> kg m-3]
+  real :: rho_restore !< The density that is used to convert piston velocities
+                      !! into salt or heat fluxes [R ~> kg m-3]
 end type
 
 ! This include declares and sets the variable "version".
@@ -184,6 +186,9 @@ subroutine SCM_CVMix_tests_surface_forcing_init(Time, G, param_file, CS)
                  "properties, or with BOUSSINSEQ false to convert some "//&
                  "parameters from vertical units of m to kg m-2.", &
                  units="kg m-3", default=1035.0, scale=US%kg_m3_to_R)
+  call get_param(param_file, mdl, "RESTORE_FLUX_RHO", CS%rho_restore, &
+                 "The density that is used to convert piston velocities into salt or heat fluxes.", &
+                 units="kg m-3", default=CS%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R)
 
 end subroutine SCM_CVMix_tests_surface_forcing_init
 
@@ -214,7 +219,11 @@ subroutine SCM_CVMix_tests_wind_forcing(sfc_state, forces, day, G, US, CS)
 
   mag_tau = 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( US%L_to_Z * mag_tau / 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
   enddo ; enddo ; endif
 
 end subroutine SCM_CVMix_tests_wind_forcing
@@ -246,7 +255,7 @@ subroutine SCM_CVMix_tests_buoyancy_forcing(sfc_state, fluxes, day, G, US, CS)
     ! therefore must convert to [Q R Z T-1 ~> W m-2] by multiplying
     ! by Rho0*Cp
     do J=Jsq,Jeq ; do i=is,ie
-      fluxes%sens(i,J) = CS%surf_HF * CS%Rho0 * fluxes%C_p
+      fluxes%sens(i,J) = CS%surf_HF * CS%rho_restore * fluxes%C_p
     enddo ; enddo
   endif
 
@@ -255,7 +264,7 @@ subroutine SCM_CVMix_tests_buoyancy_forcing(sfc_state, fluxes, day, G, US, CS)
     ! Note CVMix test inputs give evaporation in [Z T-1 ~> m s-1]
     ! This therefore must be converted to mass flux in [R Z T-1 ~> kg m-2 s-1]
     ! by multiplying by density and some unit conversion factors.
-      fluxes%evap(i,J) = CS%surf_evap * CS%Rho0
+      fluxes%evap(i,J) = CS%surf_evap * CS%rho_restore
     enddo ; enddo
   endif
 
@@ -264,7 +273,8 @@ subroutine SCM_CVMix_tests_buoyancy_forcing(sfc_state, fluxes, day, G, US, CS)
     ! Note CVMix test inputs give max sw rad in [Z C T-1 ~> m degC s-1]
     ! therefore must convert to [Q R Z T-1 ~> W m-2] by multiplying by Rho0*Cp
     ! Note diurnal cycle peaks at Noon.
-      fluxes%sw(i,J) = CS%Max_sw *  max(0.0, cos(2*PI*(time_type_to_real(DAY)/86400.0 - 0.5))) * CS%RHO0 * fluxes%C_p
+      fluxes%sw(i,J) = CS%Max_sw *  max(0.0, cos(2*PI*(time_type_to_real(DAY)/86400.0 - 0.5))) * &
+                       CS%rho_restore * fluxes%C_p
     enddo ; enddo
   endif
 
diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90
index ca383ba1f1..6f6e4da439 100644
--- a/src/user/dumbbell_surface_forcing.F90
+++ b/src/user/dumbbell_surface_forcing.F90
@@ -25,9 +25,8 @@ module dumbbell_surface_forcing
 type, public :: dumbbell_surface_forcing_CS ; private
   logical :: use_temperature !< If true, temperature and salinity are used as state variables.
   logical :: restorebuoy     !< If true, use restoring surface buoyancy forcing.
-  real :: Rho0               !< The density used in the Boussinesq approximation [R ~> kg m-3].
   real :: G_Earth            !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
-  real :: Flux_const         !< The restoring rate at the surface [Z T-1 ~> m s-1].
+  real :: Flux_const         !< The restoring rate at the surface [R Z T-1 ~> kg m-2 s-1].
 ! real :: gust_const         !< A constant unresolved background gustiness
 !                            !! that contributes to ustar [R L Z T-2 ~> Pa].
   real :: slp_amplitude      !< The amplitude of pressure loading [R L2 T-2 ~> Pa] applied
@@ -114,7 +113,7 @@ subroutine dumbbell_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
   if (CS%use_temperature .and. CS%restorebuoy) then
     do j=js,je ; do i=is,ie
       if (CS%forcing_mask(i,j)>0.) then
-        fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)) * &
+        fluxes%vprec(i,j) = - (G%mask2dT(i,j) * CS%Flux_const) * &
                 ((CS%S_restore(i,j) - sfc_state%SSS(i,j)) /  (0.5 * (CS%S_restore(i,j) + sfc_state%SSS(i,j))))
 
       endif
@@ -181,6 +180,9 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
   real :: S_surf  ! Initial surface salinity [S ~> ppt]
   real :: S_range ! Range of the initial vertical distribution of salinity [S ~> ppt]
   real :: x       ! Latitude normalized by the domain size [nondim]
+  real :: Rho0          ! The density used in the Boussinesq approximation [R ~> kg m-3]
+  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]
   integer :: i, j
   logical :: dbrotate    ! If true, rotate the domain.
 # include "version_variable.h"
@@ -202,7 +204,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
   call get_param(param_file, mdl, "G_EARTH", CS%G_Earth, &
                  "The gravitational acceleration of the Earth.", &
                  units="m s-2", default=9.80, scale=US%m_to_L**2*US%Z_to_m*US%T_to_s**2)
-  call get_param(param_file, mdl, "RHO_0", CS%Rho0, &
+  call get_param(param_file, mdl, "RHO_0", Rho0, &
                  "The mean ocean density used with BOUSSINESQ true to "//&
                  "calculate accelerations and the mass for conservation "//&
                  "properties, or with BOUSSINSEQ false to convert some "//&
@@ -233,8 +235,13 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
                  "The constant that relates the restoring surface fluxes to the relative "//&
                  "surface anomalies (akin to a piston velocity).  Note the non-MKS units.", &
                  default=0.0, units="m day-1", scale=US%m_to_Z*US%T_to_s)
-    ! Convert CS%Flux_const from m day-1 to m s-1.
-    CS%Flux_const = CS%Flux_const / 86400.0
+    call get_param(param_file, mdl, "RESTORE_FLUX_RHO", rho_restore, &
+                 "The density that is used to convert piston velocities into salt or heat "//&
+                 "fluxes with RESTORE_SALINITY or RESTORE_TEMPERATURE.", &
+                 units="kg m-3", default=Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R, &
+                 do_not_log=(CS%Flux_const==0.0))
+    ! Convert FLUXCONST from m day-1 to m s-1 and Flux_const to [R Z T-1 ~> kg m-2 s-1]
+    CS%Flux_const = rho_restore * (CS%Flux_const / 86400.0)
 
 
     allocate(CS%forcing_mask(G%isd:G%ied, G%jsd:G%jed), source=0.0)