+Rescaled the units of fluxes%heat_content_... vars

  +Rescaled the units of the 9 fluxes%heat_content_... variables to units of
[J kg-1 R Z T-1], and of tv%TempxPmE to units of [degC R Z] for greater
dimensional consistency testing and for code simplification.  All answers are
bitwise identical.

config_src/coupled_driver/MOM_surface_forcing_gfdl.F90

@@ -464,13 +464,13 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc
     endif
 
     if (associated(IOB%runoff_hflx)) then
-      fluxes%heat_content_lrunoff(i,j) = IOB%runoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
+      fluxes%heat_content_lrunoff(i,j) = kg_m2_s_conversion * IOB%runoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
       if (CS%check_no_land_fluxes) &
         call check_mask_val_consistency(IOB%runoff_hflx(i-i0,j-j0), G%mask2dT(i,j), i, j, 'runoff_hflx', G)
     endif
 
     if (associated(IOB%calving_hflx)) then
-      fluxes%heat_content_frunoff(i,j) = IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
+      fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion * IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
       if (CS%check_no_land_fluxes) &
         call check_mask_val_consistency(IOB%calving_hflx(i-i0,j-j0), G%mask2dT(i,j), i, j, 'calving_hflx', G)
     endif

config_src/ice_solo_driver/MOM_surface_forcing.F90

@@ -731,7 +731,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, CS)
       fluxes%sw(i,j)                   = fluxes%sw(i,j)               * G%mask2dT(i,j)
       fluxes%latent(i,j)               = fluxes%latent(i,j)           * G%mask2dT(i,j)
 
-      fluxes%heat_content_lrunoff(i,j) = fluxes%C_p*US%R_to_kg_m3*US%Z_to_m*US%s_to_T * &
+      fluxes%heat_content_lrunoff(i,j) = fluxes%C_p * &
                                          fluxes%lrunoff(i,j)*sfc_state%SST(i,j)
       fluxes%latent_evap_diag(i,j)     = fluxes%latent_evap_diag(i,j) * G%mask2dT(i,j)
       fluxes%latent_fprec_diag(i,j)    = -US%R_to_kg_m3*US%Z_to_m*US%s_to_T*fluxes%fprec(i,j)*hlf

config_src/mct_driver/mom_surface_forcing_mct.F90

@@ -458,7 +458,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, &
       fluxes%heat_content_frunoff(i,j) = 0.0 * G%mask2dT(i,j)
 
     if (associated(IOB%calving_hflx)) &
-         fluxes%heat_content_frunoff(i,j) = IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
+         fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion * IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
 
     ! longwave radiation, sum up and down (W/m2)
     if (associated(IOB%lw_flux)) &

config_src/nuopc_driver/mom_surface_forcing_nuopc.F90

@@ -458,7 +458,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, &
       fluxes%heat_content_lrunoff(i,j) = IOB%runoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
 
     if (associated(IOB%calving_hflx)) &
-       fluxes%heat_content_frunoff(i,j) = IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
+       fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion*IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
 
     if (associated(IOB%lw_flux)) &
          fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j)

src/core/MOM.F90

@@ -2926,7 +2926,7 @@ subroutine extract_surface_state(CS, sfc_state)
   if (allocated(sfc_state%TempxPmE) .and. associated(CS%tv%TempxPmE)) then
     !$OMP parallel do default(shared)
     do j=js,je ; do i=is,ie
-      sfc_state%TempxPmE(i,j) = CS%tv%TempxPmE(i,j)
+      sfc_state%TempxPmE(i,j) = US%R_to_kg_m3*US%Z_to_m*CS%tv%TempxPmE(i,j)
     enddo ; enddo
   endif
   if (allocated(sfc_state%internal_heat) .and. associated(CS%tv%internal_heat)) then

src/core/MOM_variables.F90

@@ -105,7 +105,7 @@ module MOM_variables
   real, dimension(:,:), pointer :: TempxPmE => NULL()
                          !<   The net inflow of water into the ocean times the
                          !! temperature at which this inflow occurs since the
-                         !! last call to calculate_surface_state [degC kg m-2].
+                         !! last call to calculate_surface_state [degC R Z ~> degC kg m-2].
                          !! This should be prescribed in the forcing fields, but
                          !! as it often is not, this is a useful heat budget diagnostic.
   real, dimension(:,:), pointer :: internal_heat => NULL()
@@ -467,7 +467,7 @@ subroutine MOM_thermovar_chksum(mesg, tv, G)
   if (associated(tv%salt_deficit)) &
     call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", G%HI)
   if (associated(tv%TempxPmE)) &
-    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI)
+    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI, scale=G%US%R_to_kg_m3*G%US%Z_to_m)
 end subroutine MOM_thermovar_chksum
 
 end module MOM_variables

src/diagnostics/MOM_diagnostics.F90

@@ -1796,7 +1796,7 @@ subroutine register_surface_diags(Time, G, IDs, diag, tv)
   IDs%id_salt_deficit = register_diag_field('ocean_model', 'salt_deficit', diag%axesT1, Time, &
          'Salt sink in ocean due to ice flux', 'psu m-2 s-1')
   IDs%id_Heat_PmE = register_diag_field('ocean_model', 'Heat_PmE', diag%axesT1, Time, &
-         'Heat flux into ocean from mass flux into ocean', 'W m-2')
+         'Heat flux into ocean from mass flux into ocean', 'W m-2', conversion=G%US%R_to_kg_m3*G%US%Z_to_m)
   IDs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, Time,&
          'Heat flux into ocean from geothermal or other internal sources', 'W m-2')
 

src/diagnostics/MOM_sum_output.F90

@@ -1012,7 +1012,7 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS)
     ! smg: new code
     ! include heat content from water transport across ocean surface
 !    if (associated(fluxes%heat_content_lprec)) then ; do j=js,je ; do i=is,ie
-!      heat_in(i,j) = heat_in(i,j) + dt*US%L_to_m**2*G%areaT(i,j) * &
+!      heat_in(i,j) = heat_in(i,j) + dt_in_T*RZL2_to_kg*G%areaT(i,j) * &
 !         (fluxes%heat_content_lprec(i,j)   + (fluxes%heat_content_fprec(i,j)   &
 !       + (fluxes%heat_content_lrunoff(i,j) + (fluxes%heat_content_frunoff(i,j) &
 !       + (fluxes%heat_content_cond(i,j)    + (fluxes%heat_content_vprec(i,j)   &
@@ -1022,7 +1022,7 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS)
     ! smg: old code
     if (associated(tv%TempxPmE)) then
       do j=js,je ; do i=is,ie
-        heat_in(i,j) = heat_in(i,j) + (C_p * US%L_to_m**2*G%areaT(i,j)) * tv%TempxPmE(i,j)
+        heat_in(i,j) = heat_in(i,j) + (C_p * RZL2_to_kg*G%areaT(i,j)) * tv%TempxPmE(i,j)
       enddo ; enddo
     elseif (associated(fluxes%evap)) then
       do j=js,je ; do i=is,ie

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -1118,11 +1118,11 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot,      &
                  (dRcv_dT(i)*(Net_heat(i) + Pen_absorbed) - &
                   dRcv_dS(i) * (netMassIn(i) * S(i,1) - Net_salt(i)))
     Conv_En(i) = 0.0 ; dKE_FC(i) = 0.0
-    if (associated(fluxes%heat_content_massin))                            &
-           fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + US%s_to_T * &
-                         T_precip * netMassIn(i) * GV%H_to_kg_m2 * fluxes%C_p * Idt
+    if (associated(fluxes%heat_content_massin)) &
+      fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
+                         T_precip * netMassIn(i) * GV%H_to_RZ * fluxes%C_p * Idt
     if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
-                         T_precip * netMassIn(i) * GV%H_to_kg_m2
+                         T_precip * netMassIn(i) * GV%H_to_RZ
   endif ; enddo
 
   ! Now do netMassOut case in this block.
@@ -1168,14 +1168,14 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot,      &
         d_eb(i,k) = d_eb(i,k) - h_evap
 
         ! smg: when resolve the A=B code, we will set
-        ! heat_content_massout = heat_content_massout - T(i,k)*h_evap*GV%H_to_kg_m2*fluxes%C_p*Idt
+        ! heat_content_massout = heat_content_massout - T(i,k)*h_evap*GV%H_to_RZ*fluxes%C_p*Idt
         ! by uncommenting the lines here.
         ! we will also then completely remove TempXpme from the model.
         if (associated(fluxes%heat_content_massout)) &
-           fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) - US%s_to_T * &
-                                      T(i,k)*h_evap*GV%H_to_kg_m2 * fluxes%C_p * Idt
+          fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) - &
+                                      T(i,k)*h_evap*GV%H_to_RZ * fluxes%C_p * Idt
         if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
-                                      T(i,k)*h_evap*GV%H_to_kg_m2
+                                      T(i,k)*h_evap*GV%H_to_RZ
 
       endif
 

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -882,8 +882,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt_in_T, fluxes, optics, nsw,
   ! Local variables
   integer, parameter :: maxGroundings = 5
   integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings)
-  real :: H_limit_fluxes, IforcingDepthScale
-  real :: Idt
+  real :: H_limit_fluxes
+  real :: IforcingDepthScale
+  real :: Idt        ! The inverse of the timestep [T-1 ~> s-1]
   real :: dThickness, dTemp, dSalt
   real :: fractionOfForcing, hOld, Ithickness
   real :: RivermixConst  ! A constant used in implementing river mixing [Pa s].
@@ -942,7 +943,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt_in_T, fluxes, optics, nsw,
   if (.not.associated(fluxes%sw)) return
 
 #define _OLD_ALG_
-  Idt = 1.0/ (US%T_to_s*dt_in_T)
+  Idt = 1.0 / dt_in_T
 
   calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS))
   calculate_buoyancy = present(SkinBuoyFlux)
@@ -1112,12 +1113,12 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt_in_T, fluxes, optics, nsw,
           ! Diagnostics of heat content associated with mass fluxes
           if (associated(fluxes%heat_content_massin))                             &
             fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) +   &
-                         T2d(i,k) * max(0.,dThickness) * GV%H_to_kg_m2 * fluxes%C_p * Idt
+                         T2d(i,k) * max(0.,dThickness) * GV%H_to_RZ * fluxes%C_p * Idt
           if (associated(fluxes%heat_content_massout))                            &
             fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
-                         T2d(i,k) * min(0.,dThickness) * GV%H_to_kg_m2 * fluxes%C_p * Idt
+                         T2d(i,k) * min(0.,dThickness) * GV%H_to_RZ * fluxes%C_p * Idt
           if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
-                         T2d(i,k) * dThickness * GV%H_to_kg_m2
+                         T2d(i,k) * dThickness * GV%H_to_RZ
 
           ! Determine the energetics of river mixing before updating the state.
           if (calculate_energetics .and. associated(fluxes%lrunoff) .and. CS%do_rivermix) then
@@ -1193,14 +1194,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt_in_T, fluxes, optics, nsw,
           dTemp = dTemp + dThickness*T2d(i,k)
 
           ! Diagnostics of heat content associated with mass fluxes
-          if (associated(fluxes%heat_content_massin))                             &
-            fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) +   &
-                         T2d(i,k) * max(0.,dThickness) * GV%H_to_kg_m2 * fluxes%C_p * Idt
-          if (associated(fluxes%heat_content_massout))                            &
+          if (associated(fluxes%heat_content_massin)) &
+            fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
+                         T2d(i,k) * max(0.,dThickness) * GV%H_to_RZ * fluxes%C_p * Idt
+          if (associated(fluxes%heat_content_massout)) &
             fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
-                         T2d(i,k) * min(0.,dThickness) * GV%H_to_kg_m2 * fluxes%C_p * Idt
+                         T2d(i,k) * min(0.,dThickness) * GV%H_to_RZ * fluxes%C_p * Idt
           if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
-                         T2d(i,k) * dThickness * GV%H_to_kg_m2
+                         T2d(i,k) * dThickness * GV%H_to_RZ
 
           ! Update state by the appropriate increment.
           hOld     = h2d(i,k)               ! Keep original thickness in hand
@@ -1304,7 +1305,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt_in_T, fluxes, optics, nsw,
 
       ! convergence of SW into a layer
       do k=1,nz ; do i=is,ie
-        CS%penSW_diag(i,j,k) = (T2d(i,k)-CS%penSW_diag(i,j,k))*h(i,j,k) * Idt * tv%C_p * GV%H_to_kg_m2
+        CS%penSW_diag(i,j,k) = (T2d(i,k)-CS%penSW_diag(i,j,k))*h(i,j,k) * US%s_to_T*Idt * tv%C_p * GV%H_to_kg_m2
       enddo ; enddo
 
       ! Perform a cumulative sum upwards from bottom to