+Rescaled specific volume units of dSV_dT args

  Rescaled the specific volume (density) units of the dSV_dT and dSV_dS
variables passed to energetic_PBL, applyBoundaryFluxesInOut, and
absorbRemainingSW for dimensional consistency testing.  Also rescaled the
dimensions of TKE returned from absorbRemainingSW.  All answers are bitwise
identical, but the units of a arguments to 3 public interfaces have changed.

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -871,10 +871,10 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
                                                !! forcing through each layer [R Z3 T-2 ~> J m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                  optional, intent(out)   :: dSV_dT !< Partial derivative of specific volume with
-                                               !! potential temperature [m3 kg-1 degC-1].
+                                               !! potential temperature [R-1 degC-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                  optional, intent(out)   :: dSV_dS !< Partial derivative of specific volume with
-                                               !! salinity [m3 kg-1 ppt-1].
+                                               !! salinity [R-1 ppt-1].
   real, dimension(SZI_(G),SZJ_(G)), &
                    optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].
 
@@ -909,8 +909,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
     h2d, &           ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
     T2d, &           ! A 2-d copy of the layer temperatures [degC]
     pen_TKE_2d, &    ! The TKE required to homogenize the heating by shortwave radiation within
-                     ! a layer [kg m-3 Z3 T-2 ~> J m-2]
-    dSV_dT_2d        ! The partial derivative of specific volume with temperature [m3 kg-1 degC-1]
+                     ! a layer [R Z3 T-2 ~> J m-2]
+    dSV_dT_2d        ! The partial derivative of specific volume with temperature [R-1 degC-1]
   real, dimension(SZI_(G),SZK_(G)+1) :: netPen
   real, dimension(max(nsw,1),SZI_(G)) :: &
     Pen_SW_bnd, &    ! The penetrative shortwave heating integrated over a timestep by band
@@ -922,9 +922,10 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
                      ! band of shortwave radation in each layer [H-1 ~> m-1 or m2 kg-1]
   real, dimension(maxGroundings) :: hGrounding
   real    :: Temp_in, Salin_in
-!  real   :: I_G_Earth ! The inverse of the gravitational acceleration with conversion factors [s2 m-1].
+! real   :: I_G_Earth ! The inverse of the gravitational acceleration with conversion factors [R m2 kg-1 s2 ~> s2 m-1]
   real    :: dt_in_T ! The time step converted to T units [T ~> s]
-  real    :: g_Hconv2
+  real    :: g_Hconv2 ! A conversion factor for use in the TKE calculation
+                      ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
   real    :: GoRho    ! g_Earth times a unit conversion factor divided by density
                       ! [Z3 m T-2 kg-1 ~> m4 s-2 kg-1]
   logical :: calculate_energetics
@@ -945,8 +946,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
   calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS))
   calculate_buoyancy = present(SkinBuoyFlux)
   if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0
-!  I_G_Earth = US%Z_to_m / (US%L_T_to_m_s**2 * GV%g_Earth)
-  g_Hconv2 = (US%m_to_Z**3 * US%T_to_s**2) * GV%H_to_Pa * GV%H_to_kg_m2*US%kg_m3_to_R
+!  I_G_Earth = US%kg_m3_to_R*US%Z_to_m / (US%L_T_to_m_s**2 * GV%g_Earth)
+  g_Hconv2 = (US%m_to_Z**3 * US%T_to_s**2) * GV%H_to_Pa * GV%H_to_kg_m2*US%kg_m3_to_R**2
 
   if (present(cTKE)) cTKE(:,:,:) = 0.0
   if (calculate_buoyancy) then
@@ -1004,7 +1005,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
           pres(i) = pres(i) + d_pres(i)
         enddo
         call calculate_specific_vol_derivs(T2d(:,k), tv%S(:,j,k), p_lay(:),&
-                 dSV_dT(:,j,k), dSV_dS(:,j,k), is, ie-is+1, tv%eqn_of_state)
+                 dSV_dT(:,j,k), dSV_dS(:,j,k), is, ie-is+1, tv%eqn_of_state, scale=US%R_to_kg_m3)
         do i=is,ie ; dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; enddo
 !        do i=is,ie
 !          dT_to_dPE(i,k) = I_G_Earth * d_pres(i) * p_lay(i) * dSV_dT(i,j,k)
@@ -1134,10 +1135,11 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
             ! rivermix_depth =  The prescribed depth over which to mix river inflow
             ! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
             ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
-            RivermixConst = -0.5*(CS%rivermix_depth*dt)*(US%m_to_Z**3 * US%T_to_s**2) * GV%Z_to_H*GV%H_to_Pa
+            RivermixConst = -0.5*(CS%rivermix_depth*dt)*(US%m_to_Z**3 * US%T_to_s**2) * &
+                            GV%Z_to_H*GV%H_to_Pa*US%kg_m3_to_R
 
-            cTKE(i,j,k) = cTKE(i,j,k) + max(0.0, RivermixConst*US%kg_m3_to_R*dSV_dS(i,j,1) * &
-                  (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1))
+            cTKE(i,j,k) = cTKE(i,j,k) + max(0.0, RivermixConst*dSV_dS(i,j,1) * &
+                  US%kg_m3_to_R*(fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1))
           endif
 
           ! Update state
@@ -1283,7 +1285,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
                              .false., .true., T2d, Pen_SW_bnd, TKE=pen_TKE_2d, dSV_dT=dSV_dT_2d)
       k = 1 ! For setting break-points.
       do k=1,nz ; do i=is,ie
-        cTKE(i,j,k) = cTKE(i,j,k) + US%kg_m3_to_R*pen_TKE_2d(i,k)
+        cTKE(i,j,k) = cTKE(i,j,k) + pen_TKE_2d(i,k)
       enddo ; enddo
     else
       call absorbRemainingSW(G, GV, US, h2d, opacityBand, nsw, optics, j, dt_in_T, H_limit_fluxes, &

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -483,8 +483,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
     hold,   &    ! layer thickness before diapycnal entrainment, and later
                  ! the initial layer thicknesses (if a mixed layer is used),
                  ! [H ~> m or kg m-2]
-    dSV_dT, &    ! The partial derivative of specific volume with temperature [m3 kg-1 degC-1]
-    dSV_dS, &    ! The partial derivative of specific volume with salinity [m3 kg-1 ppt-1].
+    dSV_dT, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
+    dSV_dS, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
     cTKE,   &    ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
     u_h,    &    ! zonal and meridional velocities at thickness points after
     v_h          ! entrainment [L T-1 ~> m s-1]
@@ -835,10 +835,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
       call hchksum(eb_t, "after applyBoundaryFluxes eb_t",G%HI,haloshift=0, scale=GV%H_to_m)
       call hchksum(ea_s, "after applyBoundaryFluxes ea_s",G%HI,haloshift=0, scale=GV%H_to_m)
       call hchksum(eb_s, "after applyBoundaryFluxes eb_s",G%HI,haloshift=0, scale=GV%H_to_m)
-      call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0, &
+      call hchksum(cTKE, "after applyBoundaryFluxes cTKE", G%HI, haloshift=0, &
                    scale=US%R_to_kg_m3*US%Z_to_m**3*US%s_to_T**2)
-      call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0)
-      call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0)
+      call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT", G%HI, haloshift=0, scale=US%kg_m3_to_R)
+      call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS", G%HI, haloshift=0, scale=US%kg_m3_to_R)
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
@@ -1267,8 +1267,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
 !    hold,   &    ! layer thickness before diapycnal entrainment, and later
                  ! the initial layer thicknesses (if a mixed layer is used),
                  ! [H ~> m or kg m-2]
-    dSV_dT, &    ! The partial derivative of specific volume with temperature [m3 kg-1 degC-1]
-    dSV_dS, &    ! The partial derivative of specific volume with salinity [m3 kg-1 ppt-1].
+    dSV_dT, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
+    dSV_dS, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
     cTKE,   &    ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
     u_h,    &    ! zonal and meridional velocities at thickness points after
     v_h          ! entrainment [L T-1 ~> m s-1]
@@ -1568,8 +1568,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
       call hchksum(eb_s, "after applyBoundaryFluxes eb_s",G%HI,haloshift=0, scale=GV%H_to_m)
       call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0, &
                    scale=US%R_to_kg_m3*US%Z_to_m**3*US%s_to_T**2)
-      call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0)
-      call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0)
+      call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0, scale=US%kg_m3_to_R)
+      call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0, scale=US%kg_m3_to_R)
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
@@ -1942,8 +1942,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
     hold,   &    ! layer thickness before diapycnal entrainment, and later
                  ! the initial layer thicknesses (if a mixed layer is used),
                  ! [H ~> m or kg m-2]
-    dSV_dT, &    ! The partial derivative of specific volume with temperature [m3 kg-1 degC-1]
-    dSV_dS, &    ! The partial derivative of specific volume with salinity [m3 kg-1 ppt-1].
+    dSV_dT, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
+    dSV_dS, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
     cTKE,   &    ! convective TKE requirements for each layer [kg m-3 Z3 T-2 ~> J m-2].
     u_h,    &    ! zonal and meridional velocities at thickness points after
     v_h          ! entrainment [L T-1 ~> m s-1]

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -254,10 +254,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                            intent(in)    :: dSV_dT !< The partial derivative of in-situ specific
                                                    !! volume with potential temperature
-                                                   !! [m3 kg-1 degC-1].
+                                                   !! [R-1 degC-1 ~> m3 kg-1 degC-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                            intent(in)    :: dSV_dS !< The partial derivative of in-situ specific
-                                                   !! volume with salinity [m3 kg-1 ppt-1].
+                                                   !! volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                            intent(in)    :: TKE_forced !< The forcing requirements to homogenize the
                                                    !! forcing that has been applied to each layer
@@ -407,7 +407,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
       h_2d(i,k) = h_3d(i,j,k) ; u_2d(i,k) = u_3d(i,j,k) ; v_2d(i,k) = v_3d(i,j,k)
       T_2d(i,k) = tv%T(i,j,k) ; S_2d(i,k) = tv%S(i,j,k)
       TKE_forced_2d(i,k) = TKE_forced(i,j,k)
-      dSV_dT_2d(i,k) = US%R_to_kg_m3*dSV_dT(i,j,k) ; dSV_dS_2d(i,k) = US%R_to_kg_m3*dSV_dS(i,j,k)
+      dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; dSV_dS_2d(i,k) = dSV_dS(i,j,k)
     enddo ; enddo
 
     !   Determine the initial mech_TKE and conv_PErel, including the energy required

src/parameterizations/vertical/MOM_opacity.F90

@@ -556,9 +556,9 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l
   real, dimension(SZI_(G)), optional, intent(inout) :: Ttot !< Depth integrated mixed layer
                                                            !! temperature [degC H ~> degC m or degC kg m-2]
   real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: dSV_dT !< The partial derivative of specific
-                                                           !! volume with temperature [m3 kg-1 degC-1].
+                                                           !! volume with temperature [R-1 degC-1].
   real, dimension(SZI_(G),SZK_(GV)), optional, intent(inout) :: TKE !< The TKE sink from mixing the heating
-                                                           !! throughout a layer [kg m-3 Z3 T-2 ~> J m-2].
+                                                           !! throughout a layer [R Z3 T-2 ~> J m-2].
 
   ! Local variables
   real, dimension(SZI_(G),SZK_(GV)) :: &
@@ -599,7 +599,7 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l
   real :: epsilon           ! A small thickness that must remain in each
                             ! layer, and which will not be subject to heating [H ~> m or kg m-2]
   real :: g_Hconv2          ! A conversion factor for use in the TKE calculation
-                            ! in units of [Z3 kg2 m-6 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
+                            ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
   logical :: SW_Remains     ! If true, some column has shortwave radiation that
                             ! was not entirely absorbed.
   logical :: TKE_calc       ! If true, calculate the implications to the
@@ -618,9 +618,9 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l
   TKE_calc = (present(TKE) .and. present(dSV_dT))
 
   if (optics%answers_2018) then
-    g_Hconv2 = (US%m_to_Z**2 * US%L_to_Z**2*GV%g_Earth * GV%H_to_kg_m2) * GV%H_to_kg_m2
+    g_Hconv2 = (US%m_to_Z**2 * US%L_to_Z**2*GV%g_Earth * GV%H_to_kg_m2) * GV%H_to_kg_m2 * US%kg_m3_to_R**2
   else
-    g_Hconv2 = US%m_to_Z**2 * US%L_to_Z**2*GV%g_Earth * GV%H_to_kg_m2**2
+    g_Hconv2 = US%m_to_Z**2 * US%L_to_Z**2*GV%g_Earth * GV%H_to_kg_m2**2 * US%kg_m3_to_R**2
   endif
 
   h_heat(:) = 0.0