From 50df270e87c40ad861c8d65cde6de9b53de6cf2f Mon Sep 17 00:00:00 2001
From: OlgaSergienko <39838355+OlgaSergienko@users.noreply.github.com>
Date: Fri, 17 Dec 2021 17:33:09 -0500
Subject: [PATCH 1/7] Ice dynamics (#35)

* In MOM_ice_shelf_dynamics.F90 changes are  made to calc_shelf_visc() and calc_shelf_driving_stress() to account for an irregular quadrilateral grid. In MOM_ice_shelf_initialize.F90 changes are made to initialize_ice_thickness_from_file() to correct masks at initialization.

* Changed indentation

* Changes are made to 'calc_shelf_visc()` to make computations of the velocity derivatives rotation-invariant. Changes in `update_ice_shelf()` utilize G%IareaT(:,:) instead of 1/G%areaT(:,:).
---
 src/ice_shelf/MOM_ice_shelf_dynamics.F90   | 118 ++++++++++++++-------
 src/ice_shelf/MOM_ice_shelf_initialize.F90 |  45 ++++----
 2 files changed, 99 insertions(+), 64 deletions(-)

diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
index df2e801613..8fb674e36c 100644
--- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90
+++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
@@ -260,9 +260,9 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS)
     allocate( CS%v_shelf(IsdB:IedB,JsdB:JedB), source=0.0 )
     allocate( CS%t_shelf(isd:ied,jsd:jed), source=-10.0 ) ! [degC]
     allocate( CS%ice_visc(isd:ied,jsd:jed), source=0.0 )
-    allocate( CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25 ) ! [Units?]
+    allocate( CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25 ) ! [Pa-3s-1]
     allocate( CS%basal_traction(isd:ied,jsd:jed), source=0.0 )
-    allocate( CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10 ) ! [Units?]
+    allocate( CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10 ) ! [Pa (m-1 s)^n_sliding]
     allocate( CS%OD_av(isd:ied,jsd:jed), source=0.0 )
     allocate( CS%ground_frac(isd:ied,jsd:jed), source=0.0 )
     allocate( CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0 )
@@ -553,8 +553,8 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
      call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE)
      call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE)
 
-     !initialize ice flow velocities from file
-     call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, &
+     !initialize ice flow characteristic (velocities, bed elevation under the grounded part, etc) from file
+     call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, &
             G, US, param_file)
      call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE)
      call pass_var(CS%bed_elev, G%domain,CENTER)
@@ -567,9 +567,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
        'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s)
     ! I think that the conversion factors for the next two diagnostics are wrong. - RWH
     CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesB1, Time, &
-       'x-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa)
+       'x-driving stress of ice', 'kPa', conversion=1.e-3*US%RL2_T2_to_Pa)
     CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesB1, Time, &
-       'y-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa)
+       'y-driving stress of ice', 'kPa', conversion=1.e-3*US%RL2_T2_to_Pa)
     CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesB1, Time, &
        'mask for u-nodes', 'none')
     CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesB1, Time, &
@@ -579,9 +579,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
     CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, &
        'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m)
     CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, &
-       'viscosity', 'm', conversion=1e-6*US%Z_to_m)
+       'vi-viscosity', 'Pa s-1 m', conversion=US%RL2_T2_to_Pa*US%L_T_to_m_s) !vertically integrated viscosity
     CS%id_taub = register_diag_field('ice_shelf_model','taub_beta',CS%diag%axesT1, Time, &
-       'taub', 'Pa yr m-1', conversion=1e-6*US%Z_to_m)
+       'taub', 'MPa', conversion=1e-6*US%RL2_T2_to_Pa)
     CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, &
        'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m)
   endif
@@ -673,7 +673,10 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled
   logical,      optional, intent(in)    :: coupled_grounding !< If true, the grounding line is
                                               !! determined by coupled ice-ocean dynamics
   logical,      optional, intent(in)    :: must_update_vel !< Always update the ice velocities if true.
-
+  real, dimension(SZDIB_(G),SZDJB_(G))  ::taud_x,taud_y   !<area-averaged driving stress [R L2T-2 ~> Pa]
+  real, dimension(SZDI_(G),SZDJ_(G))  :: ice_visc ! <area-averaged vertically integrated ice viscosity
+                                              !![R2 L2T-3 ~> Pa s-1 m]
+  real, dimension(SZDI_(G),SZDJ_(G))  :: basal_tr ! <area-averaged basal traction [R L2T-2 ~> Pa]
   integer :: iters
   logical :: update_ice_vel, coupled_GL
 
@@ -706,12 +709,24 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled
      if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag)
      if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag)
 !    if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag)
-     if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag)
-     if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag)
+     if (CS%id_taudx_shelf > 0) then
+        taud_x(:,:) = CS%taudx_shelf(:,:)*G%IareaT(:,:)
+        call post_data(CS%id_taudx_shelf,taud_x , CS%diag)
+     endif
+     if (CS%id_taudy_shelf > 0) then
+        taud_y(:,:) = CS%taudy_shelf(:,:)*G%IareaT(:,:)
+        call post_data(CS%id_taudy_shelf,taud_y , CS%diag)
+     endif
      if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag)
      if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag)
-     if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag)
-     if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag)
+     if (CS%id_visc_shelf > 0) then
+       ice_visc(:,:)=CS%ice_visc(:,:)*G%IareaT(:,:)
+       call post_data(CS%id_visc_shelf, ice_visc,CS%diag)
+     endif
+     if (CS%id_taub > 0) then
+        basal_tr(:,:) = CS%basal_traction(:,:)*G%IareaT(:,:)
+        call post_data(CS%id_taub, basal_tr,CS%diag)
+     endif
 !!
      if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag)
      if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag)
@@ -874,6 +889,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i
       if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) > 0) then
         float_cond(i,j) = 1.0
         CS%ground_frac(i,j) = 1.0
+        CS%OD_av(i,j) =0.0
       endif
     enddo
   enddo
@@ -960,7 +976,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i
 
   !! begin loop
 
-  do iter=1,100
+  do iter=1,50
 
     call ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H_node, float_cond, &
                                ISS%hmask, conv_flag, iters, time, Phi, Phisub)
@@ -1775,7 +1791,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
                          intent(inout) :: taudx  !< X-direction driving stress at q-points [kg L s-2 ~> kg m s-2]
   real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: taudy  !< Y-direction driving stress at q-points [kg L s-2 ~> kg m s-2]
-                                                  ! This will become [R L3 Z T-2 ~> kg m s-2]
+                                                  ! This will become  [R L3 Z T-2 ~> kg m s-2]
 
 ! driving stress!
 
@@ -1790,12 +1806,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
 
   real, dimension(SIZE(OD,1),SIZE(OD,2))  :: S, &     ! surface elevation [Z ~> m].
                             BASE     ! basal elevation of shelf/stream [Z ~> m].
+  real, pointer, dimension(:,:,:,:) :: Phi => NULL() ! The gradients of bilinear basis elements at Gaussian
+                                                ! quadrature points surrounding the cell vertices [m-1].
 
 
   real    :: rho, rhow, rhoi_rhow ! Ice and ocean densities [R ~> kg m-3]
   real    :: sx, sy    ! Ice shelf top slopes [Z L-1 ~> nondim]
   real    :: neumann_val ! [R Z L2 T-2 ~> kg s-2]
-  real    :: dxh, dyh  ! Local grid spacing [L ~> m]
+  real    :: dxh, dyh,Dx,Dy  ! Local grid spacing [L ~> m]
   real    :: grav      ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
 
   integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
@@ -1813,6 +1831,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
   is = iscq - 1; js = jscq - 1
   i_off = G%idg_offset ; j_off = G%jdg_offset
 
+
   rho =  CS%density_ice
   rhow = CS%density_ocean_avg
   grav = CS%g_Earth
@@ -1821,13 +1840,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
 
   ! or is this faster?
   BASE(:,:) = -CS%bed_elev(:,:) + OD(:,:)
-  S(:,:) = BASE(:,:) + ISS%h_shelf(:,:)
-
+  S(:,:) = -CS%bed_elev(:,:) + ISS%h_shelf(:,:)
   ! check whether the ice is floating or grounded
   do j=jsc-G%domain%njhalo,jec+G%domain%njhalo
     do i=isc-G%domain%nihalo,iec+G%domain%nihalo
       if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then
         S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j)
+      else
+        S(i,j)=ISS%h_shelf(i,j)-CS%bed_elev(i,j)
       endif
     enddo
   enddo
@@ -1838,7 +1858,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
       sy = 0
       dxh = G%dxT(i,j)
       dyh = G%dyT(i,j)
-
+      Dx=dxh
+      Dy=dyh
       if (ISS%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell
 
         ! calculate sx
@@ -1857,12 +1878,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
         else ! interior
           if (ISS%hmask(i+1,j) == 1) then
             cnt = cnt+1
+            Dx =dxh+ G%dxT(i+1,j)
             sx = S(i+1,j)
           else
             sx = S(i,j)
           endif
           if (ISS%hmask(i-1,j) == 1) then
             cnt = cnt+1
+            Dx =dxh+ G%dxT(i-1,j)
             sx = sx - S(i-1,j)
           else
             sx = sx - S(i,j)
@@ -1870,7 +1893,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
           if (cnt == 0) then
             sx = 0
           else
-            sx = sx / (cnt * dxh)
+            sx = sx / ( Dx)
           endif
         endif
 
@@ -1892,6 +1915,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
         else ! interior
           if (ISS%hmask(i,j+1) == 1) then
             cnt = cnt+1
+            Dy =dyh+ G%dyT(i,j+1)
             sy = S(i,j+1)
           else
             sy = S(i,j)
@@ -1899,13 +1923,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
           if (ISS%hmask(i,j-1) == 1) then
             cnt = cnt+1
             sy = sy - S(i,j-1)
+            Dy =dyh+ G%dyT(i,j-1)
           else
             sy = sy - S(i,j)
           endif
           if (cnt == 0) then
             sy = 0
           else
-            sy = sy / (cnt * dyh)
+            sy = sy / (Dy)
           endif
         endif
 
@@ -1930,10 +1955,10 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
                 taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j)
         endif
         if (CS%ground_frac(i,j) == 1) then
-!          neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2)
+!          neumann_val = (.5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2))
           neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2
         else
-          neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2
+          neumann_val = (.5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2)
         endif
 
         if ((CS%u_face_mask(I-1,j) == 2) .OR. (ISS%hmask(i-1,j) == 0) .OR. (ISS%hmask(i-1,j) == 2) ) then
@@ -1971,7 +1996,6 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
       endif
     enddo
   enddo
-
 end subroutine calc_shelf_driving_stress
 
 subroutine init_boundary_values(CS, G, time, hmask, input_flux, input_thick, new_sim)
@@ -2528,8 +2552,8 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc,
   call pass_vector(u_bdry_contr, v_bdry_contr, G%domain, TO_ALL, BGRID_NE)
 end subroutine apply_boundary_values
 
+
 !> Update depth integrated viscosity, based on horizontal strain rates, and also update the
-!! nonlinear part of the basal traction.
 subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
   type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
   type(ice_shelf_state),  intent(in)    :: ISS !< A structure with elements that describe
@@ -2540,9 +2564,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
                           intent(inout) :: u_shlf !< The zonal ice shelf velocity [L T-1 ~> m s-1].
   real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
                           intent(inout) :: v_shlf !< The meridional ice shelf velocity [L T-1 ~> m s-1].
-
+  real, pointer, dimension(:,:,:,:) :: Phi => NULL() ! The gradients of bilinear basis elements at Gaussian
+                                                ! quadrature points surrounding the cell vertices [m-1].
 ! update DEPTH_INTEGRATED viscosity, based on horizontal strain rates - this is for bilinear FEM solve
-! so there is an "upper" and "lower" bilinear viscosity
 
 ! also this subroutine updates the nonlinear part of the basal traction
 
@@ -2553,7 +2577,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
   real :: Visc_coef, n_g
   real :: ux, uy, vx, vy
   real :: eps_min, dxh, dyh ! Velocity shears [T-1 ~> s-1]
-  real, dimension(8,4)  :: Phi
+!  real, dimension(8,4)  :: Phi
   real, dimension(2) :: xquad
 !  real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1]
 
@@ -2566,6 +2590,12 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
   is = iscq - 1; js = jscq - 1
     i_off = G%idg_offset ; j_off = G%jdg_offset
 
+  allocate(Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0)
+
+  do j=jsc,jec ; do i=isc,iec
+    call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j))
+  enddo ; enddo
+
   n_g = CS%n_glen; eps_min = CS%eps_glen_min
   CS%ice_visc(:,:)=1e22
 !  Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen)
@@ -2575,21 +2605,35 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
       if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then
         Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%AGlen_visc(i,j))**(-1./CS%n_glen)
 
-        ux = ((u_shlf(I,J) + (u_shlf(I,J-1) + u_shlf(I,J+1))) - &
-                (u_shlf(I-1,J) + (u_shlf(I-1,J-1) + u_shlf(I-1,J+1)))) / (3*G%dxT(i,j))
-        vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - &
-                (v_shlf(I-1,J) + (v_shlf(I-1,J-1) + v_shlf(I-1,J+1)))) / (3*G%dxT(i,j))
-        uy = ((u_shlf(I,J) + (u_shlf(I-1,J) + u_shlf(I+1,J))) - &
-                (u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j))
-        vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - &
-                (v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j))
+      do iq=1,2 ; do jq=1,2
+
+        ux = ( (u_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq,i,j) + &
+                u_shlf(I,J) * Phi(7,2*(jq-1)+iq,i,j)) + &
+               (u_shlf(I-1,J) * Phi(5,2*(jq-1)+iq,i,j) + &
+                u_shlf(I,J-1) * Phi(3,2*(jq-1)+iq,i,j)) )
+
+        vx = ( (v_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq,i,j) + &
+               v_shlf(I,J) * Phi(7,2*(jq-1)+iq,i,j)) + &
+               (v_shlf(I-1,J) * Phi(5,2*(jq-1)+iq,i,j) + &
+               v_shlf(I,J-1) * Phi(3,2*(jq-1)+iq,i,j)) )
+
+        uy = ( (u_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq,i,j) + &
+               u_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j)) + &
+              (u_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + &
+               u_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j)) )
+
+        vy = ( (v_shlf(I-1,j-1) * Phi(2,2*(jq-1)+iq,i,j) + &
+               v_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j)) + &
+              (v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + &
+              v_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j)) )
+     enddo ; enddo
 !        CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging
         CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * &
              (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g))
       endif
     enddo
   enddo
-
+  deallocate(Phi)
 end subroutine calc_shelf_visc
 
 subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90
index 469cba39ce..7cc3c020a3 100644
--- a/src/ice_shelf/MOM_ice_shelf_initialize.F90
+++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90
@@ -149,14 +149,13 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
     enddo ; enddo
   endif
 
-  if (len_sidestress > 0.) then
     do j=jsc,jec
       do i=isc,iec
 
       ! taper ice shelf in area where there is no sidestress -
       ! but do not interfere with hmask
 
-        if (G%geoLonCv(i,j) > len_sidestress) then
+        if ((len_sidestress > 0.) .and. (G%geoLonCv(i,j) > len_sidestress)) then
           udh = exp(-(G%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j)
           if (udh <= 25.0) then
             h_shelf(i,j) = 0.0
@@ -180,7 +179,6 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
         endif
       enddo
     enddo
-  endif
 end subroutine initialize_ice_thickness_from_file
 
 !> Initialize ice shelf thickness for a channel configuration
@@ -397,13 +395,13 @@ end subroutine initialize_ice_shelf_boundary_channel
 
 
 !> Initialize ice shelf flow from file
-subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
-                                         hmask,h_shelf, G, US, PF)
-!subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,ice_visc,float_cond,&
+!subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
 !                                         hmask,h_shelf, G, US, PF)
+subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
+                                         G, US, PF)
   type(ocean_grid_type), intent(in)    :: G    !< The ocean's grid structure
   real, dimension(SZDI_(G),SZDJ_(G)), &
-                         intent(inout) :: bed_elev !< The ice shelf u velocity  [Z ~> m].
+                         intent(inout) :: bed_elev !< The bed elevation   [Z ~> m].
    real, dimension(SZIB_(G),SZJB_(G)), &
                           intent(inout) :: u_shelf !< The zonal ice shelf velocity  [L T-1 ~> m s-1].
    real, dimension(SZIB_(G),SZJB_(G)), &
@@ -412,12 +410,12 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
   real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(inout)    :: float_cond !< An array indicating where the ice
                                                 !! shelf is floating: 0 if floating, 1 if not.
-  real, dimension(SZDI_(G),SZDJ_(G)), &
-                         intent(in) :: hmask !< A mask indicating which tracer points are
-                                             !! partly or fully covered by an ice-shelf
-  real, dimension(SZDI_(G),SZDJ_(G)), &
-                         intent(in) :: h_shelf !< A mask indicating which tracer points are
+!  real, dimension(SZDI_(G),SZDJ_(G)), &
+!                         intent(in) :: hmask !< A mask indicating which tracer points are
                                              !! partly or fully covered by an ice-shelf
+!  real, dimension(SZDI_(G),SZDJ_(G)), &
+!                         intent(in) :: h_shelf !< A mask indicating which tracer points are
+!                                             !! partly or fully covered by an ice-shelf
   type(unit_scale_type), intent(in)    :: US !< A structure containing unit conversion factors
   type(param_file_type), intent(in)    :: PF !< A structure to parse for run-time parameters
 
@@ -453,10 +451,10 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
                  "The name of the thickness variable in ICE_VELOCITY_FILE.", &
                  default="viscosity")
   call get_param(PF, mdl, "BED_TOPO_FILE", bed_topo_file, &
-                 "The file from which the velocity is read.", &
+                 "The file from which the bed elevation is read.", &
                  default="ice_shelf_vel.nc")
   call get_param(PF, mdl, "BED_TOPO_VARNAME", bed_varname, &
-                 "The name of the thickness variable in ICE_VELOCITY_FILE.", &
+                 "The name of the thickness variable in ICE_INPUT_FILE.", &
                  default="depth")
   if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
        " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename))
@@ -470,15 +468,8 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
 
  filename = trim(inputdir)//trim(bed_topo_file)
  call MOM_read_data(filename,trim(bed_varname), bed_elev, G%Domain, scale=1.)
-  isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec
+!  isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec
 
-!  do j=jsc,jec
-!    do i=isc,iec
-!       if  (hmask(i,j) == 1.) then
-!               ice_visc(i,j) = ice_visc(i,j) * (G%areaT(i,j) * h_shelf(i,j))
-!       endif
-!    enddo
-!  enddo
 
 end subroutine initialize_ice_flow_from_file
 
@@ -510,7 +501,7 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask
                           intent(inout) :: hmask !< A mask indicating which tracer points are
                                               !! partly or fully covered by an ice-shelf
    real, dimension(SZDI_(G),SZDJ_(G)), &
-                          intent(inout) :: h_shelf !< Ice-shelf thickness
+                          intent(in) :: h_shelf !< Ice-shelf thickness
    type(unit_scale_type), intent(in)    :: US !< A structure containing unit conversion factors
    type(param_file_type), intent(in)    :: PF !< A structure to parse for run-time parameters
 
@@ -535,9 +526,9 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask
   call get_param(PF, mdl, "ICE_THICKNESS_FILE", icethick_file, &
                  "The file from which the ice-shelf thickness is read.", &
                  default="ice_shelf_thick.nc")
-  call get_param(PF, mdl, "ICE_THICKNESS_VARNAME", h_varname, &
-                 "The name of the thickness variable in ICE_THICKNESS_FILE.", &
-                 default="h_shelf")
+!  call get_param(PF, mdl, "ICE_THICKNESS_VARNAME", h_varname, &
+!                 "The name of the thickness variable in ICE_THICKNESS_FILE.", &
+!                 default="h_shelf")
   call get_param(PF, mdl, "ICE_THICKNESS_MASK_VARNAME", hmsk_varname, &
                  "The name of the icethickness mask variable in ICE_THICKNESS_FILE.", &
                  default="h_mask")
@@ -574,7 +565,7 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask
  call MOM_read_data(filename,trim(vmask_varname), vmask, G%Domain, position=CORNER,scale=1.)
  filename = trim(inputdir)//trim(icethick_file)
 
- call MOM_read_data(filename, trim(h_varname), h_shelf, G%Domain, scale=US%m_to_Z)
+! call MOM_read_data(filename, trim(h_varname), h_shelf, G%Domain, scale=US%m_to_Z)
  call MOM_read_data(filename,trim(hmsk_varname), hmask, G%Domain, scale=1.)
   isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec
 

From 12f29f67abdc6a8afb783e50159b1d9ffe014692 Mon Sep 17 00:00:00 2001
From: wfcooke <wfcooke@users.noreply.github.com>
Date: Fri, 17 Dec 2021 21:48:53 -0500
Subject: [PATCH 2/7] Adding temperature restore capability for SPEAR. (#36)

* Adding temperature restore capability for SPEAR.

Added parameter SPEAR_ECDA_SST_RESTORE_TFREEZE to allow activation of
sea surface salinity based modification of restoring of temperature.
The formula used is different from the Millero (default in SPEAR runs)
scheme.

* removed spaces on blank line.

* (*)Changed hard wired value to parameter defined in MOM_override

The freezing temperature came from SIS2 code. Changing the default value here to be consistent with that. (-0.054 vs -0.0539)
The salinity restoring code used the -0.0539 value also so answers may change if using that code (RESTORE_SALINITY=T)

* (*)Changed hard wired value to parameter defined in MOM_override

The freezing temperature came from SIS2 code. Changing the default value here to be consistent with that. (-0.054 vs -0.0539)
The salinity restoring code used the -0.0539 value also so answers may change if using that code (RESTORE_SALINITY=T)

* Forgot to replace the salinity masking mulitplier with the override
parameter
---
 .../FMS_cap/MOM_surface_forcing_gfdl.F90       | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

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 09ba9e1156..cab870fed4 100644
--- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
+++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
@@ -113,6 +113,8 @@ module MOM_surface_forcing_gfdl
   real    :: Flux_const                     !< Piston velocity for surface restoring [Z T-1 ~> m s-1]
   real    :: Flux_const_salt                !< Piston velocity for surface salt restoring [Z T-1 ~> m s-1]
   real    :: Flux_const_temp                !< Piston velocity for surface temp restoring [Z T-1 ~> m s-1]
+  logical :: trestore_SPEAR_ECDA            !< If true, modify restoring data wrt local SSS
+  real    :: SPEAR_dTf_dS                   !< The derivative of the freezing temperature with salinity.
   logical :: salt_restore_as_sflux          !< If true, SSS restore as salt flux instead of water flux
   logical :: adjust_net_srestore_to_zero    !< Adjust srestore to zero (for both salt_flux or vprec)
   logical :: adjust_net_srestore_by_scaling !< Adjust srestore w/o moving zero contour
@@ -346,7 +348,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     open_ocn_mask(:,:) = 1.0
     if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice
       do j=js,je ; do i=is,ie
-        if (sfc_state%SST(i,j) <= -0.0539*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
+        if (sfc_state%SST(i,j) <= CS%SPEAR_dTf_dS*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
       enddo ; enddo
     endif
     if (CS%salt_restore_as_sflux) then
@@ -400,6 +402,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   ! SST restoring logic
   if (CS%restore_temp) then
     call time_interp_external(CS%id_trestore, Time, data_restore)
+    if ( CS%trestore_SPEAR_ECDA ) then
+      do j=js,je ; do i=is,ie
+        if (abs(data_restore(i,j)+1.8)<0.0001) then
+          data_restore(i,j) = CS%SPEAR_dTf_dS*sfc_state%SSS(i,j)
+        endif
+      enddo ; enddo
+    endif
+
     do j=js,je ; do i=is,ie
       delta_sst = data_restore(i,j)- sfc_state%SST(i,j)
       delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore)
@@ -1448,7 +1458,13 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
                  "If true, read a file (temp_restore_mask) containing "//&
                  "a mask for SST restoring.", default=.false.)
 
+    call get_param(param_file, mdl, "SPEAR_ECDA_SST_RESTORE_TFREEZE", CS%trestore_SPEAR_ECDA, &
+                 "If true, modify SST restoring field using SSS state. This only modifies the "//&
+                 "restoring data that is within 0.0001degC of -1.8degC.", default=.false.)
   endif
+  call get_param(param_file, mdl, "SPEAR_DTFREEZE_DS", CS%SPEAR_dTf_dS, &
+                 "The derivative of the freezing temperature with salinity.", &
+                 units="deg C PSU-1", default=-0.054, do_not_log=.not.CS%trestore_SPEAR_ECDA)
 
   ! Optionally read tidal amplitude from input file [Z T-1 ~> m s-1] on model grid.
   ! Otherwise use default tidal amplitude for bottom frictionally-generated

From a902e75845b38369bafdc9d0526300e5e8ff3bb0 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Mon, 13 Dec 2021 06:26:03 -0500
Subject: [PATCH 3/7] +Add US args and rescale dt arg to generic tracers

  Added unit_scaling_type arguments to various routines that had previously used
a unit scaling type, but did so via the G%US pointer, to make the type
dependencies more explicit and to avoid unnecessary pointer use.  It had been
the intention to make these arguments explicit from the time they were
introduced via a pointer in the ocean_grid_type as a temporary convenience.
The construct G%US%... was replaced with US%... wherever it was possible.

  Also rescaled some local variables or corrected comments in oil_tracer.F90,
nw2_tracers.F90, and boundary_impulse_tracer.F90, and rescaled the units of
the dt argument to MOM_generic_tracer_column_physics from [s] to [T ~> s].

  All answers are bitwise identical, although there are multiple changes to
public interfaces.
---
 src/core/MOM.F90                              |  8 ++--
 src/core/MOM_boundary_update.F90              |  4 +-
 src/core/MOM_variables.F90                    | 10 ++--
 src/diagnostics/MOM_sum_output.F90            | 11 +++--
 .../MOM_state_initialization.F90              |  2 +-
 .../vertical/MOM_diabatic_driver.F90          | 48 +++++++++----------
 src/tracer/MOM_CFC_cap.F90                    |  5 +-
 src/tracer/MOM_OCMIP2_CFC.F90                 |  5 +-
 src/tracer/MOM_generic_tracer.F90             | 32 +++++++------
 src/tracer/MOM_tracer_Z_init.F90              |  2 +-
 src/tracer/MOM_tracer_flow_control.F90        | 37 +++++++-------
 src/tracer/advection_test_tracer.F90          |  5 +-
 src/tracer/boundary_impulse_tracer.F90        | 27 ++++++-----
 src/tracer/dye_example.F90                    |  5 +-
 src/tracer/ideal_age_example.F90              |  5 +-
 src/tracer/nw2_tracers.F90                    | 32 +++++++------
 src/tracer/oil_tracer.F90                     | 28 ++++++-----
 src/tracer/pseudo_salt_tracer.F90             |  5 +-
 src/tracer/tracer_example.F90                 |  5 +-
 src/user/dyed_channel_initialization.F90      |  5 +-
 src/user/shelfwave_initialization.F90         |  2 +-
 src/user/supercritical_initialization.F90     | 11 ++---
 src/user/tidal_bay_initialization.F90         |  7 +--
 23 files changed, 163 insertions(+), 138 deletions(-)

diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index db114ac3fa..ea54aece44 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -1214,7 +1214,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
     if (associated(CS%tv%T)) call hchksum(CS%tv%T, "Pre-advection T", G%HI, haloshift=1)
     if (associated(CS%tv%S)) call hchksum(CS%tv%S, "Pre-advection S", G%HI, haloshift=1)
     if (associated(CS%tv%frazil)) call hchksum(CS%tv%frazil, "Pre-advection frazil", G%HI, haloshift=0, &
-                                               scale=G%US%Q_to_J_kg*G%US%RZ_to_kg_m2)
+                                               scale=US%Q_to_J_kg*US%RZ_to_kg_m2)
     if (associated(CS%tv%salt_deficit)) call hchksum(CS%tv%salt_deficit, &
                    "Pre-advection salt deficit", G%HI, haloshift=0, scale=US%RZ_to_kg_m2)
   ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G, US)
@@ -1445,7 +1445,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
       if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1)
       if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1)
       if (associated(tv%frazil)) call hchksum(tv%frazil, "Post-diabatic frazil", G%HI, haloshift=0, &
-                                              scale=G%US%Q_to_J_kg*G%US%RZ_to_kg_m2)
+                                              scale=US%Q_to_J_kg*US%RZ_to_kg_m2)
       if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
                                "Post-diabatic salt deficit", G%HI, haloshift=0, scale=US%RZ_to_kg_m2)
     ! call MOM_thermo_chksum("Post-diabatic ", tv, G, US)
@@ -3514,7 +3514,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
                 'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
                 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), &
                 'x=',G%gridLonT(ig), 'y=',G%gridLatT(jg), &
-                'D=',CS%US%Z_to_m*(G%bathyT(i,j)+G%Z_ref),  'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), &
+                'D=',US%Z_to_m*(G%bathyT(i,j)+G%Z_ref),  'SSH=',US%Z_to_m*sfc_state%sea_lev(i,j), &
                 'SST=',sfc_state%SST(i,j), 'SSS=',sfc_state%SSS(i,j), &
                 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), &
                 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J)
@@ -3523,7 +3523,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
                 'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
                 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), &
                 'x=',G%gridLonT(i), 'y=',G%gridLatT(j), &
-                'D=',CS%US%Z_to_m*(G%bathyT(i,j)+G%Z_ref),  'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), &
+                'D=',US%Z_to_m*(G%bathyT(i,j)+G%Z_ref),  'SSH=',US%Z_to_m*sfc_state%sea_lev(i,j), &
                 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), &
                 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J)
             endif
diff --git a/src/core/MOM_boundary_update.F90 b/src/core/MOM_boundary_update.F90
index 286cec20d4..11973f8c02 100644
--- a/src/core/MOM_boundary_update.F90
+++ b/src/core/MOM_boundary_update.F90
@@ -147,13 +147,13 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time)
 ! if (CS%use_files) &
 !     call update_OBC_segment_data(G, GV, OBC, tv, h, Time)
   if (CS%use_tidal_bay) &
-      call tidal_bay_set_OBC_data(OBC, CS%tidal_bay_OBC, G, GV, h, Time)
+      call tidal_bay_set_OBC_data(OBC, CS%tidal_bay_OBC, G, GV, US, h, Time)
   if (CS%use_Kelvin)  &
       call Kelvin_set_OBC_data(OBC, CS%Kelvin_OBC_CSp, G, GV, US, h, Time)
   if (CS%use_shelfwave) &
       call shelfwave_set_OBC_data(OBC, CS%shelfwave_OBC_CSp, G, GV, US, h, Time)
   if (CS%use_dyed_channel) &
-      call dyed_channel_update_flow(OBC, CS%dyed_channel_OBC_CSp, G, GV, Time)
+      call dyed_channel_update_flow(OBC, CS%dyed_channel_OBC_CSp, G, GV, US, Time)
   if (OBC%needs_IO_for_data .or. OBC%add_tide_constituents)  &
       call update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
 
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index f61879845a..5de7ea7319 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -11,6 +11,7 @@ module MOM_variables
 use MOM_EOS,           only : EOS_type
 use MOM_error_handler, only : MOM_error, FATAL
 use MOM_grid,          only : ocean_grid_type
+use MOM_unit_scaling,  only : unit_scale_type
 use MOM_verticalGrid,  only : verticalGrid_type
 
 implicit none ; private
@@ -562,10 +563,11 @@ subroutine dealloc_BT_cont_type(BT_cont)
 end subroutine dealloc_BT_cont_type
 
 !> Diagnostic checksums on various elements of a thermo_var_ptrs type for debugging.
-subroutine MOM_thermovar_chksum(mesg, tv, G)
+subroutine MOM_thermovar_chksum(mesg, tv, G, US)
   character(len=*),      intent(in) :: mesg !< A message that appears in the checksum lines
   type(thermo_var_ptrs), intent(in) :: tv   !< A structure pointing to various thermodynamic variables
   type(ocean_grid_type), intent(in) :: G    !< The ocean's grid structure
+  type(unit_scale_type), intent(in) :: US   !< A dimensional unit scaling type
 
   ! Note that for the chksum calls to be useful for reproducing across PE
   ! counts, there must be no redundant points, so all variables use is..ie
@@ -575,11 +577,11 @@ subroutine MOM_thermovar_chksum(mesg, tv, G)
   if (associated(tv%S)) &
     call hchksum(tv%S, mesg//" tv%S", G%HI)
   if (associated(tv%frazil)) &
-    call hchksum(tv%frazil, mesg//" tv%frazil", G%HI, scale=G%US%Q_to_J_kg*G%US%RZ_to_kg_m2)
+    call hchksum(tv%frazil, mesg//" tv%frazil", G%HI, scale=US%Q_to_J_kg*US%RZ_to_kg_m2)
   if (associated(tv%salt_deficit)) &
-    call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", G%HI, scale=G%US%RZ_to_kg_m2)
+    call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", G%HI, scale=US%RZ_to_kg_m2)
   if (associated(tv%TempxPmE)) &
-    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI, scale=G%US%RZ_to_kg_m2)
+    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI, scale=US%RZ_to_kg_m2)
 end subroutine MOM_thermovar_chksum
 
 end module MOM_variables
diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90
index 3b6fb0c510..668c297658 100644
--- a/src/diagnostics/MOM_sum_output.F90
+++ b/src/diagnostics/MOM_sum_output.F90
@@ -532,7 +532,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
 
   nTr_stocks = 0
   Tr_minmax_avail(:) = .false.
-  call call_tracer_stocks(h, Tr_stocks, G, GV, tracer_CSp, stock_names=Tr_names, &
+  call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, &
                           stock_units=Tr_units, num_stocks=nTr_stocks,&
                           got_min_max=Tr_minmax_avail, global_min=Tr_min, global_max=Tr_max, &
                           xgmin=Tr_min_x, ygmin=Tr_min_y, zgmin=Tr_min_z,&
@@ -1248,7 +1248,7 @@ subroutine write_depth_list(G, US, DL, filename)
   character(len=16) :: depth_chksum, area_chksum
 
   ! All ranks are required to compute the global checksum
-  call get_depth_list_checksums(G, depth_chksum, area_chksum)
+  call get_depth_list_checksums(G, US, depth_chksum, area_chksum)
 
   if (.not.is_root_pe()) return
 
@@ -1313,7 +1313,7 @@ subroutine read_depth_list(G, US, DL, filename, require_chksum, file_matches)
       call MOM_error(WARNING, trim(var_msg) // " some diagnostics may not be reproducible.")
     endif
   else
-    call get_depth_list_checksums(G, depth_grid_chksum, area_grid_chksum)
+    call get_depth_list_checksums(G, US, depth_grid_chksum, area_grid_chksum)
 
     if ((trim(depth_grid_chksum) /= trim(depth_file_chksum)) .or. &
         (trim(area_grid_chksum) /= trim(area_file_chksum)) ) then
@@ -1360,8 +1360,9 @@ end subroutine read_depth_list
 !!
 !! Checksums are saved as hexadecimal strings, in order to avoid potential
 !! datatype issues with netCDF attributes.
-subroutine get_depth_list_checksums(G, depth_chksum, area_chksum)
+subroutine get_depth_list_checksums(G, US, depth_chksum, area_chksum)
   type(ocean_grid_type), intent(in) :: G          !< Ocean grid structure
+  type(unit_scale_type), intent(in) :: US         !< A dimensional unit scaling type
   character(len=16), intent(out) :: depth_chksum  !< Depth checksum hexstring
   character(len=16), intent(out) :: area_chksum   !< Area checksum hexstring
 
@@ -1378,7 +1379,7 @@ subroutine get_depth_list_checksums(G, depth_chksum, area_chksum)
 
   ! Area checksum
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    field(i,j) = G%mask2dT(i,j) * G%US%L_to_m**2*G%areaT(i,j)
+    field(i,j) = G%mask2dT(i,j) * US%L_to_m**2*G%areaT(i,j)
   enddo ; enddo
   write(area_chksum, '(Z16)') field_chksum(field(:,:))
 
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 2aab378b4a..8055440cce 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -620,7 +620,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
     elseif (trim(config) == "shelfwave") then
       OBC%update_OBC = .true.
     elseif (lowercase(trim(config)) == "supercritical") then
-      call supercritical_set_OBC_data(OBC, G, GV, PF)
+      call supercritical_set_OBC_data(OBC, G, GV, US, PF)
     elseif (trim(config) == "tidal_bay") then
       OBC%update_OBC = .true.
     elseif (trim(config) == "USER") then
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 77ec87b230..c123e60800 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -601,7 +601,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
   if (CS%debug) then
     call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after set_diffusivity ", tv, G)
+    call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
     call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
   endif
 
@@ -678,7 +678,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
     if (CS%debug) then
       call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
-      call MOM_thermovar_chksum("after KPP", tv, G)
+      call MOM_thermovar_chksum("after KPP", tv, G, US)
       call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
       call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
       call hchksum(CS%KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, scale=GV%H_to_m*US%s_to_T)
@@ -699,7 +699,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
     if (CS%debug) then
       call MOM_state_chksum("after KPP_applyNLT ", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0)
-      call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G)
+      call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US)
     endif
   endif ! endif for KPP
 
@@ -747,7 +747,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
 
   if (CS%debug) then
     call MOM_forcing_chksum("after calc_entrain ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after calc_entrain ", tv, G)
+    call MOM_thermovar_chksum("after calc_entrain ", tv, G, US)
     call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, US, haloshift=0)
     call hchksum(ent_s, "after calc_entrain ent_s", G%HI, haloshift=0, scale=GV%H_to_m)
   endif
@@ -842,7 +842,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
   call cpu_clock_end(id_clock_remap)
   if (CS%debug) then
     call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G)
+    call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G, US)
     call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, US, haloshift=0)
   endif
   if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)")
@@ -851,7 +851,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
   if (CS%debug) then
     call MOM_state_chksum("after negative check ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after negative check ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after negative check ", tv, G)
+    call MOM_thermovar_chksum("after negative check ", tv, G, US)
   endif
   if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)")
   if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G, GV, US)
@@ -908,7 +908,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
 
   if (CS%debug) then
     call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, US, haloshift=0)
-    call MOM_thermovar_chksum("after mixed layer ", tv, G)
+    call MOM_thermovar_chksum("after mixed layer ", tv, G, US)
   endif
 
   ! Whenever thickness changes let the diag manager know, as the
@@ -1020,7 +1020,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
     call cpu_clock_end(id_clock_sponge)
     if (CS%debug) then
       call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, US, haloshift=0)
-      call MOM_thermovar_chksum("apply_sponge ", tv, G)
+      call MOM_thermovar_chksum("apply_sponge ", tv, G, US)
     endif
   endif ! CS%use_sponge
 
@@ -1032,7 +1032,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
     call cpu_clock_end(id_clock_oda_incupd)
     if (CS%debug) then
       call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0)
-      call MOM_thermovar_chksum("apply_oda_incupd ", tv, G)
+      call MOM_thermovar_chksum("apply_oda_incupd ", tv, G, US)
     endif
   endif ! CS%use_oda_incupd
 
@@ -1185,7 +1185,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
   if (CS%debug) then
     call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after set_diffusivity ", tv, G)
+    call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
     call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
   endif
 
@@ -1253,7 +1253,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
     if (CS%debug) then
       call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
-      call MOM_thermovar_chksum("after KPP", tv, G)
+      call MOM_thermovar_chksum("after KPP", tv, G, US)
       call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
       call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
       call hchksum(CS%KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, scale=GV%H_to_m*US%s_to_T)
@@ -1274,7 +1274,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
     if (CS%debug) then
       call MOM_state_chksum("after KPP_applyNLT ", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0)
-      call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G)
+      call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US)
     endif
   endif ! endif for KPP
 
@@ -1372,7 +1372,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
   call cpu_clock_end(id_clock_remap)
   if (CS%debug) then
     call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G)
+    call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G, US)
     call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, US, haloshift=0)
   endif
   if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)")
@@ -1439,7 +1439,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
 
   if (CS%debug) then
     call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, US, haloshift=0)
-    call MOM_thermovar_chksum("after mixed layer ", tv, G)
+    call MOM_thermovar_chksum("after mixed layer ", tv, G, US)
   endif
 
   ! Whenever thickness changes let the diag manager know, as the
@@ -1526,7 +1526,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
     call cpu_clock_end(id_clock_sponge)
     if (CS%debug) then
       call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, US, haloshift=0)
-      call MOM_thermovar_chksum("apply_sponge ", tv, G)
+      call MOM_thermovar_chksum("apply_sponge ", tv, G, US)
     endif
   endif ! CS%use_sponge
 
@@ -1538,7 +1538,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
     call cpu_clock_end(id_clock_oda_incupd)
     if (CS%debug) then
       call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0)
-      call MOM_thermovar_chksum("apply_oda_incupd ", tv, G)
+      call MOM_thermovar_chksum("apply_oda_incupd ", tv, G, US)
     endif
   endif ! CS%use_oda_incupd
 
@@ -1789,7 +1789,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
   if (CS%debug) then
     call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after set_diffusivity ", tv, G)
+    call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
     call hchksum(Kd_lay, "after set_diffusivity Kd_lay", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
     call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
   endif
@@ -1863,7 +1863,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
     if (CS%debug) then
       call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
-      call MOM_thermovar_chksum("after KPP", tv, G)
+      call MOM_thermovar_chksum("after KPP", tv, G, US)
       call hchksum(Kd_lay, "after KPP Kd_lay", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
       call hchksum(Kd_Int, "after KPP Kd_Int", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
     endif
@@ -1896,7 +1896,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
     if (CS%debug) then
       call MOM_state_chksum("after KPP_applyNLT ", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0)
-      call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G)
+      call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US)
     endif
   endif ! endif for KPP
 
@@ -1934,7 +1934,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
 
   if (CS%debug) then
     call MOM_forcing_chksum("after calc_entrain ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after calc_entrain ", tv, G)
+    call MOM_thermovar_chksum("after calc_entrain ", tv, G, US)
     call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, US, haloshift=0)
     call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, scale=GV%H_to_m)
     call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, scale=GV%H_to_m)
@@ -1985,7 +1985,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
   if (CS%debug) then
     call MOM_state_chksum("after negative check ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after negative check ", fluxes, G, US, haloshift=0)
-    call MOM_thermovar_chksum("after negative check ", tv, G)
+    call MOM_thermovar_chksum("after negative check ", tv, G, US)
   endif
   if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)")
   if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G, GV, US)
@@ -2183,7 +2183,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
 
   if (CS%debug) then
     call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, US, haloshift=0)
-    call MOM_thermovar_chksum("after mixed layer ", tv, G)
+    call MOM_thermovar_chksum("after mixed layer ", tv, G, US)
     call hchksum(ea, "after mixed layer ea", G%HI, scale=GV%H_to_m)
     call hchksum(eb, "after mixed layer eb", G%HI, scale=GV%H_to_m)
   endif
@@ -2331,7 +2331,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
     call cpu_clock_end(id_clock_sponge)
     if (CS%debug) then
       call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, US, haloshift=0)
-      call MOM_thermovar_chksum("apply_sponge ", tv, G)
+      call MOM_thermovar_chksum("apply_sponge ", tv, G, US)
     endif
   endif ! CS%use_sponge
 
@@ -2342,7 +2342,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
     call cpu_clock_end(id_clock_oda_incupd)
     if (CS%debug) then
       call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0)
-      call MOM_thermovar_chksum("apply_oda_incupd ", tv, G)
+      call MOM_thermovar_chksum("apply_oda_incupd ", tv, G, US)
     endif
   endif ! CS%use_oda_incupd
 
diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90
index 79df57cc23..7296f1d469 100644
--- a/src/tracer/MOM_CFC_cap.F90
+++ b/src/tracer/MOM_CFC_cap.F90
@@ -341,13 +341,14 @@ end subroutine CFC_cap_column_physics
 !> Calculates the mass-weighted integral of all tracer stocks,
 !! returning the number of stocks it has calculated.  If the stock_index
 !! is present, only the stock corresponding to that coded index is returned.
-function CFC_cap_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function CFC_cap_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),           intent(in)    :: G      !< The ocean's grid structure.
   type(verticalGrid_type),         intent(in)    :: GV     !< The ocean's vertical grid structure.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                    intent(in)    :: h      !< Layer thicknesses [H ~> m or kg m-2].
   real, dimension(:),              intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                            !! tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),           intent(in)    :: US     !< A dimensional unit scaling type
   type(CFC_cap_CS),                pointer       :: CS     !< The control structure returned by a
                                                            !! previous call to register_CFC_cap.
   character(len=*), dimension(:),  intent(out)   :: names  !< The names of the stocks calculated.
@@ -376,7 +377,7 @@ function CFC_cap_stock(h, stocks, G, GV, CS, names, units, stock_index)
   call query_vardesc(CS%CFC12_desc, name=names(2), units=units(2), caller="CFC_cap_stock")
   units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg"
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   stocks(1) = 0.0 ; stocks(2) = 0.0
   do k=1,nz ; do j=js,je ; do i=is,ie
     mass = G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k)
diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90
index 38193a3abc..5fe55b896b 100644
--- a/src/tracer/MOM_OCMIP2_CFC.F90
+++ b/src/tracer/MOM_OCMIP2_CFC.F90
@@ -478,13 +478,14 @@ end subroutine OCMIP2_CFC_column_physics
 !> This function calculates the mass-weighted integral of all tracer stocks,
 !! returning the number of stocks it has calculated.  If the stock_index
 !! is present, only the stock corresponding to that coded index is returned.
-function OCMIP2_CFC_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function OCMIP2_CFC_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),           intent(in)    :: G      !< The ocean's grid structure.
   type(verticalGrid_type),         intent(in)    :: GV     !< The ocean's vertical grid structure.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                    intent(in)    :: h      !< Layer thicknesses [H ~> m or kg m-2].
   real, dimension(:),              intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                            !! tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),           intent(in)    :: US     !< A dimensional unit scaling type
   type(OCMIP2_CFC_CS),             pointer       :: CS     !< The control structure returned by a
                                                            !! previous call to register_OCMIP2_CFC.
   character(len=*), dimension(:),  intent(out)   :: names  !< The names of the stocks calculated.
@@ -513,7 +514,7 @@ function OCMIP2_CFC_stock(h, stocks, G, GV, CS, names, units, stock_index)
   call query_vardesc(CS%CFC12_desc, name=names(2), units=units(2), caller="OCMIP2_CFC_stock")
   units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg"
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   stocks(1) = 0.0 ; stocks(2) = 0.0
   do k=1,nz ; do j=js,je ; do i=is,ie
     mass = G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k)
diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90
index 4627d0ec80..bf9f01e266 100644
--- a/src/tracer/MOM_generic_tracer.F90
+++ b/src/tracer/MOM_generic_tracer.F90
@@ -395,7 +395,7 @@ end subroutine initialize_MOM_generic_tracer
   !! tracer physics or chemistry to the tracers from this file.
   !! CFCs are relatively simple, as they are passive tracers. with only a surface
   !! flux as a source.
-  subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, CS, tv, optics, &
+  subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, US, CS, tv, optics, &
         evap_CFL_limit, minimum_forcing_depth)
     type(ocean_grid_type),   intent(in) :: G     !< The ocean's grid structure
     type(verticalGrid_type), intent(in) :: GV    !< The ocean's vertical grid structure
@@ -412,7 +412,8 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
     type(forcing),           intent(in) :: fluxes !< A structure containing pointers to thermodynamic
                                                  !! and tracer forcing fields.
     real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Hml  !< Mixed layer depth [Z ~> m]
-    real,                    intent(in) :: dt    !< The amount of time covered by this call [s]
+    real,                    intent(in) :: dt    !< The amount of time covered by this call [T ~> s]
+    type(unit_scale_type),   intent(in) :: US    !< A dimensional unit scaling type
     type(MOM_generic_tracer_CS), pointer :: CS   !< Pointer to the control structure for this module.
     type(thermo_var_ptrs),   intent(in) :: tv    !< A structure pointing to various thermodynamic variables
     type(optics_type),       intent(in) :: optics !< The structure containing optical properties.
@@ -469,7 +470,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
           call g_tracer_get_pointer(g_tracer,g_tracer_name,'runoff_tracer_flux',runoff_tracer_flux_array)
           !nnz: Why is fluxes%river = 0?
           runoff_tracer_flux_array(:,:) = trunoff_array(:,:) * &
-                   G%US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:)
+                   US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:)
           stf_array = stf_array + runoff_tracer_flux_array
        endif
 
@@ -496,25 +497,25 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
     dz_ml(:,:) = 0.0
     do j=jsc,jec ; do i=isc,iec
       surface_field(i,j) = tv%S(i,j,1)
-      dz_ml(i,j) = G%US%Z_to_m * Hml(i,j)
+      dz_ml(i,j) = US%Z_to_m * Hml(i,j)
     enddo ; enddo
     sosga = global_area_mean(surface_field, G)
 
     !
     !Calculate tendencies (i.e., field changes at dt) from the sources / sinks
     !
-    if ((G%US%L_to_m == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. (G%US%s_to_T == 1.0)) then
+    if ((US%L_to_m == 1.0) .and. (US%RZ_to_kg_m2 == 1.0) .and. (US%s_to_T == 1.0)) then
       ! Avoid unnecessary copies when no unit conversion is needed.
       call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
                G%areaT, get_diag_time_end(CS%diag), &
                optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, &
                internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga)
     else
-      call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
-               G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), &
+      call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, US%T_to_s*dt, &
+               US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), &
                optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, &
-               internal_heat=G%US%RZ_to_kg_m2*tv%internal_heat(:,:), &
-               frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga)
+               internal_heat=US%RZ_to_kg_m2*tv%internal_heat(:,:), &
+               frunoff=US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga)
     endif
 
     ! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes
@@ -526,7 +527,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
           do k=1,nk ;do j=jsc,jec ; do i=isc,iec
             h_work(i,j,k) = h_old(i,j,k)
           enddo ; enddo ; enddo
-          call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), G%US%s_to_T*dt, &
+          call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, &
                             fluxes, h_work, evap_CFL_limit, minimum_forcing_depth)
         endif
 
@@ -544,16 +545,16 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
     ! surface source is applied and diapycnal advection and diffusion occurs.
     if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
       ! Last arg is tau which is always 1 for MOM6
-      call generic_tracer_vertdiff_G(h_work, ea, eb, dt, GV%kg_m2_to_H, GV%m_to_H, 1)
+      call generic_tracer_vertdiff_G(h_work, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1)
     else
       ! Last arg is tau which is always 1 for MOM6
-      call generic_tracer_vertdiff_G(h_old, ea, eb, dt, GV%kg_m2_to_H, GV%m_to_H, 1)
+      call generic_tracer_vertdiff_G(h_old, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1)
     endif
 
     ! Update bottom fields after vertical processes
 
     ! Second arg is tau which is always 1 for MOM6
-    call generic_tracer_update_from_bottom(dt, 1, get_diag_time_end(CS%diag))
+    call generic_tracer_update_from_bottom(US%T_to_s*dt, 1, get_diag_time_end(CS%diag))
 
     !Output diagnostics via diag_manager for all generic tracers and their fluxes
     call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1)
@@ -568,12 +569,13 @@ end subroutine MOM_generic_tracer_column_physics
   !! being requested specifically, returning the number of stocks it has
   !! calculated. If the stock_index is present, only the stock corresponding
   !! to that coded index is returned.
-  function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index)
+  function MOM_generic_tracer_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
     type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
     type(verticalGrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure
     real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
     real, dimension(:),                 intent(out)   :: stocks !< The mass-weighted integrated amount of each
                                                                 !! tracer, in kg times concentration units [kg conc].
+    type(unit_scale_type),              intent(in)    :: US     !< A dimensional unit scaling type
     type(MOM_generic_tracer_CS),        pointer       :: CS     !< Pointer to the control structure for this module.
     character(len=*), dimension(:),     intent(out)   :: names  !< The names of the stocks calculated.
     character(len=*), dimension(:),     intent(out)   :: units  !< The units of the stocks calculated.
@@ -604,7 +606,7 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_inde
 
     if (.NOT. associated(CS%g_tracer_list)) return ! No stocks.
 
-    stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+    stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
     m=1 ; g_tracer=>CS%g_tracer_list
     do
       call g_tracer_get_alias(g_tracer,names(m))
diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90
index 1be976d3f2..e8324b6043 100644
--- a/src/tracer/MOM_tracer_Z_init.F90
+++ b/src/tracer/MOM_tracer_Z_init.F90
@@ -565,7 +565,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h,
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                  intent(inout) :: temp !< potential temperature [degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
-                                 intent(inout) :: salt !< salinity [PSU]
+                                 intent(inout) :: salt !< salinity [ppt]
   real, dimension(SZK_(GV)),     intent(in)    :: R_tgt !< desired potential density [R ~> kg m-3].
   real,                          intent(in)    :: p_ref !< reference pressure [R L2 T-2 ~> Pa].
   integer,                       intent(in)    :: niter !< maximum number of iterations
diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90
index 4278594913..2ae72a3270 100644
--- a/src/tracer/MOM_tracer_flow_control.F90
+++ b/src/tracer/MOM_tracer_flow_control.F90
@@ -262,13 +262,13 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
     register_pseudo_salt_tracer(HI, GV, param_file,  CS%pseudo_salt_tracer_CSp, &
                                 tr_Reg, restart_CS)
   if (CS%use_boundary_impulse_tracer) CS%use_boundary_impulse_tracer = &
-    register_boundary_impulse_tracer(HI, GV, param_file,  CS%boundary_impulse_tracer_CSp, &
+    register_boundary_impulse_tracer(HI, GV, US, param_file,  CS%boundary_impulse_tracer_CSp, &
                                      tr_Reg, restart_CS)
   if (CS%use_dyed_obc_tracer) CS%use_dyed_obc_tracer = &
     register_dyed_obc_tracer(HI, GV, param_file, CS%dyed_obc_tracer_CSp, &
                              tr_Reg, restart_CS)
   if (CS%use_nw2_tracers) CS%use_nw2_tracers = &
-    register_nw2_tracers(HI, GV, param_file,  CS%nw2_tracers_CSp, tr_Reg, restart_CS)
+    register_nw2_tracers(HI, GV, US, param_file,  CS%nw2_tracers_CSp, tr_Reg, restart_CS)
 
 end subroutine call_tracer_register
 
@@ -346,7 +346,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag
     call initialize_pseudo_salt_tracer(restart, day, G, GV, h, diag, OBC, CS%pseudo_salt_tracer_CSp, &
                                 sponge_CSp, tv)
   if (CS%use_boundary_impulse_tracer) &
-    call initialize_boundary_impulse_tracer(restart, day, G, GV, h, diag, OBC, CS%boundary_impulse_tracer_CSp, &
+    call initialize_boundary_impulse_tracer(restart, day, G, GV, US, h, diag, OBC, CS%boundary_impulse_tracer_CSp, &
                                 sponge_CSp, tv)
   if (CS%use_dyed_obc_tracer) &
     call initialize_dyed_obc_tracer(restart, day, G, GV, h, diag, OBC, CS%dyed_obc_tracer_CSp)
@@ -495,8 +495,8 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV,
       if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//&
             "has not been written to permit dimensionsal rescaling.  Set all 4 of the "//&
             "[QRZT]_RESCALE_POWER parameters to 0.")
-      call MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, US%T_to_s*dt, &
-                                             G, GV, CS%MOM_generic_tracer_CSp, tv, optics, &
+      call MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, &
+                                             G, GV, US, CS%MOM_generic_tracer_CSp, tv, optics, &
                                              evap_CFL_limit=evap_CFL_limit, &
                                              minimum_forcing_depth=minimum_forcing_depth)
     endif
@@ -555,8 +555,8 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV,
       if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//&
             "has not been written to permit dimensionsal rescaling.  Set all 4 of the "//&
             "[QRZT]_RESCALE_POWER parameters to 0.")
-      call MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, US%T_to_s*dt, &
-                                     G, GV, CS%MOM_generic_tracer_CSp, tv, optics)
+      call MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, &
+                                     G, GV, US, CS%MOM_generic_tracer_CSp, tv, optics)
     endif
     if (CS%use_pseudo_salt_tracer) &
       call pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
@@ -575,7 +575,7 @@ end subroutine call_tracer_column_fns
 
 !> This subroutine calls all registered tracer packages to enable them to
 !! add to the surface state returned to the coupler. These routines are optional.
-subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_units, &
+subroutine call_tracer_stocks(h, stock_values, G, GV, US, CS, stock_names, stock_units, &
                               num_stocks, stock_index, got_min_max, global_min, global_max, &
                               xgmin, ygmin, zgmin, xgmax, ygmax, zgmax)
   type(ocean_grid_type),          intent(in)  :: G           !< The ocean's grid structure.
@@ -584,6 +584,7 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni
                                   intent(in)  :: h           !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),             intent(out) :: stock_values !< The integrated amounts of a tracer
                              !! on the current PE, usually in kg x concentration [kg conc].
+  type(unit_scale_type),          intent(in)  :: US          !< A dimensional unit scaling type
   type(tracer_flow_control_CS),   pointer     :: CS          !< The control structure returned by a
                                                              !! previous call to
                                                              !! call_tracer_register.
@@ -624,7 +625,7 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni
 
 !  Add other user-provided calls here.
   if (CS%use_USER_tracer_example) then
-    ns = USER_tracer_stock(h, values, G, GV, CS%USER_tracer_example_CSp, &
+    ns = USER_tracer_stock(h, values, G, GV, US, CS%USER_tracer_example_CSp, &
                            names, units, stock_index)
     call store_stocks("tracer_example", ns, names, units, values, index, stock_values, &
                        set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
@@ -636,44 +637,44 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni
 !                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
 ! endif
   if (CS%use_ideal_age) then
-    ns = ideal_age_stock(h, values, G, GV, CS%ideal_age_tracer_CSp, &
+    ns = ideal_age_stock(h, values, G, GV, US, CS%ideal_age_tracer_CSp, &
                          names, units, stock_index)
     call store_stocks("ideal_age_example", ns, names, units, values, index, &
            stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
   if (CS%use_regional_dyes) then
-    ns = dye_stock(h, values, G, GV, CS%dye_tracer_CSp, &
+    ns = dye_stock(h, values, G, GV, US, CS%dye_tracer_CSp, &
                          names, units, stock_index)
     call store_stocks("regional_dyes", ns, names, units, values, index, &
            stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
   if (CS%use_oil) then
-    ns = oil_stock(h, values, G, GV, CS%oil_tracer_CSp, &
+    ns = oil_stock(h, values, G, GV, US, CS%oil_tracer_CSp, &
                          names, units, stock_index)
     call store_stocks("oil_tracer", ns, names, units, values, index, &
            stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
   if (CS%use_OCMIP2_CFC) then
-    ns = OCMIP2_CFC_stock(h, values, G, GV, CS%OCMIP2_CFC_CSp, names, units, stock_index)
+    ns = OCMIP2_CFC_stock(h, values, G, GV, US, CS%OCMIP2_CFC_CSp, names, units, stock_index)
     call store_stocks("MOM_OCMIP2_CFC", ns, names, units, values, index, stock_values, &
                        set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_CFC_cap) then
-    ns = CFC_cap_stock(h, values, G, GV, CS%CFC_cap_CSp, names, units, stock_index)
+    ns = CFC_cap_stock(h, values, G, GV, US, CS%CFC_cap_CSp, names, units, stock_index)
     call store_stocks("MOM_CFC_cap", ns, names, units, values, index, stock_values, &
                        set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_advection_test_tracer) then
-    ns = advection_test_stock( h, values, G, GV, CS%advection_test_tracer_CSp, &
+    ns = advection_test_stock( h, values, G, GV, US, CS%advection_test_tracer_CSp, &
                          names, units, stock_index )
     call store_stocks("advection_test_tracer", ns, names, units, values, index, &
            stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_MOM_generic_tracer) then
-    ns = MOM_generic_tracer_stock(h, values, G, GV, CS%MOM_generic_tracer_CSp, &
+    ns = MOM_generic_tracer_stock(h, values, G, GV, US, CS%MOM_generic_tracer_CSp, &
                                    names, units, stock_index)
     call store_stocks("MOM_generic_tracer", ns, names, units, values, index, stock_values, &
                        set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
@@ -684,14 +685,14 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni
 
   endif
   if (CS%use_pseudo_salt_tracer) then
-    ns = pseudo_salt_stock(h, values, G, GV, CS%pseudo_salt_tracer_CSp, &
+    ns = pseudo_salt_stock(h, values, G, GV, US, CS%pseudo_salt_tracer_CSp, &
                          names, units, stock_index)
     call store_stocks("pseudo_salt_tracer", ns, names, units, values, index, &
            stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_boundary_impulse_tracer) then
-    ns = boundary_impulse_stock(h, values, G, GV, CS%boundary_impulse_tracer_CSp, &
+    ns = boundary_impulse_stock(h, values, G, GV, US, CS%boundary_impulse_tracer_CSp, &
                          names, units, stock_index)
     call store_stocks("boundary_impulse_tracer", ns, names, units, values, index, &
            stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90
index b713803182..8fdb525b4a 100644
--- a/src/tracer/advection_test_tracer.F90
+++ b/src/tracer/advection_test_tracer.F90
@@ -344,12 +344,13 @@ end subroutine advection_test_tracer_surface_state
 
 !> Calculate the mass-weighted integral of all tracer stocks, returning the number of stocks it has calculated.
 !!  If the stock_index is present, only the stock corresponding to that coded index is returned.
-function advection_test_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function advection_test_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),              intent(in)    :: G      !< The ocean's grid structure
   type(verticalGrid_type),            intent(in)    :: GV     !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h   !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),                 intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                               !! tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),              intent(in)    :: US     !< A dimensional unit scaling type
   type(advection_test_tracer_CS),     pointer       :: CS     !< The control structure returned by a previous
                                                               !! call to register_advection_test_tracer.
   character(len=*), dimension(:),     intent(out)   :: names  !< the names of the stocks calculated.
@@ -373,7 +374,7 @@ function advection_test_stock(h, stocks, G, GV, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="advection_test_stock")
     stocks(m) = 0.0
diff --git a/src/tracer/boundary_impulse_tracer.F90 b/src/tracer/boundary_impulse_tracer.F90
index 18e9b8dc8e..ea60a09608 100644
--- a/src/tracer/boundary_impulse_tracer.F90
+++ b/src/tracer/boundary_impulse_tracer.F90
@@ -46,9 +46,8 @@ module boundary_impulse_tracer
 
   integer :: nkml !< Number of layers in mixed layer
   real, dimension(NTR_MAX)  :: land_val = -1.0 !< A value to use to fill in tracers over land
-  real :: kw_eff !< An effective piston velocity used to flux tracer out at the surface
   real :: remaining_source_time !< How much longer (same units as the timestep) to
-                                !! inject the tracer at the surface [s]
+                                !! inject the tracer at the surface [T ~> s]
 
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                    !! regulate the timing of diagnostic output.
@@ -60,9 +59,10 @@ module boundary_impulse_tracer
 contains
 
 !> Read in runtime options and add boundary impulse tracer to tracer registry
-function register_boundary_impulse_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
+function register_boundary_impulse_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
   type(hor_index_type),             intent(in   ) :: HI   !< A horizontal index type structure
   type(verticalGrid_type),          intent(in   ) :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),            intent(in   ) :: US   !< A dimensional unit scaling type
   type(param_file_type),            intent(in   ) :: param_file !< A structure to parse for run-time parameters
   type(boundary_impulse_tracer_CS), pointer       :: CS   !< The control structure returned by a previous
                                                           !! call to register_boundary_impulse_tracer.
@@ -79,7 +79,7 @@ function register_boundary_impulse_tracer(HI, GV, param_file, CS, tr_Reg, restar
   character(len=48)  :: flux_units ! The units for tracer fluxes, usually
                             ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
   ! This include declares and sets the variable "version".
-#include "version_variable.h"
+# include "version_variable.h"
   real, pointer :: tr_ptr(:,:,:) => NULL()
   real, pointer :: rem_time_ptr => NULL()
   logical :: register_boundary_impulse_tracer
@@ -99,7 +99,7 @@ function register_boundary_impulse_tracer(HI, GV, param_file, CS, tr_Reg, restar
                  "Length of time for the boundary tracer to be injected "//&
                  "into the mixed layer. After this time has elapsed, the "//&
                  "surface becomes a sink for the boundary impulse tracer.", &
-                 default=31536000.0)
+                 default=31536000.0, scale=US%s_to_T)
   call get_param(param_file, mdl, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, &
                  "If true, tracers may go through the initialization code "//&
                  "if they are not found in the restart files.  Otherwise "//&
@@ -145,13 +145,14 @@ function register_boundary_impulse_tracer(HI, GV, param_file, CS, tr_Reg, restar
 end function register_boundary_impulse_tracer
 
 !> Initialize tracer from restart or set to 1 at surface to initialize
-subroutine initialize_boundary_impulse_tracer(restart, day, G, GV, h, diag, OBC, CS, &
+subroutine initialize_boundary_impulse_tracer(restart, day, G, GV, US, h, diag, OBC, CS, &
                                   sponge_CSp, tv)
   logical,                            intent(in) :: restart !< .true. if the fields have already
                                                          !! been read from a restart file.
   type(time_type),            target, intent(in) :: day  !< Time of the start of the run.
   type(ocean_grid_type),              intent(in) :: G    !< The ocean's grid structure
   type(verticalGrid_type),            intent(in) :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),              intent(in) :: US   !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                       intent(in) :: h    !< Layer thicknesses [H ~> m or kg m-2]
   type(diag_ctrl),            target, intent(in) :: diag !< A structure that is used to regulate
@@ -186,14 +187,17 @@ subroutine initialize_boundary_impulse_tracer(restart, day, G, GV, h, diag, OBC,
 
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=name, caller="initialize_boundary_impulse_tracer")
-    if ((.not.restart) .or. (.not. &
-        query_initialized(CS%tr(:,:,:,m), name, CS%restart_CSp))) then
+    if ((.not.restart) .or. (.not. query_initialized(CS%tr(:,:,:,m), name, CS%restart_CSp))) then
       do k=1,CS%nkml ; do j=jsd,jed ; do i=isd,ied
         CS%tr(i,j,k,m) = 1.0
       enddo ; enddo ; enddo
     endif
   enddo ! Tracer loop
 
+  if (restart .and. (US%s_to_T_restart /= 0.0) .and. (US%s_to_T /= US%s_to_T_restart) ) then
+    CS%remaining_source_time = (US%s_to_T / US%s_to_T_restart) * CS%remaining_source_time
+  endif
+
   if (associated(OBC)) then
   ! Steal from updated DOME in the fullness of time.
   endif
@@ -268,7 +272,7 @@ subroutine boundary_impulse_tracer_column_physics(h_old, h_new, ea, eb, fluxes,
       do k=1,CS%nkml ; do j=js,je ; do i=is,ie
         CS%tr(i,j,k,m) = 1.0
       enddo ; enddo ; enddo
-      CS%remaining_source_time = CS%remaining_source_time-US%T_to_s*dt
+      CS%remaining_source_time = CS%remaining_source_time-dt
     else
       do k=1,CS%nkml ; do j=js,je ; do i=is,ie
         CS%tr(i,j,k,m) = 0.0
@@ -283,12 +287,13 @@ end subroutine boundary_impulse_tracer_column_physics
 !> This function calculates the mass-weighted integral of the boundary impulse,
 !! tracer stocks returning the number of stocks it has calculated.  If the stock_index
 !! is present, only the stock corresponding to that coded index is returned.
-function boundary_impulse_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function boundary_impulse_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),                    intent(in   ) :: G    !< The ocean's grid structure
   type(verticalGrid_type),                  intent(in   ) :: GV   !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in   ) :: h    !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),                       intent(  out) :: stocks !< the mass-weighted integrated amount of each
                                                                   !! tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),                    intent(in   ) :: US   !< A dimensional unit scaling type
   type(boundary_impulse_tracer_CS),         pointer       :: CS   !< The control structure returned by a previous
                                                                   !! call to register_boundary_impulse_tracer.
   character(len=*), dimension(:),           intent(  out) :: names  !< The names of the stocks calculated.
@@ -317,7 +322,7 @@ function boundary_impulse_stock(h, stocks, G, GV, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,1
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="boundary_impulse_stock")
     units(m) = trim(units(m))//" kg"
diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90
index 91806bb94e..dca01e974a 100644
--- a/src/tracer/dye_example.F90
+++ b/src/tracer/dye_example.F90
@@ -325,12 +325,13 @@ end subroutine dye_tracer_column_physics
 !> This function calculates the mass-weighted integral of all tracer stocks,
 !! returning the number of stocks it has calculated.  If the stock_index
 !! is present, only the stock corresponding to that coded index is returned.
-function dye_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function dye_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h  !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),                 intent(out)   :: stocks !< the mass-weighted integrated amount of
                                                             !! each tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),              intent(in)    :: US   !< A dimensional unit scaling type
   type(dye_tracer_CS),                pointer       :: CS   !< The control structure returned by a
                                                             !! previous call to register_dye_tracer.
   character(len=*), dimension(:),     intent(out)   :: names !< the names of the stocks calculated.
@@ -356,7 +357,7 @@ function dye_stock(h, stocks, G, GV, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="dye_stock")
     units(m) = trim(units(m))//" kg"
diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90
index ca47a8ca1d..d5c813b3d0 100644
--- a/src/tracer/ideal_age_example.F90
+++ b/src/tracer/ideal_age_example.F90
@@ -369,13 +369,14 @@ end subroutine ideal_age_tracer_column_physics
 
 !> Calculates the mass-weighted integral of all tracer stocks, returning the number of stocks it
 !! has calculated.  If stock_index is present, only the stock corresponding to that coded index is found.
-function ideal_age_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function ideal_age_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                       intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),                 intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                             !! tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),              intent(in)    :: US   !< A dimensional unit scaling type
   type(ideal_age_tracer_CS),          pointer       :: CS   !< The control structure returned by a previous
                                                             !! call to register_ideal_age_tracer.
   character(len=*), dimension(:),     intent(out)   :: names  !< the names of the stocks calculated.
@@ -400,7 +401,7 @@ function ideal_age_stock(h, stocks, G, GV, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="ideal_age_stock")
     units(m) = trim(units(m))//" kg"
diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
index fcb9f3e854..0e66ebbcf3 100644
--- a/src/tracer/nw2_tracers.F90
+++ b/src/tracer/nw2_tracers.F90
@@ -33,7 +33,8 @@ module nw2_tracers
   type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
   type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry
   real, pointer :: tr(:,:,:,:) => NULL()   !< The array of tracers used in this package, in g m-3?
-  real, allocatable , dimension(:) :: restore_rate !< The exponential growth rate for restoration value [year-1].
+  real, allocatable , dimension(:) :: restore_rate !< The rate at which the tracer is damped toward
+                                             !! its target profile [T-1 ~> s-1]
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                              !! regulate the timing of diagnostic output.
   type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart controls structure
@@ -42,9 +43,10 @@ module nw2_tracers
 contains
 
 !> Register the NW2 tracer fields to be used with MOM.
-logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS)
+logical function register_nw2_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
   type(hor_index_type),       intent(in) :: HI   !< A horizontal index type structure
   type(verticalGrid_type),    intent(in) :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),      intent(in) :: US   !< A dimensional unit scaling type
   type(param_file_type),      intent(in) :: param_file !< A structure to parse for run-time parameters
   type(nw2_tracers_CS),       pointer    :: CS !< The control structure returned by a previous
                                                !! call to register_nw2_tracer.
@@ -62,7 +64,7 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
   logical :: do_nw2
   integer :: isd, ied, jsd, jed, nz, m, ig
   integer :: n_groups ! Number of groups of three tracers (i.e. # tracers/3)
-  real, allocatable, dimension(:) :: timescale_in_days
+  real, allocatable, dimension(:) :: timescale_in_days ! Damping timescale [days]
   type(vardesc) :: tr_desc ! Descriptions and metadata for the tracers
   isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
 
@@ -100,7 +102,7 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
     call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, &
                          registry_diags=.true., restart_CS=restart_CS, mandatory=.false.)
     ig = int( (m+2)/3 ) ! maps (1,2,3)->1, (4,5,6)->2, ...
-    CS%restore_rate(m) = 1.0 / ( timescale_in_days(ig) * 86400.0 )
+    CS%restore_rate(m) = 1.0 / ( timescale_in_days(ig) * 86400.0*US%s_to_T )
   enddo
 
   CS%tr_Reg => tr_Reg
@@ -125,8 +127,8 @@ subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS)
   type(nw2_tracers_CS),               pointer    :: CS !< The control structure returned by a previous
                                                        !! call to register_nw2_tracer.
   ! Local variables
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights
-  real :: rscl ! z* scaling factor
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights [Z ~> m]
+  real :: rscl ! z* scaling factor [nondim]
   character(len=8)  :: var_name ! The variable's name.
   integer :: i, j, k, m
 
@@ -206,11 +208,11 @@ subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
 ! The arguments to this subroutine are redundant in that
 !     h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
   ! Local variables
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights [Z ~> m]
   integer :: i, j, k, m
-  real :: dt_x_rate ! dt * restoring rate
-  real :: rscl ! z* scaling factor
+  real :: dt_x_rate ! dt * restoring rate [nondim]
+  real :: rscl ! z* scaling factor [nondim]
   real :: target_value ! tracer value
 
 ! if (.not.associated(CS)) return
@@ -253,8 +255,8 @@ subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
   endif
 
   do m=1,CS%ntr
-    dt_x_rate = ( dt * CS%restore_rate(m) ) * US%T_to_s
-!$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate)
+    dt_x_rate = dt * CS%restore_rate(m)
+    !$OMP parallel do default(shared) private(target_value)
     do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
       target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k)
       CS%tr(i,j,k,m) = CS%tr(i,j,k,m) + G%mask2dT(i,j) * dt_x_rate * ( target_value - CS%tr(i,j,k,m) )
@@ -270,13 +272,13 @@ real function nw2_tracer_dist(m, G, GV, eta, i, j, k)
   type(ocean_grid_type),   intent(in) :: G   !< The ocean's grid structure
   type(verticalGrid_type), intent(in) :: GV  !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),0:SZK_(G)), &
-                           intent(in) :: eta !< Interface position [m]
+                           intent(in) :: eta !< Interface position [Z ~> m]
   integer, intent(in) :: i !< Cell index i
   integer, intent(in) :: j !< Cell index j
   integer, intent(in) :: k !< Layer index k
   ! Local variables
-  real :: pi ! 3.1415...
-  real :: x, y, z ! non-dimensional positions
+  real :: pi ! 3.1415... [nondim]
+  real :: x, y, z ! non-dimensional relative positions [nondim]
   pi = 2.*acos(0.)
   x = ( G%geolonT(i,j) - G%west_lon ) / G%len_lon ! 0 ... 1
   y = -G%geolatT(i,j) / G%south_lat ! -1 ... 1
diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90
index 862209a688..6f690ab760 100644
--- a/src/tracer/oil_tracer.F90
+++ b/src/tracer/oil_tracer.F90
@@ -51,7 +51,6 @@ module oil_tracer
   real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this subroutine, in g m-3?
   real, dimension(NTR_MAX) :: IC_val = 0.0    !< The (uniform) initial condition value.
   real, dimension(NTR_MAX) :: land_val = -1.0 !< The value of tr used where land is masked out.
-  real, dimension(NTR_MAX) :: oil_decay_days  !< Decay time scale of oil [days]
   real, dimension(NTR_MAX) :: oil_decay_rate  !< Decay rate of oil [T-1 ~> s-1] calculated from oil_decay_days
   integer, dimension(NTR_MAX) :: oil_source_k !< Layer of source
   logical :: oil_may_reinit  !< If true, oil tracers may be reset by the initialization code
@@ -83,7 +82,8 @@ function register_oil_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
   ! Local variables
   character(len=40)  :: mdl = "oil_tracer" ! This module's name.
 ! This include declares and sets the variable "version".
-#include "version_variable.h"
+# include "version_variable.h"
+  real, dimension(NTR_MAX) :: oil_decay_days  !< Decay time scale of oil [days]
   character(len=200) :: inputdir ! The directory where the input files are.
   character(len=48)  :: var_name ! The variable's name.
   character(len=3)   :: name_tag ! String for creating identifying oils
@@ -136,7 +136,7 @@ function register_oil_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
   call get_param(param_file, mdl, "OIL_SOURCE_RATE", CS%oil_source_rate, &
                  "The rate of oil injection.", &
                  units="kg s-1", scale=US%T_to_s, default=1.0)
-  call get_param(param_file, mdl, "OIL_DECAY_DAYS", CS%oil_decay_days, &
+  call get_param(param_file, mdl, "OIL_DECAY_DAYS", oil_decay_days, &
                  "The decay timescale in days (if positive), or no decay "//&
                  "if 0, or use the temperature dependent decay rate of "//&
                  "Adcroft et al. (GRL, 2010) if negative.", units="days", &
@@ -156,9 +156,9 @@ function register_oil_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
       CS%ntr = CS%ntr + 1
       CS%tr_desc(m) = var_desc("oil"//trim(name_tag), "kg m-3", "Oil Tracer", caller=mdl)
       CS%IC_val(m) = 0.0
-      if (CS%oil_decay_days(m)>0.) then
-        CS%oil_decay_rate(m) = 1. / (86400.0*US%s_to_T * CS%oil_decay_days(m))
-      elseif (CS%oil_decay_days(m)<0.) then
+      if (oil_decay_days(m) > 0.) then
+        CS%oil_decay_rate(m) = 1. / (86400.0*US%s_to_T * oil_decay_days(m))
+      elseif (oil_decay_days(m) < 0.) then
         CS%oil_decay_rate(m) = -1.
       endif
     endif
@@ -326,9 +326,12 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
 
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
-  real :: Isecs_per_year = 1.0 / (365.0*86400.0)
+  real :: Isecs_per_year = 1.0 / (365.0*86400.0) ! Conversion factor from seconds to year [year s-1]
   real :: vol_scale ! A conversion factor for volumes into m3 [m3 H-1 L-2 ~> 1 or m3 kg-1]
-  real :: year, h_total, ldecay
+  real :: year      ! Time in fractional years [years]
+  real :: h_total   ! A running sum of thicknesses [H ~> m or kg m-2]
+  real :: decay_timescale ! Chemical decay timescale for oil [T ~> s]
+  real :: ldecay    ! Chemical decay rate of oil [T-1 ~> s-1]
   integer :: i, j, k, is, ie, js, je, nz, m, k_max
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
@@ -360,8 +363,8 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
       if (CS%oil_decay_rate(m)>0.) then
         CS%tr(i,j,k,m) = G%mask2dT(i,j)*max(1. - dt*CS%oil_decay_rate(m),0.)*CS%tr(i,j,k,m) ! Safest
       elseif (CS%oil_decay_rate(m)<0.) then
-        ldecay = 12.*(3.0**(-(tv%T(i,j,k)-20.)/10.)) ! Timescale [days]
-        ldecay = 1. / (86400.*US%s_to_T * ldecay) ! Rate [T-1 ~> s-1]
+        decay_timescale = (12.*(3.0**(-(tv%T(i,j,k)-20.)/10.))) * (86400.*US%s_to_T) ! Timescale [s ~> T]
+        ldecay = 1. / decay_timescale ! Rate [T-1 ~> s-1]
         CS%tr(i,j,k,m) = G%mask2dT(i,j)*max(1. - dt*ldecay,0.)*CS%tr(i,j,k,m)
       endif
     enddo ; enddo ; enddo
@@ -399,12 +402,13 @@ end subroutine oil_tracer_column_physics
 
 !> Calculate the mass-weighted integral of the oil tracer stocks, returning the number of stocks it
 !! has calculated.  If the stock_index is present, only the stock corresponding to that coded index is returned.
-function oil_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function oil_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h    !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),                 intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                               !! tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),              intent(in)    :: US   !< A dimensional unit scaling type
   type(oil_tracer_CS),                pointer       :: CS   !< The control structure returned by a previous
                                                             !! call to register_oil_tracer.
   character(len=*), dimension(:),     intent(out)   :: names  !< the names of the stocks calculated.
@@ -429,7 +433,7 @@ function oil_stock(h, stocks, G, GV, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="oil_stock")
     units(m) = trim(units(m))//" kg"
diff --git a/src/tracer/pseudo_salt_tracer.F90 b/src/tracer/pseudo_salt_tracer.F90
index 94ee126a59..c441e519be 100644
--- a/src/tracer/pseudo_salt_tracer.F90
+++ b/src/tracer/pseudo_salt_tracer.F90
@@ -253,12 +253,13 @@ end subroutine pseudo_salt_tracer_column_physics
 
 !> Calculates the mass-weighted integral of all tracer stocks, returning the number of stocks it has
 !! calculated.  If the stock_index is present, only the stock corresponding to that coded index is returned.
-function pseudo_salt_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function pseudo_salt_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),              intent(in)    :: G      !< The ocean's grid structure
   type(verticalGrid_type),            intent(in)    :: GV     !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h  !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),                 intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                               !! tracer, in kg times concentration units [kg conc]
+  type(unit_scale_type),              intent(in)    :: US     !< A dimensional unit scaling type
   type(pseudo_salt_tracer_CS),        pointer       :: CS     !< The control structure returned by a previous
                                                               !! call to register_pseudo_salt_tracer
   character(len=*), dimension(:),     intent(out)   :: names  !< The names of the stocks calculated
@@ -284,7 +285,7 @@ function pseudo_salt_stock(h, stocks, G, GV, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   call query_vardesc(CS%tr_desc, name=names(1), units=units(1), caller="pseudo_salt_stock")
   units(1) = trim(units(1))//" kg"
   stocks(1) = 0.0
diff --git a/src/tracer/tracer_example.F90 b/src/tracer/tracer_example.F90
index 10551ea247..a41f0ab76d 100644
--- a/src/tracer/tracer_example.F90
+++ b/src/tracer/tracer_example.F90
@@ -358,13 +358,14 @@ end subroutine tracer_column_physics
 !> This function calculates the mass-weighted integral of all tracer stocks,
 !! returning the number of stocks it has calculated.  If the stock_index
 !! is present, only the stock corresponding to that coded index is returned.
-function USER_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index)
+function USER_tracer_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                       intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(:),                 intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                               !! tracer, in kg times concentration units [kg conc].
+  type(unit_scale_type),              intent(in)    :: US     !< A dimensional unit scaling type
   type(USER_tracer_example_CS),       pointer       :: CS     !< The control structure returned by a
                                                               !! previous call to register_USER_tracer.
   character(len=*), dimension(:),     intent(out)   :: names  !< The names of the stocks calculated.
@@ -389,7 +390,7 @@ function USER_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2
+  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,NTR
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="USER_tracer_stock")
     units(m) = trim(units(m))//" kg"
diff --git a/src/user/dyed_channel_initialization.F90 b/src/user/dyed_channel_initialization.F90
index 317ed4ac21..ff98f16529 100644
--- a/src/user/dyed_channel_initialization.F90
+++ b/src/user/dyed_channel_initialization.F90
@@ -133,13 +133,14 @@ subroutine dyed_channel_set_OBC_tracer_data(OBC, G, GV, param_file, tr_Reg)
 end subroutine dyed_channel_set_OBC_tracer_data
 
 !> This subroutine updates the long-channel flow
-subroutine dyed_channel_update_flow(OBC, CS, G, GV, Time)
+subroutine dyed_channel_update_flow(OBC, CS, G, GV, US, Time)
   type(ocean_OBC_type),       pointer    :: OBC !< This open boundary condition type specifies
                                                 !! whether, where, and what open boundary
                                                 !! conditions are used.
   type(dyed_channel_OBC_CS),  pointer    :: CS  !< Dyed channel control structure.
   type(ocean_grid_type),      intent(in) :: G   !< The ocean's grid structure.
   type(verticalGrid_type),    intent(in) :: GV  !< The ocean's vertical grid structure.
+  type(unit_scale_type),      intent(in) :: US  !< A dimensional unit scaling type
   type(time_type),            intent(in) :: Time !< model time.
   ! Local variables
   character(len=40)  :: mdl = "dyed_channel_update_flow" ! This subroutine's name.
@@ -154,7 +155,7 @@ subroutine dyed_channel_update_flow(OBC, CS, G, GV, Time)
   if (.not.associated(OBC)) call MOM_error(FATAL, 'dyed_channel_initialization.F90: '// &
         'dyed_channel_update_flow() was called but OBC type was not initialized!')
 
-  time_sec = G%US%s_to_T * time_type_to_real(Time)
+  time_sec = US%s_to_T * time_type_to_real(Time)
   PI = 4.0*atan(1.0)
 
   do l=1, OBC%number_of_segments
diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90
index 2c84a6040c..840f0bf3ed 100644
--- a/src/user/shelfwave_initialization.F90
+++ b/src/user/shelfwave_initialization.F90
@@ -158,7 +158,7 @@ subroutine shelfwave_set_OBC_data(OBC, CS, G, GV, US, h, Time)
   time_sec = US%s_to_T*time_type_to_real(Time)
   omega = CS%omega
   alpha = CS%alpha
-  my_amp = 1.0*G%US%m_s_to_L_T
+  my_amp = 1.0*US%m_s_to_L_T
   jj = CS%jj
   kk = CS%kk
   ll = CS%ll
diff --git a/src/user/supercritical_initialization.F90 b/src/user/supercritical_initialization.F90
index 12a31f3a75..b4ceb1905d 100644
--- a/src/user/supercritical_initialization.F90
+++ b/src/user/supercritical_initialization.F90
@@ -8,8 +8,9 @@ module supercritical_initialization
 use MOM_file_parser,    only : get_param, log_version, param_file_type
 use MOM_grid,           only : ocean_grid_type
 use MOM_open_boundary,  only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE, OBC_segment_type
-use MOM_verticalGrid,   only : verticalGrid_type
 use MOM_time_manager,   only : time_type, time_type_to_real
+use MOM_unit_scaling,   only : unit_scale_type
+use MOM_verticalGrid,   only : verticalGrid_type
 
 implicit none ; private
 
@@ -17,18 +18,16 @@ module supercritical_initialization
 
 public supercritical_set_OBC_data
 
-! This include declares and sets the variable "version".
-#include "version_variable.h"
-
 contains
 
 !> This subroutine sets the properties of flow at open boundary conditions.
-subroutine supercritical_set_OBC_data(OBC, G, GV, param_file)
+subroutine supercritical_set_OBC_data(OBC, G, GV, US, param_file)
   type(ocean_OBC_type),    pointer    :: OBC  !< This open boundary condition type specifies
                                               !! whether, where, and what open boundary
                                               !! conditions are used.
   type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
   type(verticalGrid_type), intent(in) :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
   type(param_file_type),   intent(in) :: param_file !< Parameter file structure
   ! Local variables
   character(len=40)  :: mdl = "supercritical_set_OBC_data" ! This subroutine's name.
@@ -42,7 +41,7 @@ subroutine supercritical_set_OBC_data(OBC, G, GV, param_file)
 
   call get_param(param_file, mdl, "SUPERCRITICAL_ZONAL_FLOW", zonal_flow, &
                  "Constant zonal flow imposed at upstream open boundary.", &
-                 units="m/s", default=8.57, scale=G%US%m_s_to_L_T)
+                 units="m/s", default=8.57, scale=US%m_s_to_L_T)
 
   do l=1, OBC%number_of_segments
     segment => OBC%segment(l)
diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90
index 51772e2f9f..2438b4115a 100644
--- a/src/user/tidal_bay_initialization.F90
+++ b/src/user/tidal_bay_initialization.F90
@@ -51,13 +51,14 @@ function register_tidal_bay_OBC(param_file, CS, US, OBC_Reg)
 end function register_tidal_bay_OBC
 
 !> This subroutine sets the properties of flow at open boundary conditions.
-subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, h, Time)
+subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time)
   type(ocean_OBC_type),    pointer    :: OBC  !< This open boundary condition type specifies
                                               !! whether, where, and what open boundary
                                               !! conditions are used.
   type(tidal_bay_OBC_CS),  intent(in) :: CS   !< tidal bay control structure.
   type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
   type(verticalGrid_type), intent(in) :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
   type(time_type),         intent(in) :: Time !< model time.
 
@@ -84,7 +85,7 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, h, Time)
 
   allocate(my_area(1:1,js:je))
 
-  flux_scale = GV%H_to_m*G%US%L_to_m
+  flux_scale = GV%H_to_m*US%L_to_m
 
   time_sec = time_type_to_real(Time)
   cff_eta = 0.1*GV%m_to_H * sin(2.0*PI*time_sec/(12.0*3600.0))
@@ -108,7 +109,7 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, h, Time)
 
     if (.not. segment%on_pe) cycle
 
-    segment%normal_vel_bt(:,:) = my_flux / (G%US%m_to_Z*G%US%m_to_L*total_area)
+    segment%normal_vel_bt(:,:) = my_flux / (US%m_to_Z*US%m_to_L*total_area)
     segment%eta(:,:) = cff_eta
 
   enddo ! end segment loop

From 9cb93044ef06a929171b5bd326a7c5c06aa72316 Mon Sep 17 00:00:00 2001
From: jiandewang <jiande.wang@noaa.gov>
Date: Mon, 20 Dec 2021 12:47:19 -0500
Subject: [PATCH 4/7] EMC stochastic candidate 20211028 (#1538)

The stochastic physics feature has been added in MOM6. The following are from Phil Pegion:

The ocean stochastic physics has been re-coded such that there is a wrapper in config_src/external/OCEAN_stochastic_phyiscs that contains the calls to the external stochastic_physics repository. This has been added to support non-UFS applications of MOM6 where the stochastic_physics repository is not part of the build. The init and run procedures are called from src/core/MOM.F90. I have also created a new control structure stochastic_CS, which contains the logical variables, and random patterns which are then passed into src/parameterizations/vertical/MOM_diabadic_driver.F90 and src/parameterizations/vertical/MOM_energetic-PBL.F90.

The writing of the ocean stochastic restarts sit in config_src/nuopc_cap/mom_cap.F90

Co-authored-by: pjpegion <Philip.Pegion@noaa.gov>
---
 config_src/drivers/nuopc_cap/mom_cap.F90      |  10 +-
 .../nuopc_cap/mom_ocean_model_nuopc.F90       |  26 +++-
 .../stochastic_physics/stochastic_physics.F90 |  68 +++++++++
 src/core/MOM.F90                              |  10 +-
 .../stochastic/MOM_stochastics.F90            | 144 ++++++++++++++++++
 .../vertical/MOM_diabatic_driver.F90          |  70 ++++++++-
 .../vertical/MOM_energetic_PBL.F90            |  45 ++++--
 7 files changed, 350 insertions(+), 23 deletions(-)
 create mode 100644 config_src/external/stochastic_physics/stochastic_physics.F90
 create mode 100644 src/parameterizations/stochastic/MOM_stochastics.F90

diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90
index ee498f4184..652f9e5b47 100644
--- a/config_src/drivers/nuopc_cap/mom_cap.F90
+++ b/config_src/drivers/nuopc_cap/mom_cap.F90
@@ -97,6 +97,7 @@ module MOM_cap_mod
 use NUOPC_Model, only: model_label_SetRunClock    => label_SetRunClock
 use NUOPC_Model, only: model_label_Finalize       => label_Finalize
 use NUOPC_Model, only: SetVM
+
 !$use omp_lib             , only : omp_set_num_threads
 
 implicit none; private
@@ -1524,7 +1525,7 @@ subroutine ModelAdvance(gcomp, rc)
   integer                                :: nc
   type(ESMF_Time)                        :: MyTime
   integer                                :: seconds, day, year, month, hour, minute
-  character(ESMF_MAXSTR)                 :: restartname, cvalue
+  character(ESMF_MAXSTR)                 :: restartname, cvalue, stoch_restartname
   character(240)                         :: msgString
   character(ESMF_MAXSTR)                 :: casename
   integer                                :: iostat
@@ -1738,14 +1739,19 @@ subroutine ModelAdvance(gcomp, rc)
         ! write the final restart without a timestamp
         if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then
            write(restartname,'(A)')"MOM.res"
+           write(stoch_restartname,'(A)')"ocn_stoch.res.nc"
         else
            write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2)') &
                 "MOM.res.", year, month, day, hour, minute, seconds
+           write(stoch_restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') &
+                "ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc"
         endif
         call ESMF_LogWrite("MOM_cap: Writing restart :  "//trim(restartname), ESMF_LOGMSG_INFO)
 
         ! write restart file(s)
-        call ocean_model_restart(ocean_state, restartname=restartname)
+        call ocean_model_restart(ocean_state, restartname=restartname, &
+                                stoch_restartname=stoch_restartname)
+
      endif
 
      if (is_root_pe()) then
diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
index aab909e56e..448f23140e 100644
--- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
@@ -62,6 +62,7 @@ module MOM_ocean_model_nuopc
 use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum
 use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS
 use MOM_surface_forcing_nuopc, only : forcing_save_restart
+use get_stochy_pattern_mod,  only : write_stoch_restart_ocn
 use iso_fortran_env,           only : int64
 
 #include <MOM_memory.h>
@@ -176,6 +177,10 @@ module MOM_ocean_model_nuopc
                               !! steps can span multiple coupled time steps.
   logical :: diabatic_first   !< If true, apply diabatic and thermodynamic
                               !! processes before time stepping the dynamics.
+  logical :: do_sppt         !< If true, stochastically perturb the diabatic and
+                             !! write restarts
+  logical :: pert_epbl       !< If true, then randomly perturb the KE dissipation and
+                             !! genration termsand write restarts
 
   real :: eps_omesh           !< Max allowable difference between ESMF mesh and MOM6
                               !! domain coordinates
@@ -425,6 +430,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   endif
 
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
+! get number of processors and PE list for stocasthci physics initialization
+  call get_param(param_file, mdl, "DO_SPPT", OS%do_sppt, &
+                 "If true, then stochastically perturb the thermodynamic "//&
+                 "tendencies of T,S, and h.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
+  call get_param(param_file, mdl, "PERT_EPBL", OS%pert_epbl, &
+                 "If true, then stochastically perturb the kinetic energy "//&
+                 "production and dissipation terms.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
 
   call close_param_file(param_file)
   call diag_mediator_close_registration(OS%diag)
@@ -686,7 +702,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 end subroutine update_ocean_model
 
 !> This subroutine writes out the ocean model restart file.
-subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files)
+subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, num_rest_files)
   type(ocean_state_type),     pointer    :: OS !< A pointer to the structure containing the
                                                !! internal ocean state being saved to a restart file
   character(len=*), optional, intent(in) :: timestamp !< An optional timestamp string that should be
@@ -694,6 +710,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files)
   character(len=*), optional, intent(in) :: restartname !< Name of restart file to use
                                                !! This option distinguishes the cesm interface from the
                                                !! non-cesm interface
+  character(len=*), optional, intent(in) :: stoch_restartname !< Name of restart file to use
+                                               !! This option distinguishes the cesm interface from the
+                                               !! non-cesm interface
   integer, optional, intent(out)         :: num_rest_files !< number of restart files written
 
   if (.not.MOM_state_is_synchronized(OS%MOM_CSp)) &
@@ -733,6 +752,11 @@ subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files)
         endif
      endif
   endif
+  if (present(stoch_restartname)) then
+      if (OS%do_sppt .OR. OS%pert_epbl) then
+         call write_stoch_restart_ocn('RESTART/'//trim(stoch_restartname))
+     endif
+  endif
 
 end subroutine ocean_model_restart
 ! </SUBROUTINE> NAME="ocean_model_restart"
diff --git a/config_src/external/stochastic_physics/stochastic_physics.F90 b/config_src/external/stochastic_physics/stochastic_physics.F90
new file mode 100644
index 0000000000..df62aa1591
--- /dev/null
+++ b/config_src/external/stochastic_physics/stochastic_physics.F90
@@ -0,0 +1,68 @@
+! The are stubs for ocean stochastic physics
+! the fully functional code is available at
+! http://github.com/noaa-psd/stochastic_physics
+module stochastic_physics
+
+implicit none
+
+private
+
+public :: init_stochastic_physics_ocn
+public :: run_stochastic_physics_ocn
+
+contains
+
+!!!!!!!!!!!!!!!!!!!!
+subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_in,do_sppt_in, &
+                                       mpiroot, mpicomm, iret)
+implicit none
+real,intent(in)       :: delt !< timestep in seconds between calls to run_stochastic_physics_ocn
+integer,intent(in)    :: nx   !< number of gridpoints in the x-direction of the compute grid
+integer,intent(in)    :: ny   !< number of gridpoints in the y-direction of the compute grid
+integer,intent(in)    :: nz   !< number of gridpoints in the z-direction of the compute grid
+real,intent(in)       :: geoLonT(nx,ny) !< Longitude in degrees
+real,intent(in)       :: geoLatT(nx,ny) !< Latitude in degrees
+logical,intent(in)    :: pert_epbl_in !< logical flag, if true generate random pattern for ePBL perturbations
+logical,intent(in)    :: do_sppt_in   !< logical flag, if true generate random pattern for SPPT perturbations
+integer,intent(in)    :: mpiroot !< root processor
+integer,intent(in)    :: mpicomm !< mpi communicator
+integer, intent(out)  :: iret    !< return code
+
+iret=0
+if (pert_epbl_in .EQV. .true. ) then
+   print*,'pert_epbl needs to be false if using the stub'
+   iret=-1
+endif
+if (do_sppt_in.EQV. .true. ) then
+   print*,'do_sppt needs to be false if using the stub'
+   iret=-1
+endif
+return
+end subroutine init_stochastic_physics_ocn
+
+subroutine run_stochastic_physics_ocn(sppt_wts,t_rp1,t_rp2)
+implicit none
+real, intent(inout) :: sppt_wts(:,:) !< array containing random weights for SPPT range [0,2]
+real, intent(inout) :: t_rp1(:,:)    !< array containing random weights for ePBL
+                                     !! perturbations (KE generation) range [0,2]
+real, intent(inout) :: t_rp2(:,:)    !< array containing random weights for ePBL
+                                     !! perturbations (KE dissipation) range [0,2]
+return
+end subroutine run_stochastic_physics_ocn
+
+end module stochastic_physics
+
+module get_stochy_pattern_mod
+
+private
+
+public  :: write_stoch_restart_ocn
+
+contains
+subroutine write_stoch_restart_ocn(sfile)
+
+character(len=*) :: sfile   !< name of restart file
+return
+end subroutine write_stoch_restart_ocn
+
+end module get_stochy_pattern_mod
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index 1f1492f2e5..4072cf54a0 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -59,6 +59,7 @@ module MOM
 use MOM_coord_initialization,  only : MOM_initialize_coord
 use MOM_diabatic_driver,       only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member
 use MOM_diabatic_driver,       only : adiabatic, adiabatic_driver_init, diabatic_driver_end
+use MOM_stochastics,           only : stochastics_init, update_stochastics, stochastic_CS
 use MOM_diagnostics,           only : calculate_diagnostic_fields, MOM_diagnostics_init
 use MOM_diagnostics,           only : register_transport_diags, post_transport_diagnostics
 use MOM_diagnostics,           only : register_surface_diags, write_static_fields
@@ -396,6 +397,7 @@ module MOM
   type(ODA_CS), pointer :: odaCS => NULL() !< a pointer to the control structure for handling
                                 !! ensemble model state vectors and data assimilation
                                 !! increments and priors
+  type(stochastic_CS), pointer :: stoch_CS => NULL() !< a pointer to the stochastics control structure
 end type MOM_control_struct
 
 public initialize_MOM, finish_MOM_initialization, MOM_end
@@ -640,6 +642,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
       call disable_averaging(CS%diag)
     endif
   endif
+  ! advance the random pattern if stochastic physics is active
+  if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl) call update_stochastics(CS%stoch_CS)
 
   if (do_dyn) then
     if (G%nonblocking_updates) &
@@ -790,6 +794,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
         enddo ; enddo
       endif
 
+
       call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, &
                              dt_therm_here, bbl_time_int, CS, &
                              Time_local, Waves=Waves)
@@ -1342,7 +1347,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
     call cpu_clock_begin(id_clock_diabatic)
 
     call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, dtdia, &
-                  Time_end_thermo, G, GV, US, CS%diabatic_CSp, OBC=CS%OBC, Waves=Waves)
+                  Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS,OBC=CS%OBC, Waves=Waves)
     fluxes%fluxes_used = .true.
 
     if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)")
@@ -2834,6 +2839,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     call init_oda(Time, G, GV, CS%diag, CS%odaCS)
   endif
 
+  ! initialize stochastic physics
+  call stochastics_init(CS%dt_therm, CS%G, CS%GV, CS%stoch_CS, param_file, diag, Time)
+
   !### This could perhaps go here instead of in finish_MOM_initialization?
   ! call fix_restart_scaling(GV)
   ! call fix_restart_unit_scaling(US)
diff --git a/src/parameterizations/stochastic/MOM_stochastics.F90 b/src/parameterizations/stochastic/MOM_stochastics.F90
new file mode 100644
index 0000000000..21a22a222e
--- /dev/null
+++ b/src/parameterizations/stochastic/MOM_stochastics.F90
@@ -0,0 +1,144 @@
+!> Top-level module for the MOM6 ocean model in coupled mode.
+module MOM_stochastics
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+! This is the top level module for the MOM6 ocean model.  It contains routines
+! for initialization, update, and writing restart of stochastic physics. This
+! particular version wraps all of the calls for MOM6 in the calls that had
+! been used for MOM4.
+!
+use MOM_diag_mediator,       only : register_diag_field, diag_ctrl, time_type
+use MOM_grid,                only : ocean_grid_type
+use MOM_verticalGrid,        only : verticalGrid_type
+use MOM_error_handler,       only : MOM_error, FATAL, WARNING, is_root_pe
+use MOM_error_handler,       only : callTree_enter, callTree_leave
+use MOM_file_parser,         only : get_param, log_version, close_param_file, param_file_type
+use mpp_domains_mod,         only : domain2d, mpp_get_layout, mpp_get_global_domain
+use mpp_domains_mod,         only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
+use MOM_domains,             only : root_PE,num_PEs
+use MOM_coms,                only : Get_PElist
+use stochastic_physics,      only : init_stochastic_physics_ocn, run_stochastic_physics_ocn
+
+#include <MOM_memory.h>
+
+implicit none ; private
+
+public stochastics_init, update_stochastics
+
+!> This control structure holds parameters for the MOM_stochastics module
+type, public:: stochastic_CS
+  logical :: do_sppt         !< If true, stochastically perturb the diabatic
+  logical :: pert_epbl       !< If true, then randomly perturb the KE dissipation and genration terms
+  integer :: id_sppt_wts  = -1 !< Diagnostic id for SPPT
+  integer :: id_epbl1_wts=-1 !< Diagnostic id for epbl generation perturbation
+  integer :: id_epbl2_wts=-1 !< Diagnostic id for epbl dissipation perturbation
+  ! stochastic patterns
+  real, allocatable :: sppt_wts(:,:)  !< Random pattern for ocean SPPT
+                                     !! tendencies with a number between 0 and 2
+  real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation
+  real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation
+  type(diag_ctrl), pointer :: diag   !< structure used to regulate timing of diagnostic output
+  type(time_type), pointer :: Time !< Pointer to model time (needed for sponges)
+end type stochastic_CS
+
+contains
+
+!!   This subroutine initializes the stochastics physics control structure.
+subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time)
+  real, intent(in)                     :: dt       !< time step [T ~> s]
+  type(ocean_grid_type),   intent(in)  :: grid     !< horizontal grid information
+  type(verticalGrid_type), intent(in)  :: GV       !< vertical grid structure
+  type(stochastic_CS), pointer,     intent(inout):: CS !< stochastic control structure
+  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
+  type(diag_ctrl), target, intent(inout) :: diag             !< structure to regulate diagnostic output
+  type(time_type), target                :: Time             !< model time
+  ! Local variables
+  integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean
+  integer :: mom_comm          ! list of pes for this instance of the ocean
+  integer :: num_procs         ! number of processors to pass to stochastic physics
+  integer :: iret              ! return code from stochastic physics
+  integer :: me                !  my pe
+  integer :: pe_zero           !  root pe
+  integer :: nx                ! number of x-points including halo
+  integer :: ny                ! number of x-points including halo
+
+! This include declares and sets the variable "version".
+#include "version_variable.h"
+  character(len=40)  :: mdl = "ocean_stochastics_init"  ! This module's name.
+
+  call callTree_enter("ocean_model_stochastic_init(), MOM_stochastics.F90")
+  if (associated(CS)) then
+    call MOM_error(WARNING, "MOM_stochastics_init called with an "// &
+                            "associated control structure.")
+    return
+  else ; allocate(CS) ; endif
+
+  CS%diag => diag
+  CS%Time => Time
+
+  ! Read all relevant parameters and write them to the model log.
+  call log_version(param_file, mdl, version, "")
+
+! get number of processors and PE list for stocasthci physics initialization
+  call get_param(param_file, mdl, "DO_SPPT", CS%do_sppt, &
+                 "If true, then stochastically perturb the thermodynamic "//&
+                 "tendemcies of T,S, amd h.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
+  call get_param(param_file, mdl, "PERT_EPBL", CS%pert_epbl, &
+                 "If true, then stochastically perturb the kinetic energy "//&
+                 "production and dissipation terms.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
+  if (CS%do_sppt .OR. CS%pert_epbl) then
+     num_procs=num_PEs()
+     allocate(pelist(num_procs))
+     call Get_PElist(pelist,commID = mom_comm)
+     pe_zero=root_PE()
+     nx = grid%ied - grid%isd + 1
+     ny = grid%jed - grid%jsd + 1
+     call init_stochastic_physics_ocn(dt,grid%geoLonT,grid%geoLatT,nx,ny,GV%ke, &
+                                      CS%pert_epbl,CS%do_sppt,pe_zero,mom_comm,iret)
+     if (iret/=0)  then
+         call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed")
+         return
+     endif
+
+     if (CS%do_sppt) allocate(CS%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
+     if (CS%pert_epbl) then
+       allocate(CS%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
+       allocate(CS%epbl2_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
+     endif
+  endif
+  CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesT1, Time, &
+       'random pattern for sppt', 'None')
+  CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesT1, Time, &
+      'random pattern for KE generation', 'None')
+  CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesT1, Time, &
+      'random pattern for KE dissipation', 'None')
+
+  if (is_root_pe()) &
+    write(*,'(/12x,a/)') '=== COMPLETED MOM STOCHASTIC INITIALIZATION ====='
+
+  call callTree_leave("ocean_model_init(")
+  return
+end subroutine stochastics_init
+
+!> update_ocean_model uses the forcing in Ice_ocean_boundary to advance the
+!! ocean model's state from the input value of Ocean_state (which must be for
+!! time time_start_update) for a time interval of Ocean_coupling_time_step,
+!! returning the publicly visible ocean surface properties in Ocean_sfc and
+!! storing the new ocean properties in Ocean_state.
+subroutine update_stochastics(CS)
+  type(stochastic_CS),      intent(inout) :: CS        !< diabatic control structure
+  call callTree_enter("update_stochastics(), MOM_stochastics.F90")
+
+! update stochastic physics patterns before running next time-step
+  call run_stochastic_physics_ocn(CS%sppt_wts,CS%epbl1_wts,CS%epbl2_wts)
+
+  return
+end subroutine update_stochastics
+
+end module MOM_stochastics
+
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index a546bcdec0..c9df559583 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -68,6 +68,7 @@ module MOM_diabatic_driver
 use MOM_verticalGrid,        only : verticalGrid_type, get_thickness_units
 use MOM_wave_speed,          only : wave_speeds
 use MOM_wave_interface,      only : wave_parameters_CS
+use MOM_stochastics,         only : stochastic_CS
 
 implicit none ; private
 
@@ -263,7 +264,7 @@ module MOM_diabatic_driver
 !>  This subroutine imposes the diapycnal mass fluxes and the
 !!  accompanying diapycnal advection of momentum and tracers.
 subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
-                    G, GV, US, CS, OBC, Waves)
+                    G, GV, US, CS, stoch_CS, OBC, Waves)
   type(ocean_grid_type),                      intent(inout) :: G        !< ocean grid structure
   type(verticalGrid_type),                    intent(in)    :: GV       !< ocean vertical grid structure
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u        !< zonal velocity [L T-1 ~> m s-1]
@@ -283,6 +284,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
   type(time_type),                            intent(in)    :: Time_end !< Time at the end of the interval
   type(unit_scale_type),                      intent(in)    :: US       !< A dimensional unit scaling type
   type(diabatic_CS),                          pointer       :: CS       !< module control structure
+  type(stochastic_CS),                        pointer       :: stoch_CS !< stochastic control structure
   type(ocean_OBC_type),             optional, pointer       :: OBC      !< Open boundaries control structure.
   type(Wave_parameters_CS),         optional, pointer       :: Waves    !< Surface gravity waves
 
@@ -295,6 +297,28 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
   integer :: i, j, k, m, is, ie, js, je, nz
   logical :: showCallTree ! If true, show the call tree
 
+  real, allocatable, dimension(:,:,:)    :: h_in  ! thickness before thermodynamics
+  real, allocatable, dimension(:,:,:)    :: t_in  ! temperature before thermodynamics
+  real, allocatable, dimension(:,:,:)    :: s_in  ! salinity before thermodynamics
+  real :: t_tend,s_tend,h_tend                  ! holder for tendencey needed for SPPT
+  real :: t_pert,s_pert,h_pert                  ! holder for perturbations needed for SPPT
+
+  if (G%ke == 1) return
+
+   ! save  copy of the date for SPPT if active
+  if (stoch_CS%do_sppt) then
+    allocate(h_in(G%isd:G%ied, G%jsd:G%jed,G%ke))
+    allocate(t_in(G%isd:G%ied, G%jsd:G%jed,G%ke))
+    allocate(s_in(G%isd:G%ied, G%jsd:G%jed,G%ke))
+    h_in(:,:,:)=h(:,:,:)
+    t_in(:,:,:)=tv%T(:,:,:)
+    s_in(:,:,:)=tv%S(:,:,:)
+
+    if (stoch_CS%id_sppt_wts > 0) then
+      call post_data(stoch_CS%id_sppt_wts, stoch_CS%sppt_wts, CS%diag)
+    endif
+  endif
+
   if (GV%ke == 1) return
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
@@ -378,10 +402,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
 
   if (CS%useALEalgorithm .and. CS%use_legacy_diabatic) then
     call diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
-                      G, GV, US, CS, Waves)
+                      G, GV, US, CS, stoch_CS, Waves)
   elseif (CS%useALEalgorithm) then
     call diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
-                      G, GV, US, CS, Waves)
+                      G, GV, US, CS, stoch_CS, Waves)
   else
     call layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
                           G, GV, US, CS, Waves)
@@ -447,13 +471,41 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
 
   if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US)
 
+  if (stoch_CS%do_sppt) then
+  ! perturb diabatic tendecies
+    do k=1,nz
+      do j=js,je
+        do i=is,ie
+          h_tend = (h(i,j,k)-h_in(i,j,k))*stoch_CS%sppt_wts(i,j)
+          t_tend = (tv%T(i,j,k)-t_in(i,j,k))*stoch_CS%sppt_wts(i,j)
+          s_tend = (tv%S(i,j,k)-s_in(i,j,k))*stoch_CS%sppt_wts(i,j)
+          h_pert=h_tend+h_in(i,j,k)
+          t_pert=t_tend+t_in(i,j,k)
+          s_pert=s_tend+s_in(i,j,k)
+          if (h_pert > GV%Angstrom_H) then
+            h(i,j,k) = h_pert
+          else
+            h(i,j,k) = GV%Angstrom_H
+          endif
+          tv%T(i,j,k) = t_pert
+          if (s_pert > 0.0) then
+            tv%S(i,j,k) = s_pert
+          endif
+        enddo
+      enddo
+    enddo
+    deallocate(h_in)
+    deallocate(t_in)
+    deallocate(s_in)
+  endif
+
 end subroutine diabatic
 
 
 !> Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use
 !! with an ALE algorithm.  This version uses an older set of algorithms compared with diabatic_ALE.
 subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
-                           G, GV, US, CS, Waves)
+                           G, GV, US, CS, stoch_CS, Waves)
   type(ocean_grid_type),                      intent(inout) :: G        !< ocean grid structure
   type(verticalGrid_type),                    intent(in)    :: GV       !< ocean vertical grid structure
   type(unit_scale_type),                      intent(in)    :: US       !< A dimensional unit scaling type
@@ -473,6 +525,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
   real,                                       intent(in)    :: dt       !< time increment [T ~> s]
   type(time_type),                            intent(in)    :: Time_end !< Time at the end of the interval
   type(diabatic_CS),                          pointer       :: CS       !< module control structure
+  type(stochastic_CS),                        pointer       :: stoch_CS !< stochastic control structure
   type(Wave_parameters_CS),         optional, pointer       :: Waves    !< Surface gravity waves
 
   ! local variables
@@ -774,7 +827,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
     call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
-                       CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
+                      CS%energetic_PBL_CSp, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, &
+                      waves=waves)
 
     if (associated(Hml)) then
       call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US)
@@ -1037,7 +1091,7 @@ end subroutine diabatic_ALE_legacy
 !>  This subroutine imposes the diapycnal mass fluxes and the
 !!  accompanying diapycnal advection of momentum and tracers.
 subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
-                        G, GV, US, CS, Waves)
+                        G, GV, US, CS, stoch_CS, Waves)
   type(ocean_grid_type),                      intent(inout) :: G        !< ocean grid structure
   type(verticalGrid_type),                    intent(in)    :: GV       !< ocean vertical grid structure
   type(unit_scale_type),                      intent(in)    :: US       !< A dimensional unit scaling type
@@ -1057,6 +1111,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
   real,                                       intent(in)    :: dt       !< time increment [T ~> s]
   type(time_type),                            intent(in)    :: Time_end !< Time at the end of the interval
   type(diabatic_CS),                          pointer       :: CS       !< module control structure
+  type(stochastic_CS),                        pointer       :: stoch_CS !< stochastic control structure
   type(Wave_parameters_CS),         optional, pointer       :: Waves    !< Surface gravity waves
 
   ! local variables
@@ -1309,7 +1364,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
     call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
-                       CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
+                       CS%energetic_PBL_CSp, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, &
+                       waves=waves)
 
     if (associated(Hml)) then
       call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US)
diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90
index 5a9e67bfd9..351f55984f 100644
--- a/src/parameterizations/vertical/MOM_energetic_PBL.F90
+++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90
@@ -17,6 +17,7 @@ module MOM_energetic_PBL
 use MOM_variables, only : thermo_var_ptrs
 use MOM_verticalGrid, only : verticalGrid_type
 use MOM_wave_interface, only: wave_parameters_CS, Get_Langmuir_Number
+use MOM_stochastics,         only : stochastic_CS
 
 implicit none ; private
 
@@ -168,7 +169,6 @@ module MOM_energetic_PBL
 
   real, allocatable, dimension(:,:) :: &
     ML_depth            !< The mixed layer depth determined by active mixing in ePBL [Z ~> m].
-
   ! These are terms in the mixed layer TKE budget, all in [R Z3 T-3 ~> W m-2 = kg s-3].
   real, allocatable, dimension(:,:) :: &
     diag_TKE_wind, &   !< The wind source of TKE [R Z3 T-3 ~> W m-2].
@@ -245,9 +245,9 @@ module MOM_energetic_PBL
 !!  mixed layer model.  It assumes that heating, cooling and freshwater fluxes
 !!  have already been applied.  All calculations are done implicitly, and there
 !!  is no stability limit on the time step.
-subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS, &
-                         dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, last_call, &
-                         dT_expected, dS_expected, Waves )
+subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS,  &
+                         stoch_CS, dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, &
+                         last_call, dT_expected, dS_expected, Waves)
   type(ocean_grid_type),   intent(inout) :: G      !< The ocean's grid structure.
   type(verticalGrid_type), intent(in)    :: GV     !< The ocean's vertical grid structure.
   type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
@@ -300,6 +300,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
                                                    !! diffusivities are applied [ppt].
   type(wave_parameters_CS), &
                  optional, pointer       :: Waves  !< Wave CS
+  type(stochastic_CS),  pointer       :: stoch_CS  !< The control structure returned by a previous
 
 !    This subroutine determines the diffusivities from the integrated energetics
 !  mixed layer model.  It assumes that heating, cooling and freshwater fluxes
@@ -456,11 +457,17 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
       ! Perhaps provide a first guess for MLD based on a stored previous value.
       MLD_io = -1.0
       if (CS%MLD_iteration_guess .and. (CS%ML_Depth(i,j) > 0.0))  MLD_io = CS%ML_Depth(i,j)
-
-      call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
-                       u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, &
-                       US, CS, eCD, dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j)
-
+      if (stoch_CS%pert_epbl) then ! stochastics are active
+        call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
+                         u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, &
+                         US, CS, eCD, dt_diag=dt_diag, Waves=Waves, G=G, &
+                         epbl1_wt=stoch_CS%epbl1_wts(i,j),epbl2_wt=stoch_CS%epbl2_wts(i,j), &
+                         i=i, j=j)
+      else
+        call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
+                         u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, &
+                         US, CS, eCD, dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j)
+      endif
 
       ! Copy the diffusivities to a 2-d array.
       do K=1,nz+1
@@ -531,6 +538,12 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
     if (CS%id_LA > 0)       call post_data(CS%id_LA, CS%LA, CS%diag)
     if (CS%id_LA_MOD > 0)   call post_data(CS%id_LA_MOD, CS%LA_MOD, CS%diag)
     if (CS%id_MSTAR_LT > 0) call post_data(CS%id_MSTAR_LT, CS%MSTAR_LT, CS%diag)
+    ! only write random patterns if running with stochastic physics, otherwise the
+    ! array is unallocated and will give an error
+    if  (stoch_CS%pert_epbl) then
+      if (stoch_CS%id_epbl1_wts > 0)    call post_data(stoch_CS%id_epbl1_wts, stoch_CS%epbl1_wts, CS%diag)
+      if (stoch_CS%id_epbl2_wts > 0)    call post_data(stoch_CS%id_epbl2_wts, stoch_CS%epbl2_wts, CS%diag)
+    endif
   endif
 
   if (debug) deallocate(eCD%dT_expect, eCD%dS_expect)
@@ -543,7 +556,7 @@ end subroutine energetic_PBL
 !!  mixed layer model for a single column of water.
 subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, &
                        u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, &
-                       dt_diag, Waves, G, i, j)
+                       dt_diag, Waves, G, epbl1_wt, epbl2_wt, i, j)
   type(verticalGrid_type), intent(in)    :: GV     !< The ocean's vertical grid structure.
   type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
   real, dimension(SZK_(GV)), intent(in)  :: h      !< Layer thicknesses [H ~> m or kg m-2].
@@ -587,6 +600,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
                  optional, pointer       :: Waves  !< Wave CS for Langmuir turbulence
   type(ocean_grid_type), &
                  optional, intent(inout) :: G      !< The ocean's grid structure.
+  real,          optional, intent(in)    :: epbl1_wt !< random number to perturb KE generation
+  real,          optional, intent(in)    :: epbl2_wt !< random number to perturb KE dissipation
   integer,       optional, intent(in)    :: i      !< The i-index to work on (used for Waves)
   integer,       optional, intent(in)    :: j      !< The i-index to work on (used for Waves)
 
@@ -881,6 +896,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
       else
         mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3)
       endif
+      ! stochastically pertrub mech_TKE in the UFS
+      if (present(epbl1_wt)) mech_TKE=mech_TKE*epbl1_wt
 
       if (CS%TKE_diagnostics) then
         eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0
@@ -963,7 +980,12 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
         if (Idecay_len_TKE > 0.0) exp_kh = exp(-h(k-1)*Idecay_len_TKE)
         if (CS%TKE_diagnostics) &
           eCD%dTKE_mech_decay = eCD%dTKE_mech_decay + (exp_kh-1.0) * mech_TKE * I_dtdiag
-        mech_TKE = mech_TKE * exp_kh
+        if (present(epbl2_wt)) then ! perturb the TKE destruction
+           mech_TKE = mech_TKE * (1+(exp_kh-1) * epbl2_wt)
+        else
+           mech_TKE = mech_TKE * exp_kh
+        endif
+        !if ( i .eq. 10 .and. j .eq. 10 .and. k .eq. nz) print*,'mech TKE', mech_TKE
 
         !   Accumulate any convectively released potential energy to contribute
         ! to wstar and to drive penetrating convection.
@@ -2373,7 +2395,6 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS)
       Time, 'Velocity Scale that is used.', 'm s-1', conversion=US%Z_to_m*US%s_to_T)
   CS%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, &
       Time, 'Total mstar that is used.', 'nondim')
-
   if (CS%use_LT) then
     CS%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, &
         Time, 'Langmuir number.', 'nondim')

From 986bc8c0b9391222fe84c95c1db3fe8d1743d342 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Mon, 13 Dec 2021 20:19:21 -0500
Subject: [PATCH 5/7] Corrected the unit documentation for 31 variables

  Corrected the documentation of the units for 31 variables in various modules.
All answers and output are bitwise identical.
---
 src/core/MOM_PressureForce_FV.F90                |  3 ++-
 src/core/MOM_continuity_PPM.F90                  |  4 ++--
 src/core/MOM_dynamics_split_RK2.F90              |  4 ++--
 src/diagnostics/MOM_wave_structure.F90           | 16 ++++++++--------
 src/initialization/MOM_shared_initialization.F90 |  2 +-
 .../vertical/MOM_bulk_mixed_layer.F90            |  4 ++--
 .../vertical/MOM_diabatic_aux.F90                |  2 +-
 .../vertical/MOM_set_viscosity.F90               |  4 ++--
 src/tracer/MOM_lateral_boundary_diffusion.F90    |  2 +-
 src/tracer/MOM_tracer_registry.F90               |  3 ++-
 src/user/Idealized_Hurricane.F90                 |  2 +-
 src/user/MOM_controlled_forcing.F90              |  5 +++--
 src/user/SCM_CVMix_tests.F90                     |  6 +++---
 src/user/baroclinic_zone_initialization.F90      |  2 +-
 14 files changed, 31 insertions(+), 28 deletions(-)

diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90
index 5ead019717..1666b4a97e 100644
--- a/src/core/MOM_PressureForce_FV.F90
+++ b/src/core/MOM_PressureForce_FV.F90
@@ -153,7 +153,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
                         ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1].
   real :: H_to_RL2_T2   ! A factor to convert from thickness units (H) to pressure
                         ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].
-!  real :: oneatm = 101325.0  ! 1 atm in [Pa] = [kg m-1 s-2]
+!  real :: oneatm       ! 1 standard atmosphere of pressure in [R L2 T-2 ~> Pa]
   real, parameter :: C1_6 = 1.0/6.0
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
@@ -187,6 +187,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
       p(i,j,1) = p_atm(i,j)
     enddo ; enddo
   else
+    ! oneatm = 101325.0 * US%kg_m3_to_R * US%m_s_to_L_T**2 ! 1 atm scaled to [R L2 T-2 ~> Pa]
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       p(i,j,1) = 0.0 ! or oneatm
diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90
index 17d2f830c0..7bc67e2fdf 100644
--- a/src/core/MOM_continuity_PPM.F90
+++ b/src/core/MOM_continuity_PPM.F90
@@ -924,8 +924,8 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0,
     FAmt_0, &     ! test velocities [H L ~> m2 or kg m-1].
     uhtot_L, &    ! The summed transport with the westerly (uhtot_L) and
     uhtot_R       ! and easterly (uhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
-  real :: FA_0    ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m].
-  real :: FA_avg  ! The average effective face area [L H ~> m2 or kg m], nominally given by
+  real :: FA_0    ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m-1].
+  real :: FA_avg  ! The average effective face area [L H ~> m2 or kg m-1], nominally given by
                   ! the realized transport divided by the barotropic velocity.
   real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]. This
                        ! limiting is necessary to keep the inverse of visc_rem
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index a762da7f33..f1a2b2469d 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -337,8 +337,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
 
   real, pointer, dimension(:,:,:) :: &
     ! These pointers are used to alter which fields are passed to btstep with various options:
-    u_ptr => NULL(), &   ! A pointer to a zonal velocity [L T-1]
-    v_ptr => NULL(), &   ! A pointer to a meridional velocity [L T-1]
+    u_ptr => NULL(), &   ! A pointer to a zonal velocity [L T-1 ~> m s-1]
+    v_ptr => NULL(), &   ! A pointer to a meridional velocity [L T-1 ~> m s-1]
     uh_ptr => NULL(), &  ! A pointer to a zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
     vh_ptr => NULL(), &  ! A pointer to a meridional volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
     ! These pointers are just used as shorthand for CS%u_av, CS%v_av, and CS%h_av.
diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90
index 2dd272d409..c833e973c5 100644
--- a/src/diagnostics/MOM_wave_structure.F90
+++ b/src/diagnostics/MOM_wave_structure.F90
@@ -171,11 +171,11 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
   real, dimension(SZK_(GV))   :: dz           !< thicknesses of merged layers (same as Hc I hope) [Z ~> m]
   ! real, dimension(SZK_(GV)+1) :: dWdz_profile !< profile of dW/dz
   real                        :: w2avg   !< average of squared vertical velocity structure funtion [Z ~> m]
-  real                        :: int_dwdz2
-  real                        :: int_w2
-  real                        :: int_N2w2
-  real                        :: KE_term !< terms in vertically averaged energy equation
-  real                        :: PE_term !< terms in vertically averaged energy equation
+  real                        :: int_dwdz2 !< Vertical integral of the square of u_strct [Z ~> m]
+  real                        :: int_w2   !< Vertical integral of the square of w_strct [Z ~> m]
+  real                        :: int_N2w2 !< Vertical integral of N2 [Z T-2 ~> m s-2]
+  real                        :: KE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2]
+  real                        :: PE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2]
   real                        :: W0      !< A vertical velocity magnitude [Z T-1 ~> m s-1]
   real                        :: gp_unscaled !< A version of gprime rescaled to [L T-2 ~> m s-2].
   real, dimension(SZK_(GV)-1) :: lam_z   !< product of eigen value and gprime(k); one value for each
@@ -183,8 +183,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
   real, dimension(SZK_(GV)-1) :: a_diag, b_diag, c_diag
                                          !< diagonals of tridiagonal matrix; one value for each
                                          !< interface (excluding surface and bottom)
-  real, dimension(SZK_(GV)-1) :: e_guess !< guess at eigen vector with unit amplitde (for TDMA)
-  real, dimension(SZK_(GV)-1) :: e_itt   !< improved guess at eigen vector (from TDMA)
+  real, dimension(SZK_(GV)-1) :: e_guess !< guess at eigen vector with unit amplitude (for TDMA) [nondim]
+  real, dimension(SZK_(GV)-1) :: e_itt   !< improved guess at eigen vector (from TDMA) [nondim]
   real    :: Pi
   integer :: kc
   integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop
@@ -523,7 +523,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
 
             ! Back-calculate amplitude from energy equation
             if (present(En) .and. (freq**2*Kmag2 > 0.0)) then
-              ! Units here are [R
+              ! Units here are [R Z ~> kg m-2]
               KE_term = 0.25*GV%Rho0*( ((freq**2 + f2) / (freq**2*Kmag2))*int_dwdz2 + int_w2 )
               PE_term = 0.25*GV%Rho0*( int_N2w2 / freq**2 )
               if (En(i,j) >= 0.0) then
diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90
index bb5a84033b..fc5ceaf3e4 100644
--- a/src/initialization/MOM_shared_initialization.F90
+++ b/src/initialization/MOM_shared_initialization.F90
@@ -814,7 +814,7 @@ subroutine reset_face_lengths_list(G, param_file, US)
   real, allocatable, dimension(:) :: &
     Dmin_u, Dmax_u, Davg_u   ! Porous barrier monomial fit params [m]
   real, allocatable, dimension(:) :: &
-    Dmin_v, Dmax_v, Davg_v
+    Dmin_v, Dmax_v, Davg_v   ! Porous barrier monomial fit params [m]
   real    :: lat, lon     ! The latitude and longitude of a point.
   real    :: len_lon      ! The periodic range of longitudes, usually 360 degrees.
   real    :: len_lat      ! The range of latitudes, usually 180 degrees.
diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
index 046329523d..dd160c300c 100644
--- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
+++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
@@ -489,7 +489,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
     ! netMassInOut = water [H ~> m or kg m-2] added/removed via surface fluxes
     ! netMassOut   = water [H ~> m or kg m-2] removed via evaporating surface fluxes
     ! net_heat     = heat via surface fluxes [degC H ~> degC m or degC kg m-2]
-    ! net_salt     = salt via surface fluxes [ppt H ~> dppt m or gSalt m-2]
+    ! net_salt     = salt via surface fluxes [ppt H ~> ppt m or gSalt m-2]
     ! Pen_SW_bnd   = components to penetrative shortwave radiation
     call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
                   CS%H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
@@ -1527,7 +1527,7 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, &
   real :: TKE_full_ent  ! The TKE remaining if a layer is fully entrained
                         ! [Z L2 T-2 ~> m3 s-2].
   real :: dRL       ! Work required to mix water from the next layer
-                    ! across the mixed layer [L2 T-2 ~> L2 s-2].
+                    ! across the mixed layer [L2 T-2 ~> m2 s-2].
   real :: Pen_En_Contrib  ! Penetrating SW contributions to the changes in
                           ! TKE, divided by layer thickness in m [L2 T-2 ~> m2 s-2].
   real :: Cpen1     ! A temporary variable [L2 T-2 ~> m2 s-2].
diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90
index 312d114dde..c421b3a0f7 100644
--- a/src/parameterizations/vertical/MOM_diabatic_aux.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90
@@ -1171,7 +1171,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
     !                netMassOut < 0 means mass leaves ocean.
     ! netHeat      = heat via surface fluxes [degC H ~> degC m or degC kg m-2], excluding the part
     !                contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
-    ! netSalt      = surface salt fluxes [ppt H ~> dppt m or gSalt m-2]
+    ! netSalt      = surface salt fluxes [ppt H ~> ppt m or gSalt m-2]
     ! Pen_SW_bnd   = components to penetrative shortwave radiation split according to bands.
     !                This field provides that portion of SW from atmosphere that in fact
     !                enters to the ocean and participates in pentrative SW heating.
diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90
index 861a8957c1..f7d4b0cc0d 100644
--- a/src/parameterizations/vertical/MOM_set_viscosity.F90
+++ b/src/parameterizations/vertical/MOM_set_viscosity.F90
@@ -200,7 +200,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
   real :: hwtot            ! Sum of the thicknesses used to calculate
                            ! the near-bottom velocity magnitude [H ~> m or kg m-2].
   real :: hutot            ! Running sum of thicknesses times the
-                           ! velocity magnitudes [H T T-1 ~> m2 s-1 or kg m-1 s-1].
+                           ! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
   real :: Thtot            ! Running sum of thickness times temperature [degC H ~> degC m or degC kg m-2].
   real :: Shtot            ! Running sum of thickness times salinity [ppt H ~> ppt m or ppt kg m-2].
   real :: hweight          ! The thickness of a layer that is within Hbbl
@@ -597,7 +597,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
     ! When stratification dominates h_N<<h_f, and vice versa.
     do i=is,ie ; if (do_i(i)) then
       ! The 400.0 in this expression is the square of a Ci introduced in KW99, eq. 2.22.
-      ustarsq = Rho0x400_G * ustar(i)**2 ! Note not in units of u*^2 but [H R ~> kg m-2 or kg m-5]
+      ustarsq = Rho0x400_G * ustar(i)**2 ! Note not in units of u*^2 but [H R ~> kg m-2 or kg2 m-5]
       htot = 0.0
 
       ! Calculate the thickness of a stratification limited BBL ignoring rotation:
diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90
index 19d40f2db1..c11bc9856c 100644
--- a/src/tracer/MOM_lateral_boundary_diffusion.F90
+++ b/src/tracer/MOM_lateral_boundary_diffusion.F90
@@ -596,7 +596,7 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_
   real, allocatable :: F_layer_z(:)  !< Diffusive flux at U/V-point in the ztop grid    [H L2 conc ~> m3 conc]
   real              :: h_vel(ke)     !< Thicknesses at u- and v-points in the native grid
                                      !! The harmonic mean is used to avoid zero values      [H ~> m or kg m-2]
-  real    :: khtr_avg                !< Thickness-weighted diffusivity at the velocity-point [L2 T-1 ~> m s-1]
+  real    :: khtr_avg                !< Thickness-weighted diffusivity at the velocity-point [L2 T-1 ~> m2 s-1]
                                      !! This is just to remind developers that khtr_avg should be
                                      !! computed once khtr is 3D.
   real    :: htot                    !< Total column thickness                              [H ~> m or kg m-2]
diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90
index c729231927..2c77df3e74 100644
--- a/src/tracer/MOM_tracer_registry.F90
+++ b/src/tracer/MOM_tracer_registry.F90
@@ -84,7 +84,8 @@ module MOM_tracer_registry
 !  real, dimension(:,:,:), pointer :: diff_cont_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
 !                                                              !! [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1]
 !  real, dimension(:,:,:), pointer :: diff_conc_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
-!                                                              !! expressed as a change in concentration [conc T-1]
+!                                                              !! expressed as a change in concentration
+!                                                              !! [conc T-1 ~> conc s-1]
   real, dimension(:,:,:), pointer :: t_prev         => NULL() !< tracer concentration array at a previous
                                                               !! timestep used for diagnostics [conc]
   real, dimension(:,:,:), pointer :: Trxh_prev      => NULL() !< layer integrated tracer concentration array
diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90
index 7182fc364a..707a0972f9 100644
--- a/src/user/Idealized_Hurricane.F90
+++ b/src/user/Idealized_Hurricane.F90
@@ -595,7 +595,7 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C
   U_TS = CS%hurr_translation_spd*0.5*cos(transdir)
   V_TS = CS%hurr_translation_spd*0.5*sin(transdir)
 
-  ! Set the surface wind stresses, in [Pa]. A positive taux
+  ! Set the surface wind stresses, in [R L Z T-2 ~> Pa]. A positive taux
   ! accelerates the ocean to the (pseudo-)east.
   !   The i-loop extends to is-1 so that taux can be used later in the
   ! calculation of ustar - otherwise the lower bound would be Isq.
diff --git a/src/user/MOM_controlled_forcing.F90 b/src/user/MOM_controlled_forcing.F90
index f783a271a6..d136d58a19 100644
--- a/src/user/MOM_controlled_forcing.F90
+++ b/src/user/MOM_controlled_forcing.F90
@@ -46,7 +46,7 @@ module MOM_controlled_forcing
   real    :: lam_prec       !< A constant of proportionality between SSS anomalies
                             !! (normalised by mean SSS) and precipitation [R Z T-1 ~> kg m-2 s-1]
   real    :: lam_cyc_heat   !< A constant of proportionality between cyclical SST
-                            !! anomalies and corrective heat fluxes [W m-2 degC-1]
+                            !! anomalies and corrective heat fluxes [Q R Z T-1 degC-1 ~> W m-2 degC-1]
   real    :: lam_cyc_prec   !< A constant of proportionality between cyclical SSS
                             !! anomalies (normalised by mean SSS) and corrective
                             !! precipitation [R Z T-1 ~> kg m-2 s-1]
@@ -270,7 +270,8 @@ subroutine apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, virt_heat, virt_prec
     ! Accumulate the average anomalies for this period.
     dt_wt = wt_per1 * dt
     CS%avg_time(m_mid) = CS%avg_time(m_mid) + dt_wt
-    ! These loops temporarily change the units of the CS%avg_ variables to [degC s] or [ppt s].
+    ! These loops temporarily change the units of the CS%avg_ variables to [degC T ~> degC s]
+    ! or [ppt T ~> ppt s].
     do j=js,je ; do i=is,ie
       CS%avg_SST_anom(i,j,m_mid) = CS%avg_SST_anom(i,j,m_mid) + &
                                    dt_wt * G%mask2dT(i,j) * SST_anom(i,j)
diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90
index 1fbc7a2b62..5bbe65b8d8 100644
--- a/src/user/SCM_CVMix_tests.F90
+++ b/src/user/SCM_CVMix_tests.F90
@@ -36,8 +36,8 @@ module SCM_CVMix_tests
   logical :: UseHeatFlux    !< True to use heat flux
   logical :: UseEvaporation !< True to use evaporation
   logical :: UseDiurnalSW   !< True to use diurnal sw radiation
-  real :: tau_x !< (Constant) Wind stress, X [Pa]
-  real :: tau_y !< (Constant) Wind stress, Y [Pa]
+  real :: tau_x !< (Constant) Wind stress, X [R L Z T-2 ~> Pa]
+  real :: tau_y !< (Constant) Wind stress, Y [R L Z T-2 ~> Pa]
   real :: surf_HF !< (Constant) Heat flux [degC Z T-1 ~> m degC s-1]
   real :: surf_evap !< (Constant) Evaporation rate [Z T-1 ~> m s-1]
   real :: Max_sw !< maximum of diurnal sw radiation [degC Z T-1 ~> degC m s-1]
@@ -56,7 +56,7 @@ subroutine SCM_CVMix_tests_TS_init(T, S, h, G, GV, US, param_file, just_read)
   type(ocean_grid_type),                     intent(in)  :: G  !< Grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< Vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T  !< Potential temperature [degC]
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [psu]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [ppt]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h  !< Layer thickness [H ~> m or kg m-2]
   type(unit_scale_type),                     intent(in)  :: US !< A dimensional unit scaling type
   type(param_file_type),                     intent(in)  :: param_file !< Input parameter structure
diff --git a/src/user/baroclinic_zone_initialization.F90 b/src/user/baroclinic_zone_initialization.F90
index 1555f4ecad..a214012541 100644
--- a/src/user/baroclinic_zone_initialization.F90
+++ b/src/user/baroclinic_zone_initialization.F90
@@ -36,7 +36,7 @@ subroutine bcz_params(G, GV, US, param_file, S_ref, dSdz, delta_S, dSdx, T_ref,
   real,                    intent(out) :: S_ref      !< Reference salinity [ppt]
   real,                    intent(out) :: dSdz       !< Salinity stratification [ppt Z-1 ~> ppt m-1]
   real,                    intent(out) :: delta_S    !< Salinity difference across baroclinic zone [ppt]
-  real,                    intent(out) :: dSdx       !< Linear salinity gradient [ppt m-1]
+  real,                    intent(out) :: dSdx       !< Linear salinity gradient [ppt G%xaxis_units-1]
   real,                    intent(out) :: T_ref      !< Reference temperature [degC]
   real,                    intent(out) :: dTdz       !< Temperature stratification [degC Z-1 ~> degC m-1]
   real,                    intent(out) :: delta_T    !< Temperature difference across baroclinic zone [degC]

From d2442461abdaf32cb9d737d3e798fcdea597b239 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Tue, 14 Dec 2021 13:25:16 -0500
Subject: [PATCH 6/7] +Rescale tides and ramp-up times

  Rescaled the dimensions of the tidal amplitudes and frequencies used
internally in calc_tidal_forcing() and ramp-up times used by update_OBC_ramp()
and updateCFLtruncationValue() for nearly complete dimensional consistency
testing.  New unit_scale_type arguments were added to 5 routines, in the case of
calc_tidal_forcing() replacing a previous optional argument that was always
being used.  One overly short internal variable, "N", was renamed "nodelon" to
make its purpose clearer and easier to search for.  All answers are bitwise
identical, but there are changes to the argument lists of 5 routines.
---
 src/core/MOM_PressureForce_FV.F90             |   4 +-
 src/core/MOM_PressureForce_Montgomery.F90     |   4 +-
 src/core/MOM_dynamics_split_RK2.F90           |  10 +-
 src/core/MOM_dynamics_unsplit.F90             |   2 +-
 src/core/MOM_dynamics_unsplit_RK2.F90         |   2 +-
 src/core/MOM_open_boundary.F90                |  25 +--
 .../lateral/MOM_tidal_forcing.F90             | 158 ++++++++++--------
 .../vertical/MOM_vert_friction.F90            |  34 ++--
 8 files changed, 127 insertions(+), 112 deletions(-)

diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90
index 1666b4a97e..30b2e90d1a 100644
--- a/src/core/MOM_PressureForce_FV.F90
+++ b/src/core/MOM_PressureForce_FV.F90
@@ -306,7 +306,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref
     enddo ; enddo
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       za(i,j) = za(i,j) - GV%g_Earth * e_tidal(i,j)
@@ -574,7 +574,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
       enddo ; enddo
     enddo
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
   endif
 
 !    Here layer interface heights, e, are calculated.
diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90
index a827fb12d0..18ea07b313 100644
--- a/src/core/MOM_PressureForce_Montgomery.F90
+++ b/src/core/MOM_PressureForce_Montgomery.F90
@@ -203,7 +203,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
       enddo ; enddo ; enddo
     endif
 
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       geopot_bot(i,j) = -GV%g_Earth*(e_tidal(i,j) + G%bathyT(i,j))
@@ -451,7 +451,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
       enddo ; enddo
     enddo
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
   endif
 
 !    Here layer interface heights, e, are calculated.
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index f1a2b2469d..68b844562f 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -374,7 +374,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   enddo
 
   ! Update CFL truncation value as function of time
-  call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp)
+  call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp, US)
 
   if (CS%debug) then
     call MOM_state_chksum("Start predictor ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
@@ -395,7 +395,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     if (CS%debug_OBC) call open_boundary_test_extern_h(G, GV, CS%OBC, h)
 
     ! Update OBC ramp value as function of time
-    call update_OBC_ramp(Time_local, CS%OBC)
+    call update_OBC_ramp(Time_local, CS%OBC, US)
 
     do k=1,nz ; do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB
       u_old_rad_OBC(I,j,k) = u_av(I,j,k)
@@ -1207,20 +1207,20 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp)
+  if (use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
                           CS%tides_CSp)
   call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp)
   call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, &
                      ntrunc, CS%vertvisc_CSp)
   CS%set_visc_CSp => set_visc
-  call updateCFLtruncationValue(Time, CS%vertvisc_CSp, &
+  call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, &
                                 activate=is_new_run(restart_CS) )
 
   if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp
   if (associated(OBC)) then
     CS%OBC => OBC
-    if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, &
+    if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, &
                                 activate=is_new_run(restart_CS) )
   endif
   if (associated(update_OBC_CSp)) CS%update_OBC_CSp => update_OBC_CSp
diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90
index 88a11e071c..fcc4c3d49b 100644
--- a/src/core/MOM_dynamics_unsplit.F90
+++ b/src/core/MOM_dynamics_unsplit.F90
@@ -666,7 +666,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp)
+  if (use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
                           CS%tides_CSp)
   call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc)
diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90
index 26bd00aaf5..694d88f2ea 100644
--- a/src/core/MOM_dynamics_unsplit_RK2.F90
+++ b/src/core/MOM_dynamics_unsplit_RK2.F90
@@ -628,7 +628,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp)
+  if (use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
                           CS%tides_CSp)
   call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc)
diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90
index 6d8696216a..2c3f016005 100644
--- a/src/core/MOM_open_boundary.F90
+++ b/src/core/MOM_open_boundary.F90
@@ -264,7 +264,7 @@ module MOM_open_boundary
   logical :: add_tide_constituents = .false.          !< If true, add tidal constituents to the boundary elevation
                                                       !! and velocity. Will be set to true if n_tide_constituents > 0.
   character(len=2), allocatable, dimension(:) :: tide_names  !< Names of tidal constituents to add to the boundary data.
-  real, allocatable, dimension(:) :: tide_frequencies        !< Angular frequencies of chosen tidal constituents [s-1].
+  real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal constituents [T-1 ~> s-1].
   real, allocatable, dimension(:) :: tide_eq_phases   !< Equilibrium phases of chosen tidal constituents [rad].
   real, allocatable, dimension(:) :: tide_fn          !< Amplitude modulation of boundary tides by nodal cycle [nondim].
   real, allocatable, dimension(:) :: tide_un          !< Phase modulation of boundary tides by nodal cycle [rad].
@@ -305,8 +305,8 @@ module MOM_open_boundary
                    !! the independence of the OBCs to this external data [L T-1 ~> m s-1].
   logical :: ramp = .false.                 !< If True, ramp from zero to the external values for SSH.
   logical :: ramping_is_activated = .false. !< True if the ramping has been initialized
-  real :: ramp_timescale                    !< If ramp is True, use this timescale for ramping [s].
-  real :: trunc_ramp_time                   !< If ramp is True, time after which ramp is done [s].
+  real :: ramp_timescale                    !< If ramp is True, use this timescale for ramping [T ~> s].
+  real :: trunc_ramp_time                   !< If ramp is True, time after which ramp is done [T ~> s].
   real :: ramp_value                        !< If ramp is True, where we are on the ramp from
                                             !! zero to one [nondim].
   type(time_type) :: ramp_start_time        !< Time when model was started.
@@ -627,7 +627,7 @@ subroutine open_boundary_config(G, US, param_file, OBC)
        "Symmetric memory must be used when using Flather OBCs.")
   ! Need to do this last, because it depends on time_interp_external_init having already been called
   if (OBC%add_tide_constituents) then
-    call initialize_obc_tides(OBC, param_file)
+    call initialize_obc_tides(OBC, US, param_file)
     ! Tide update is done within update_OBC_segment_data, so this should be true if tides are included.
     OBC%update_OBC = .true.
   endif
@@ -948,8 +948,9 @@ subroutine initialize_segment_data(G, OBC, PF)
 
 end subroutine initialize_segment_data
 
-subroutine initialize_obc_tides(OBC, param_file)
+subroutine initialize_obc_tides(OBC, US, param_file)
   type(ocean_OBC_type), intent(inout) :: OBC  !< Open boundary control structure
+  type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
   type(param_file_type), intent(in) :: param_file !< Parameter file handle
   integer, dimension(3) :: tide_ref_date      !< Reference date (t = 0) for tidal forcing (year, month, day).
   integer, dimension(3) :: nodal_ref_date     !< Date to calculate nodal modulation for (year, month, day).
@@ -1022,7 +1023,8 @@ subroutine initialize_obc_tides(OBC, param_file)
         "Frequency of the "//trim(OBC%tide_names(c))//" tidal constituent. "//&
         "This is only used if TIDES and TIDE_"//trim(OBC%tide_names(c))// &
         " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(OBC%tide_names(c))//&
-        " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=tidal_frequency(trim(OBC%tide_names(c))))
+        " is in OBC_TIDE_CONSTITUENTS.", &
+        units="s-1", default=tidal_frequency(trim(OBC%tide_names(c))), scale=US%T_to_s)
 
     ! Find equilibrium phase if needed
     if (OBC%add_eq_phase) then
@@ -3727,7 +3729,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
   real :: tidal_elev  ! Interpolated tidal elevation at the OBC points [m]
   real, allocatable :: normal_trans_bt(:,:) ! barotropic transport [H L2 T-1 ~> m3 s-1]
   integer :: turns    ! Number of index quarter turns
-  real :: time_delta  ! Time since tidal reference date [s]
+  real :: time_delta  ! Time since tidal reference date [T ~> s]
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
@@ -3738,7 +3740,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
 
   if (.not. associated(OBC)) return
 
-  if (OBC%add_tide_constituents) time_delta = time_type_to_real(Time - OBC%time_ref)
+  if (OBC%add_tide_constituents) time_delta = US%s_to_T * time_type_to_real(Time - OBC%time_ref)
 
   do n = 1, OBC%number_of_segments
     segment => OBC%segment(n)
@@ -4336,14 +4338,15 @@ end subroutine update_OBC_segment_data
 !> Update the OBC ramp value as a function of time.
 !! If called with the optional argument activate=.true., record the
 !! value of Time as the beginning of the ramp period.
-subroutine update_OBC_ramp(Time, OBC, activate)
+subroutine update_OBC_ramp(Time, OBC, US, activate)
   type(time_type), target, intent(in)    :: Time     !< Current model time
   type(ocean_OBC_type),    intent(inout) :: OBC      !< Open boundary structure
+  type(unit_scale_type),   intent(in)    :: US       !< A dimensional unit scaling type
   logical, optional,       intent(in)    :: activate !< Specify whether to record the value of
                                                      !! Time as the beginning of the ramp period
 
   ! Local variables
-  real :: deltaTime ! The time since start of ramping [s]
+  real :: deltaTime ! The time since start of ramping [T ~> s]
   real :: wghtA     ! A temporary variable used to set OBC%ramp_value [nondim]
   character(len=12) :: msg
 
@@ -4359,7 +4362,7 @@ subroutine update_OBC_ramp(Time, OBC, activate)
     endif
   endif
   if (.not.OBC%ramping_is_activated) return
-  deltaTime = max( 0., time_type_to_real( Time - OBC%ramp_start_time ) )
+  deltaTime = max( 0., US%s_to_T*time_type_to_real( Time - OBC%ramp_start_time ) )
   if (deltaTime >= OBC%trunc_ramp_time) then
     OBC%ramp_value = 1.0
     OBC%ramp = .false. ! This turns off ramping after this call
diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90
index b8d5c44098..cc4517a473 100644
--- a/src/parameterizations/lateral/MOM_tidal_forcing.F90
+++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90
@@ -11,6 +11,7 @@ module MOM_tidal_forcing
 use MOM_grid,          only : ocean_grid_type
 use MOM_io,            only : field_exists, file_exists, MOM_read_data
 use MOM_time_manager,  only : set_date, time_type, time_type_to_real, operator(-)
+use MOM_unit_scaling,  only : unit_scale_type
 
 implicit none ; private
 
@@ -47,12 +48,12 @@ module MOM_tidal_forcing
                       !! astronomical/equilibrium argument.
   real    :: sal_scalar !< The constant of proportionality between sea surface
                       !! height (really it should be bottom pressure) anomalies
-                      !! and bottom geopotential anomalies.
+                      !! and bottom geopotential anomalies [nondim].
   integer :: nc       !< The number of tidal constituents in use.
   real, dimension(MAX_CONSTITUENTS) :: &
-    freq, &           !< The frequency of a tidal constituent [s-1].
-    phase0, &         !< The phase of a tidal constituent at time 0, in radians.
-    amp, &            !< The amplitude of a tidal constituent at time 0 [m].
+    freq, &           !< The frequency of a tidal constituent [T-1 ~> s-1].
+    phase0, &         !< The phase of a tidal constituent at time 0 [rad].
+    amp, &            !< The amplitude of a tidal constituent at time 0 [Z ~> m].
     love_no           !< The Love number of a tidal constituent at time 0 [nondim].
   integer :: struct(MAX_CONSTITUENTS) !< An encoded spatial structure for each constituent
   character (len=16) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent
@@ -62,13 +63,13 @@ module MOM_tidal_forcing
                                    !! tidal phases at t = 0.
   real, allocatable :: &
     sin_struct(:,:,:), &    !< The sine and cosine based structures that can
-    cos_struct(:,:,:), &    !< be associated with the astronomical forcing.
+    cos_struct(:,:,:), &    !< be associated with the astronomical forcing [nondim].
     cosphasesal(:,:,:), &   !< The cosine and sine of the phase of the
     sinphasesal(:,:,:), &   !< self-attraction and loading amphidromes.
-    ampsal(:,:,:), &        !< The amplitude of the SAL [m].
+    ampsal(:,:,:), &        !< The amplitude of the SAL [Z ~> m].
     cosphase_prev(:,:,:), & !< The cosine and sine of the phase of the
     sinphase_prev(:,:,:), & !< amphidromes in the previous tidal solutions.
-    amp_prev(:,:,:)         !< The amplitude of the previous tidal solution [m].
+    amp_prev(:,:,:)         !< The amplitude of the previous tidal solution [Z ~> m].
 end type tidal_forcing_CS
 
 integer :: id_clock_tides !< CPU clock for tides
@@ -87,8 +88,9 @@ module MOM_tidal_forcing
 subroutine astro_longitudes_init(time_ref, longitudes)
   type(time_type), intent(in) :: time_ref            !> Time to calculate longitudes for.
   type(astro_longitudes), intent(out) :: longitudes  !> Lunar and solar longitudes at time_ref.
-  real :: D, T                                       !> Date offsets
-  real, parameter :: PI = 4.0 * atan(1.0)            !> 3.14159...
+  real :: D                                          !> Time since the reference date [days]
+  real :: T                                          !> Time in Julian centuries [centuries]
+  real, parameter :: PI = 4.0 * atan(1.0)            !> 3.14159... [nondim]
   ! Find date at time_ref in days since 1900-01-01
   D = time_type_to_real(time_ref - set_date(1900, 1, 1)) / (24.0 * 3600.0)
   ! Time since 1900-01-01 in Julian centuries
@@ -176,44 +178,45 @@ end function tidal_frequency
 !> Find amplitude (f) and phase (u) modulation of tidal constituents by the 18.6
 !! year nodal cycle. Values here follow Table I.6 in Kowalik and Luick,
 !! "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019.
-subroutine nodal_fu(constit, N, fn, un)
-  character (len=2), intent(in) :: constit              !> Tidal constituent to find modulation for.
-  real, intent(in) :: N                                 !> Longitude of ascending node [rad].
-                                                        !! Calculate using astro_longitudes_init.
-  real, parameter :: RADIANS = 4.0 * atan(1.0) / 180.0  !> Converts degrees to radians.
-  real, intent(out) :: &
-    fn, & !> Amplitude modulation [nondim]
-    un    !> Phase modulation [rad]
+subroutine nodal_fu(constit, nodelon, fn, un)
+  character (len=2), intent(in)  :: constit !> Tidal constituent to find modulation for.
+  real,              intent(in)  :: nodelon !> Longitude of ascending node [rad], which
+                                            !! can be calculated using astro_longitudes_init.
+  real,              intent(out) :: fn      !> Amplitude modulation [nondim]
+  real,              intent(out) :: un      !> Phase modulation [rad]
+
+  real, parameter :: RADIANS = 4.0 * atan(1.0) / 180.0  !> Converts degrees to radians [nondim]
+
   select case (constit)
     case ("M2")
-      fn = 1.0 - 0.037 * cos(N)
-      un = -2.1 * RADIANS * sin(N)
+      fn = 1.0 - 0.037 * cos(nodelon)
+      un = -2.1 * RADIANS * sin(nodelon)
     case ("S2")
       fn = 1.0  ! Solar S2 has no amplitude modulation.
       un = 0.0  ! S2 has no phase modulation.
     case ("N2")
-      fn = 1.0 - 0.037 * cos(N)
-      un = -2.1 * RADIANS * sin(N)
+      fn = 1.0 - 0.037 * cos(nodelon)
+      un = -2.1 * RADIANS * sin(nodelon)
     case ("K2")
-      fn = 1.024 + 0.286 * cos(N)
-      un = -17.7 * RADIANS * sin(N)
+      fn = 1.024 + 0.286 * cos(nodelon)
+      un = -17.7 * RADIANS * sin(nodelon)
     case ("K1")
-      fn = 1.006 + 0.115 * cos(N)
-      un = -8.9 * RADIANS * sin(N)
+      fn = 1.006 + 0.115 * cos(nodelon)
+      un = -8.9 * RADIANS * sin(nodelon)
     case ("O1")
-      fn = 1.009 + 0.187 * cos(N)
-      un = 10.8 * RADIANS * sin(N)
+      fn = 1.009 + 0.187 * cos(nodelon)
+      un = 10.8 * RADIANS * sin(nodelon)
     case ("P1")
       fn = 1.0  ! P1 has no amplitude modulation.
       un = 0.0  ! P1 has no phase modulation.
     case ("Q1")
-      fn = 1.009 + 0.187 * cos(N)
-      un = 10.8 * RADIANS * sin(N)
+      fn = 1.009 + 0.187 * cos(nodelon)
+      un = 10.8 * RADIANS * sin(nodelon)
     case ("MF")
-      fn = 1.043 + 0.414 * cos(N)
-      un = -23.7 * RADIANS * sin(N)
+      fn = 1.043 + 0.414 * cos(nodelon)
+      un = -23.7 * RADIANS * sin(nodelon)
     case ("MM")
-      fn = 1.0 - 0.130 * cos(N)
+      fn = 1.0 - 0.130 * cos(nodelon)
       un = 0.0  ! MM has no phase modulation.
     case default
       call MOM_error(FATAL, "nodal_fu: unrecognized constituent")
@@ -226,10 +229,11 @@ end subroutine nodal_fu
 !! while fields like the background viscosities are 2-D arrays.
 !! ALLOC is a macro defined in MOM_memory.h for allocate or nothing with
 !! static memory.
-subroutine tidal_forcing_init(Time, G, param_file, CS)
-  type(time_type),       intent(in)    :: Time !< The current model time.
-  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure.
-  type(param_file_type), intent(in)    :: param_file !< A structure to parse for run-time parameters.
+subroutine tidal_forcing_init(Time, G, US, param_file, CS)
+  type(time_type),        intent(in)    :: Time !< The current model time.
+  type(ocean_grid_type),  intent(inout) :: G    !< The ocean's grid structure.
+  type(unit_scale_type),  intent(in)    :: US   !< A dimensional unit scaling type
+  type(param_file_type),  intent(in)    :: param_file !< A structure to parse for run-time parameters.
   type(tidal_forcing_CS), intent(inout) :: CS   !< Tidal forcing control struct
 
   ! Local variables
@@ -237,15 +241,18 @@ subroutine tidal_forcing_init(Time, G, param_file, CS)
     phase, &          ! The phase of some tidal constituent.
     lat_rad, lon_rad  ! Latitudes and longitudes of h-points in radians.
   real :: deg_to_rad
-  real, dimension(MAX_CONSTITUENTS) :: freq_def, phase0_def, amp_def, love_def
+  real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [s-1]
+  real, dimension(MAX_CONSTITUENTS) :: phase0_def ! Default reference phase for each tidal constituent [rad]
+  real, dimension(MAX_CONSTITUENTS) :: amp_def  ! Default amplitude for each tidal constituent [m]
+  real, dimension(MAX_CONSTITUENTS) :: love_def ! Default love number for each constituent [nondim]
   integer, dimension(3) :: tide_ref_date !< Reference date (t = 0) for tidal forcing.
   logical :: use_const  ! True if a constituent is being used.
   logical :: use_M2, use_S2, use_N2, use_K2, use_K1, use_O1, use_P1, use_Q1
   logical :: use_MF, use_MM
   logical :: tides      ! True if a tidal forcing is to be used.
   logical :: FAIL_IF_MISSING = .true.
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   character(len=40)  :: mdl = "MOM_tidal_forcing" ! This module's name.
   character(len=128) :: mesg
   character(len=200) :: tidal_input_files(4*MAX_CONSTITUENTS)
@@ -389,68 +396,68 @@ subroutine tidal_forcing_init(Time, G, param_file, CS)
     endif
     CS%time_ref = set_date(tide_ref_date(1), tide_ref_date(2), tide_ref_date(3))
   endif
-  ! Set the parameters for all components that are in use.
-  ! Initialize reference time for tides and
-  ! find relevant lunar and solar longitudes at the reference time.
+
+  ! Initialize reference time for tides and find relevant lunar and solar
+  ! longitudes at the reference time.
   if (CS%use_eq_phase) call astro_longitudes_init(CS%time_ref, CS%tidal_longitudes)
+
+  ! Set the parameters for all components that are in use.
   c=0
   if (use_M2) then
     c=c+1 ; CS%const_name(c) = "M2" ; CS%struct(c) = 2
-    CS%love_no(c) = 0.693 ; CS%amp(c) = 0.242334
+    CS%love_no(c) = 0.693 ; amp_def(c) = 0.242334 ! Default amplitude in m.
   endif
 
   if (use_S2) then
     c=c+1 ; CS%const_name(c) = "S2" ; CS%struct(c) = 2
-    CS%love_no(c) = 0.693 ; CS%amp(c) = 0.112743
+    CS%love_no(c) = 0.693 ; amp_def(c) = 0.112743 ! Default amplitude in m.
   endif
 
   if (use_N2) then
     c=c+1 ; CS%const_name(c) = "N2" ; CS%struct(c) = 2
-    CS%love_no(c) = 0.693 ; CS%amp(c) = 0.046397
+    CS%love_no(c) = 0.693 ; amp_def(c) = 0.046397 ! Default amplitude in m.
   endif
 
   if (use_K2) then
     c=c+1 ; CS%const_name(c) = "K2" ; CS%struct(c) = 2
-    CS%love_no(c) = 0.693 ; CS%amp(c) = 0.030684
+    CS%love_no(c) = 0.693 ; amp_def(c) = 0.030684 ! Default amplitude in m.
   endif
 
   if (use_K1) then
     c=c+1 ; CS%const_name(c) = "K1" ; CS%struct(c) = 1
-    CS%love_no(c) = 0.736 ; CS%amp(c) = 0.141565
+    CS%love_no(c) = 0.736 ; amp_def(c) = 0.141565 ! Default amplitude in m.
   endif
 
   if (use_O1) then
     c=c+1 ; CS%const_name(c) = "O1" ; CS%struct(c) = 1
-    CS%love_no(c) = 0.695 ; CS%amp(c) = 0.100661
+    CS%love_no(c) = 0.695 ; amp_def(c) = 0.100661 ! Default amplitude in m.
   endif
 
   if (use_P1) then
     c=c+1 ; CS%const_name(c) = "P1" ; CS%struct(c) = 1
-    CS%love_no(c) = 0.706 ; CS%amp(c) = 0.046848
+    CS%love_no(c) = 0.706 ; amp_def(c) = 0.046848 ! Default amplitude in m.
   endif
 
   if (use_Q1) then
     c=c+1 ; CS%const_name(c) = "Q1" ; CS%struct(c) = 1
-    CS%love_no(c) = 0.695 ; CS%amp(c) = 0.019273
+    CS%love_no(c) = 0.695 ; amp_def(c) = 0.019273 ! Default amplitude in m.
   endif
 
   if (use_MF) then
     c=c+1 ; CS%const_name(c) = "MF" ; CS%struct(c) = 3
-    CS%love_no(c) = 0.693 ; CS%amp(c) = 0.042041
+    CS%love_no(c) = 0.693 ; amp_def(c) = 0.042041 ! Default amplitude in m.
   endif
 
   if (use_MM) then
     c=c+1 ; CS%const_name(c) = "MM" ; CS%struct(c) = 3
-    CS%love_no(c) = 0.693 ; CS%amp(c) = 0.022191
+    CS%love_no(c) = 0.693 ; amp_def(c) = 0.022191 ! Default amplitude in m.
   endif
 
   ! Set defaults for all included constituents
   ! and things that can be set by functions
   do c=1,nc
-    CS%freq(c) = tidal_frequency(CS%const_name(c))
-    freq_def(c) = CS%freq(c)
+    freq_def(c) = tidal_frequency(CS%const_name(c))
     love_def(c) = CS%love_no(c)
-    amp_def(c) = CS%amp(c)
     CS%phase0(c) = 0.0
     if (CS%use_eq_phase) then
       phase0_def(c) = eq_phase(CS%const_name(c), CS%tidal_longitudes)
@@ -467,11 +474,11 @@ subroutine tidal_forcing_init(Time, G, param_file, CS)
                    "Frequency of the "//trim(CS%const_name(c))//" tidal constituent. "//&
                    "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// &
                    " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(CS%const_name(c))// &
-                   " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=freq_def(c))
+                   " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=freq_def(c), scale=US%T_to_s)
     call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_AMP", CS%amp(c), &
                    "Amplitude of the "//trim(CS%const_name(c))//" tidal constituent. "//&
                    "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// &
-                   " are true.", units="m", default=amp_def(c))
+                   " are true.", units="m", default=amp_def(c), scale=US%m_to_Z)
     call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_PHASE_T0", CS%phase0(c), &
                    "Phase of the "//trim(CS%const_name(c))//" tidal constituent at time 0. "//&
                    "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// &
@@ -484,8 +491,9 @@ subroutine tidal_forcing_init(Time, G, param_file, CS)
     allocate(CS%ampsal(isd:ied,jsd:jed,nc))
     do c=1,nc
       ! Read variables with names like PHASE_SAL_M2 and AMP_SAL_M2.
-      call find_in_files(tidal_input_files,"PHASE_SAL_"//trim(CS%const_name(c)),phase,G)
-      call find_in_files(tidal_input_files,"AMP_SAL_"//trim(CS%const_name(c)),CS%ampsal(:,:,c),G)
+      call find_in_files(tidal_input_files, "PHASE_SAL_"//trim(CS%const_name(c)), phase, G)
+      call find_in_files(tidal_input_files, "AMP_SAL_"//trim(CS%const_name(c)), CS%ampsal(:,:,c), &
+                         G, scale=US%m_to_Z)
       call pass_var(phase,           G%domain,complete=.false.)
       call pass_var(CS%ampsal(:,:,c),G%domain,complete=.true.)
       do j=js-1,je+1 ; do i=is-1,ie+1
@@ -501,8 +509,9 @@ subroutine tidal_forcing_init(Time, G, param_file, CS)
     allocate(CS%amp_prev(isd:ied,jsd:jed,nc))
     do c=1,nc
       ! Read variables with names like PHASE_PREV_M2 and AMP_PREV_M2.
-      call find_in_files(tidal_input_files,"PHASE_PREV_"//trim(CS%const_name(c)),phase,G)
-      call find_in_files(tidal_input_files,"AMP_PREV_"//trim(CS%const_name(c)),CS%amp_prev(:,:,c),G)
+      call find_in_files(tidal_input_files, "PHASE_PREV_"//trim(CS%const_name(c)), phase, G)
+      call find_in_files(tidal_input_files, "AMP_PREV_"//trim(CS%const_name(c)), CS%amp_prev(:,:,c), &
+                         G, scale=US%m_to_Z)
       call pass_var(phase,             G%domain,complete=.false.)
       call pass_var(CS%amp_prev(:,:,c),G%domain,complete=.true.)
       do j=js-1,je+1 ; do i=is-1,ie+1
@@ -518,18 +527,19 @@ end subroutine tidal_forcing_init
 
 !> This subroutine finds a named variable in a list of files and reads its
 !! values into a domain-decomposed 2-d array
-subroutine find_in_files(filenames, varname, array, G)
+subroutine find_in_files(filenames, varname, array, G, scale)
   character(len=*), dimension(:),   intent(in)  :: filenames !< The names of the files to search for the named variable
   character(len=*),                 intent(in)  :: varname   !< The name of the variable to read
   type(ocean_grid_type),            intent(in)  :: G         !< The ocean's grid structure
   real, dimension(SZI_(G),SZJ_(G)), intent(out) :: array     !< The array to fill with the data
+  real,                   optional, intent(in)  :: scale     !< A factor by which to rescale the array.
   ! Local variables
   integer :: nf
 
   do nf=1,size(filenames)
     if (LEN_TRIM(filenames(nf)) == 0) cycle
     if (field_exists(filenames(nf), varname, MOM_domain=G%Domain)) then
-      call MOM_read_data(filenames(nf), varname, array, G%Domain)
+      call MOM_read_data(filenames(nf), varname, array, G%Domain, scale=scale)
       return
     endif
   enddo
@@ -571,22 +581,22 @@ end subroutine tidal_forcing_sensitivity
 !! height.  For now, eta and eta_tidal are both geopotential heights in depth
 !! units, but probably the input for eta should really be replaced with the
 !! column mass anomalies.
-subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z)
+subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, US, CS)
   type(ocean_grid_type),            intent(in)  :: G         !< The ocean's grid structure.
   type(time_type),                  intent(in)  :: Time      !< The time for the caluculation.
   real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: eta       !< The sea surface height anomaly from
                                                              !! a time-mean geoid [Z ~> m].
   real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_tidal !< The tidal forcing geopotential height
                                                              !! anomalies [Z ~> m].
+  type(unit_scale_type),            intent(in)  :: US        !< A dimensional unit scaling type
   type(tidal_forcing_CS),           intent(in)  :: CS        !< The control structure returned by a
                                                              !! previous call to tidal_forcing_init.
-  real,                             intent(in)  :: m_to_Z    !< A scaling factor from m to the units of eta.
 
   ! Local variables
-  real :: now       ! The relative time in seconds.
-  real :: amp_cosomegat, amp_sinomegat
-  real :: cosomegat, sinomegat
-  real :: eta_prop  ! The nondimenional constant of proportionality beteen eta and eta_tidal.
+  real :: now       ! The relative time compared with the tidal reference [T ~> s]
+  real :: amp_cosomegat, amp_sinomegat ! The tidal amplitudes times the components of phase [Z ~> m]
+  real :: cosomegat, sinomegat ! The components of the phase [nondim]
+  real :: eta_prop  ! The nondimenional constant of proportionality beteen eta and eta_tidal [nondim]
   integer :: i, j, c, m, is, ie, js, je, Isq, Ieq, Jsq, Jeq
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
@@ -598,7 +608,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z)
     return
   endif
 
-  now = time_type_to_real(Time - cs%time_ref)
+  now = US%s_to_T * time_type_to_real(Time - cs%time_ref)
 
   if (CS%USE_SAL_SCALAR .and. CS%USE_PREV_TIDES) then
     eta_prop = 2.0*CS%SAL_SCALAR
@@ -614,8 +624,8 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z)
 
   do c=1,CS%nc
     m = CS%struct(c)
-    amp_cosomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c))
-    amp_sinomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c))
+    amp_cosomegat = CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c))
+    amp_sinomegat = CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c))
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       eta_tidal(i,j) = eta_tidal(i,j) + (amp_cosomegat*CS%cos_struct(i,j,m) + &
                                          amp_sinomegat*CS%sin_struct(i,j,m))
@@ -626,7 +636,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z)
     cosomegat = cos(CS%freq(c)*now)
     sinomegat = sin(CS%freq(c)*now)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-      eta_tidal(i,j) = eta_tidal(i,j) + m_to_Z*CS%ampsal(i,j,c) * &
+      eta_tidal(i,j) = eta_tidal(i,j) + CS%ampsal(i,j,c) * &
            (cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c))
     enddo ; enddo
   enddo ; endif
@@ -635,7 +645,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z)
     cosomegat = cos(CS%freq(c)*now)
     sinomegat = sin(CS%freq(c)*now)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-      eta_tidal(i,j) = eta_tidal(i,j) - m_to_Z*CS%SAL_SCALAR*CS%amp_prev(i,j,c) * &
+      eta_tidal(i,j) = eta_tidal(i,j) - CS%SAL_SCALAR*CS%amp_prev(i,j,c) * &
           (cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c))
     enddo ; enddo
   enddo ; endif
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index adac9e83f4..d384500c3d 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -56,12 +56,12 @@ module MOM_vert_friction
                              !! absolute velocities.
   real    :: CFL_trunc       !< Velocity components will be truncated when they
                              !! are large enough that the corresponding CFL number
-                             !! exceeds this value, nondim.
+                             !! exceeds this value [nondim].
   real    :: CFL_report      !< The value of the CFL number that will cause the
-                             !! accelerations to be reported, nondim.  CFL_report
+                             !! accelerations to be reported [nondim].  CFL_report
                              !! will often equal CFL_trunc.
   real    :: truncRampTime   !< The time-scale over which to ramp up the value of
-                             !! CFL_trunc from CFL_truncS to CFL_truncE
+                             !! CFL_trunc from CFL_truncS to CFL_truncE [T ~> s]
   real    :: CFL_truncS      !< The start value of CFL_trunc
   real    :: CFL_truncE      !< The end/target value of CFL_trunc
   logical :: CFLrampingIsActivated = .false. !< True if the ramping has been initialized
@@ -105,7 +105,7 @@ module MOM_vert_friction
                             !! thickness for viscosity.
   logical :: answers_2018   !< If true, use the order of arithmetic and expressions that recover the
                             !! answers from the end of 2018.  Otherwise, use expressions that do not
-                            !! use an arbitary and hard-coded maximum viscous coupling coefficient
+                            !! use an arbitrary and hard-coded maximum viscous coupling coefficient
                             !! between layers.
   logical :: debug          !< If true, write verbose checksums for debugging purposes.
   integer :: nkml           !< The number of layers in the mixed layer.
@@ -533,7 +533,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   endif
 
   if (associated(ADp%du_dt_str) .and.  associated(ADp%dv_dt_str)) then
-    ! Diagnostics for thickness x wind stress acclerations
+    ! Diagnostics for thickness x wind stress accelerations
     if (CS%id_h_du_dt_str > 0) call post_product_u(CS%id_h_du_dt_str, ADp%du_dt_str, ADp%diag_hu, G, nz, CS%diag)
     if (CS%id_h_dv_dt_str > 0) call post_product_v(CS%id_h_dv_dt_str, ADp%dv_dt_str, ADp%diag_hv, G, nz, CS%diag)
 
@@ -555,11 +555,11 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
   type(vertvisc_type),   intent(in)   :: visc !< Viscosities and bottom drag
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                          intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a
-                                              !! barotopic acceleration that a layer experiences after
+                                              !! barotropic acceleration that a layer experiences after
                                               !! viscosity is applied in the zonal direction [nondim]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                          intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a
-                                              !! barotopic acceleration that a layer experiences after
+                                              !! barotropic acceleration that a layer experiences after
                                               !! viscosity is applied in the meridional direction [nondim]
   real,                  intent(in)    :: dt  !< Time increment [T ~> s]
   type(unit_scale_type), intent(in)    :: US  !< A dimensional unit scaling type
@@ -692,7 +692,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
     a_shelf, &  ! The drag coefficients across interfaces in water columns under
                 ! ice shelves [Z T-1 ~> m s-1].
     z_i         ! An estimate of each interface's height above the bottom,
-                ! normalized by the bottom boundary layer thickness, nondim.
+                ! normalized by the bottom boundary layer thickness [nondim]
   real, dimension(SZIB_(G)) :: &
     kv_bbl, &     ! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
     bbl_thick, &  ! The bottom boundary layer thickness [H ~> m or kg m-2].
@@ -715,10 +715,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
                   ! than Hbbl into the interior.
   real :: topfn   ! A function which goes from 1 at the top to 0 much more
                   ! than Htbl into the interior.
-  real :: z2      ! The distance from the bottom, normalized by Hbbl, nondim.
+  real :: z2      ! The distance from the bottom, normalized by Hbbl [nondim]
   real :: z2_wt   ! A nondimensional (0-1) weight used when calculating z2.
   real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2].
-  real :: a_cpl_max  ! The maximum drag doefficient across interfaces, set so that it will be
+  real :: a_cpl_max  ! The maximum drag coefficient across interfaces, set so that it will be
                      ! representable as a 32-bit float in MKS units  [Z T-1 ~> m s-1]
   real :: h_neglect  ! A thickness that is so small it is usually lost
                      ! in roundoff and can be neglected [H ~> m or kg m-2].
@@ -1193,7 +1193,7 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i,
   h_neglect = GV%H_subroundoff
 
   if (CS%answers_2018) then
-    !   The maximum coupling coefficent was originally introduced to avoid
+    !   The maximum coupling coefficient was originally introduced to avoid
     ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
     ! sets the maximum coupling coefficient increment to 1e10 m per timestep.
     I_amax = (1.0e-10*US%Z_to_m) * dt
@@ -1759,7 +1759,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", CS%truncRampTime, &
                  "The time over which the CFL truncation value is ramped "//&
                  "up at the beginning of the run.", &
-                 units="s", default=0.)
+                 units="s", default=0., scale=US%s_to_T)
   CS%CFL_truncE = CS%CFL_trunc
   call get_param(param_file, mdl, "CFL_TRUNCATE_START", CS%CFL_truncS, &
                  "The start value of the truncation CFL number used when "//&
@@ -1937,14 +1937,16 @@ end subroutine vertvisc_init
 !> Update the CFL truncation value as a function of time.
 !! If called with the optional argument activate=.true., record the
 !! value of Time as the beginning of the ramp period.
-subroutine updateCFLtruncationValue(Time, CS, activate)
+subroutine updateCFLtruncationValue(Time, CS, US, activate)
   type(time_type), target, intent(in)    :: Time     !< Current model time
   type(vertvisc_CS),       pointer       :: CS       !< Vertical viscosity control structure
+  type(unit_scale_type),   intent(in)    :: US       !< A dimensional unit scaling type
   logical, optional,       intent(in)    :: activate !< Specify whether to record the value of
                                                      !! Time as the beginning of the ramp period
 
   ! Local variables
-  real :: deltaTime, wghtA
+  real :: deltaTime ! The time since CS%rampStartTime [T ~> s], which may be negative.
+  real :: wghtA     ! The relative weight of the final value [nondim]
   character(len=12) :: msg
 
   if (CS%truncRampTime==0.) return ! This indicates to ramping is turned off
@@ -1958,7 +1960,7 @@ subroutine updateCFLtruncationValue(Time, CS, activate)
     endif
   endif
   if (.not.CS%CFLrampingIsActivated) return
-  deltaTime = max( 0., time_type_to_real( Time - CS%rampStartTime ) )
+  deltaTime = max( 0., US%s_to_T*time_type_to_real( Time - CS%rampStartTime ) )
   if (deltaTime >= CS%truncRampTime) then
     CS%CFL_trunc = CS%CFL_truncE
     CS%truncRampTime = 0. ! This turns off ramping after this call
@@ -1966,7 +1968,7 @@ subroutine updateCFLtruncationValue(Time, CS, activate)
     wghtA = min( 1., deltaTime / CS%truncRampTime ) ! Linear profile in time
     !wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time
     !wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile
-    wghtA = 1. - ( (1. - wghtA)**2 ) ! Convert linear profiel to nverted parabolic profile
+    wghtA = 1. - ( (1. - wghtA)**2 ) ! Convert linear profile to inverted parabolic profile
     CS%CFL_trunc = CS%CFL_truncS + wghtA * ( CS%CFL_truncE - CS%CFL_truncS )
   endif
   write(msg(1:12),'(es12.3)') CS%CFL_trunc

From dad675ad84c966adf942c0cbb40993d589e347a6 Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <adcroft@users.noreply.github.com>
Date: Wed, 22 Dec 2021 16:01:25 -0500
Subject: [PATCH 7/7] Fix badge URL for codecov

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index d041a47daf..46774baaf0 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
 [![Read The Docs Status](https://readthedocs.org/projects/mom6/badge/?badge=latest)](http://mom6.readthedocs.io/)
-[![codecov](https://codecov.io/gh/NOAA-GFDL/MOM6/branch/dev%2Fmaster/graph/badge.svg)](https://codecov.io/gh/NOAA-GFDL/MOM6)
+[![codecov](https://codecov.io/gh/NOAA-GFDL/MOM6/branch/dev/gfdl/graph/badge.svg?token=uF8SVydCdp)](https://codecov.io/gh/NOAA-GFDL/MOM6)
 
 # MOM6