+Thickness-based diffusivity arguments

  Rescaled diapycnal diffusivities passed as arguments in non-Boussinesq mode,
to be equivalent to the thermal conductivity divided by the heat capacity,
analogously to the difference between a kinematic viscosity and a dynamic
viscosity, so that the new units are [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

  This includes changing the units of 4 arguments to set_diffusivity;
3 arguments each to calculate_bkgnd_mixing, add_drag_diffusivity,
add_LOTW_BBL_diffusivity, user_change_diff, calculate_tidal_mixing and
add_int_tide_diffusivity; 2 arguments to KPP_calculate, calculate_CVMix_conv,
compute_ddiff_coeffs, differential_diffuse_T_S, entrainment_diffusive,
double_diffusion, add_MLrad_diffusivity, and calculate_CVMix_tidal; and one
argument to energetic_PBL.

  The units of 36 internal variables were also changed, as were a total of
29 elements in various opaque types, including 8 elements in bkgnd_mixing_cs,
2 in diabatic_CC, 3 in tidal_mixing_diags type, 1 in user_change_diff_CS,
9 in set_diffusivity_CS type, and 6 elements in diffusivity_diags.

  Two new internal variables were added, and several redundant GV%H_to_Z
conversion factors were also cancelled out, some using that
GV%H_to_Z*GV%Rho0 = GV%H_to_RZ.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.

src/parameterizations/vertical/MOM_CVMix_KPP.F90

@@ -536,7 +536,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
       'Heat diffusivity due to KPP, as calculated by [CVMix] KPP', &
       'm2/s', conversion=US%Z2_T_to_m2_s)
   CS%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, Time, &
-      'Diffusivity passed to KPP', 'm2/s', conversion=US%Z2_T_to_m2_s)
+      'Diffusivity passed to KPP', 'm2/s', conversion=GV%HZ_T_to_m2_s)
   CS%id_Ks_KPP = register_diag_field('ocean_model', 'KPP_Ksalt', diag%axesTi, Time, &
       'Salt diffusivity due to KPP, as calculated by [CVMix] KPP', &
       'm2/s', conversion=US%Z2_T_to_m2_s)
@@ -610,10 +610,10 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kt  !< (in)  Vertical diffusivity of heat w/o KPP
                                                                     !! (out) Vertical diffusivity including KPP
-                                                                    !!       [Z2 T-1 ~> m2 s-1]
+                                                                    !!       [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Ks  !< (in)  Vertical diffusivity of salt w/o KPP
                                                                     !! (out) Vertical diffusivity including KPP
-                                                                    !!       [Z2 T-1 ~> m2 s-1]
+                                                                    !!       [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kv  !< (in)  Vertical viscosity w/o KPP
                                                                     !! (out) Vertical viscosity including KPP
                                                                     !!       [H Z T-1 ~> m2 s-1 or Pa s]
@@ -650,8 +650,8 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
     call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m)
     call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0, scale=US%Z_to_m*US%s_to_T)
     call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0, scale=US%L_to_m**2*US%s_to_T**3)
-    call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s)
-    call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s)
+    call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0, scale=GV%HZ_T_to_m2_s)
   endif
 
   nonLocalTrans(:,:) = 0.0
@@ -711,8 +711,8 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
          Kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt [m2 s-1]
          Kviscosity(:)     = 0. ! Viscosity [m2 s-1]
       else
-         Kdiffusivity(:,1) = US%Z2_T_to_m2_s * Kt(i,j,:)
-         Kdiffusivity(:,2) = US%Z2_T_to_m2_s * Ks(i,j,:)
+         Kdiffusivity(:,1) = GV%HZ_T_to_m2_s * Kt(i,j,:)
+         Kdiffusivity(:,2) = GV%HZ_T_to_m2_s * Ks(i,j,:)
          Kviscosity(:) = GV%HZ_T_to_m2_s * Kv(i,j,:)
       endif
 
@@ -860,15 +860,15 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
       if (.not. CS%passiveMode) then
         if (CS%KPPisAdditive) then
           do k=1, GV%ke+1
-            Kt(i,j,k) = Kt(i,j,k) + US%m2_s_to_Z2_T * Kdiffusivity(k,1)
-            Ks(i,j,k) = Ks(i,j,k) + US%m2_s_to_Z2_T * Kdiffusivity(k,2)
+            Kt(i,j,k) = Kt(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,1)
+            Ks(i,j,k) = Ks(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,2)
             Kv(i,j,k) = Kv(i,j,k) + GV%m2_s_to_HZ_T * Kviscosity(k)
             if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = GV%H_to_Z*Kv(i,j,k)
           enddo
         else ! KPP replaces prior diffusivity when former is non-zero
           do k=1, GV%ke+1
-            if (Kdiffusivity(k,1) /= 0.) Kt(i,j,k) = US%m2_s_to_Z2_T * Kdiffusivity(k,1)
-            if (Kdiffusivity(k,2) /= 0.) Ks(i,j,k) = US%m2_s_to_Z2_T * Kdiffusivity(k,2)
+            if (Kdiffusivity(k,1) /= 0.) Kt(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,1)
+            if (Kdiffusivity(k,2) /= 0.) Ks(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,2)
             if (Kviscosity(k) /= 0.) Kv(i,j,k) = GV%m2_s_to_HZ_T * Kviscosity(k)
             if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = GV%H_to_Z*Kv(i,j,k)
           enddo
@@ -883,8 +883,8 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
   call cpu_clock_end(id_clock_KPP_calc)
 
   if (CS%debug) then
-    call hchksum(Kt, "KPP out: Kt", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
-    call hchksum(Ks, "KPP out: Ks", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
+    call hchksum(Kt, "KPP out: Kt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(Ks, "KPP out: Ks", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
   endif
 
   ! send diagnostics to post_data

src/parameterizations/vertical/MOM_CVMix_conv.F90

@@ -143,15 +143,16 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
   type(CVMix_conv_cs),                       intent(in)  :: CS !< CVMix convection control structure
   real, dimension(SZI_(G),SZJ_(G)),          intent(in)  :: hbl !< Depth of ocean boundary layer [Z ~> m]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
-                                             intent(inout) :: Kd !< Diapycnal diffusivity at each interface that
-                                                                 !! will be incremented here [Z2 T-1 ~> m2 s-1].
+                                             intent(inout) :: Kd !< Diapycnal diffusivity at each interface
+                                                                 !! that will be incremented here
+                                                                 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
-                                             intent(inout) :: KV !< Viscosity at each interface that will be
+                                             intent(inout) :: Kv !< Viscosity at each interface that will be
                                                                  !! incremented here [H Z T-1 ~> m2 s-1 or Pa s]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                                    optional, intent(inout) :: Kd_aux !< A second diapycnal diffusivity at each
                                                                  !! interface that will also be incremented
-                                                                 !! here [Z2 T-1 ~> m2 s-1].
+                                                                 !! here [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
 
   ! local variables
   real, dimension(SZK_(GV)) :: rho_lwr !< Adiabatic Water Density [kg m-3], this is a dummy
@@ -238,12 +239,12 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
 
       ! Increment the diffusivity outside of the boundary layer.
       do K=max(1,kOBL+1),GV%ke+1
-        Kd(i,j,K) = Kd(i,j,K) + US%m2_s_to_Z2_T * kd_col(K)
+        Kd(i,j,K) = Kd(i,j,K) + GV%m2_s_to_HZ_T * kd_col(K)
       enddo
       if (present(Kd_aux)) then
         ! Increment the other diffusivity outside of the boundary layer.
         do K=max(1,kOBL+1),GV%ke+1
-          Kd_aux(i,j,K) = Kd_aux(i,j,K) + US%m2_s_to_Z2_T * kd_col(K)
+          Kd_aux(i,j,K) = Kd_aux(i,j,K) + GV%m2_s_to_HZ_T * kd_col(K)
         enddo
       endif
 
@@ -277,7 +278,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
     !   call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
     ! if (CS%id_kv_conv > 0) &
     !   call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
-    call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
+    call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
     call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
   endif
 

src/parameterizations/vertical/MOM_CVMix_ddiff.F90

@@ -150,9 +150,11 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho)
   integer,                                    intent(in)    :: j    !< Meridional grid index to work on.
   ! Kd_T and Kd_S are intent inout because only one j-row is set here, but they are essentially outputs.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kd_T !< Interface double diffusion diapycnal
-                                                                    !! diffusivity for temp [Z2 T-1 ~> m2 s-1].
+                                                                    !! diffusivity for temperature
+                                                                    !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kd_S !< Interface double diffusion diapycnal
-                                                                    !! diffusivity for salt [Z2 T-1 ~> m2 s-1].
+                                                                    !! diffusivity for salinity
+                                                                    !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
   type(CVMix_ddiff_cs),                       pointer       :: CS   !< The control structure returned
                                                                     !! by a previous call to CVMix_ddiff_init.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
@@ -254,8 +256,8 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho)
                             nlev=GV%ke,    &
                             max_nlev=GV%ke)
     do K=1,GV%ke+1
-      Kd_T(i,j,K) = US%m2_s_to_Z2_T * Kd1_T(K)
-      Kd_S(i,j,K) = US%m2_s_to_Z2_T * Kd1_S(K)
+      Kd_T(i,j,K) = GV%m2_s_to_HZ_T * Kd1_T(K)
+      Kd_S(i,j,K) = GV%m2_s_to_HZ_T * Kd1_S(K)
     enddo
 
     ! Do not apply mixing due to convection within the boundary layer