+Refactor totalTandS to work with scaled variables

Hallberg-NOAA · Hallberg-NOAA · commit 15d9b37eea34 · 2025-04-08T10:00:48.000-04:00
Refactored totalTandS and totalStuff to optionally work with scaled
variables, by adding an optional unit_scale_type argument and a optional
argument specifying the unscaling of thickness to totalTandS and an optional
unscale argument to totalStuff.  The comments describing the units of 19
variables were modified to reflect the various units that might be used.  All
solutions are bitwise identical, and output is unchanged when dimensional
rescaling is not being used, but the debugging output can now be unaltered by
the use of dimensional rescaling.  There are new optional arguments to two
publicly visible routines.

  This change has been tested via calls to totalTandS added to step_MOM_thermo,
but because totalTandS is only intended for debugging, these testing calls are
commented out.  I am uncertain whether to ultimately retain these comments to
illustrate the use of totalTandS or whether to delete them before this PR is
merged in, but retaining them for now seems like they may help the PR review
process.
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
@@ -5,7 +5,7 @@ module MOM
 
 ! Infrastructure modules
 use MOM_array_transform,      only : rotate_array, rotate_vector
-use MOM_debugging,            only : MOM_debugging_init, hchksum, uvchksum
+use MOM_debugging,            only : MOM_debugging_init, hchksum, uvchksum, totalTandS
 use MOM_debugging,            only : check_redundant, query_debugging_checks
 use MOM_checksum_packages,    only : MOM_thermo_chksum, MOM_state_chksum
 use MOM_checksum_packages,    only : MOM_accel_chksum, MOM_surface_chksum
@@ -1720,6 +1720,13 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
 
   call disable_averaging(CS%diag)
 
+! This works in general:
+!  if (associated(tv%T)) &
+!    call totalTandS(G%HI, h, G%areaT, tv%T, tv%S, "End of step_MOM", US, GV%H_to_mks)
+! This works only if there is no rescaling being used:
+!  if (associated(tv%T)) &
+!    call totalTandS(G%HI, h, G%areaT, tv%T, tv%S, "End of step_MOM")
+
   if (showCallTree) call callTree_leave("step_MOM_thermo(), MOM.F90")
 
 end subroutine step_MOM_thermo
diff --git a/src/diagnostics/MOM_debugging.F90 b/src/diagnostics/MOM_debugging.F90
@@ -8,16 +8,18 @@ module MOM_debugging
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
-use MOM_checksums, only : hchksum, Bchksum, qchksum, uvchksum, hchksum_pair
-use MOM_checksums, only : is_NaN, chksum, MOM_checksums_init
-use MOM_coms, only : PE_here, root_PE, num_PEs
-use MOM_coms, only : min_across_PEs, max_across_PEs, reproducing_sum
-use MOM_domains, only : pass_vector, pass_var, pe_here
-use MOM_domains, only : BGRID_NE, AGRID, To_All, Scalar_Pair
+use MOM_checksums,     only : hchksum, Bchksum, qchksum, uvchksum, hchksum_pair
+use MOM_checksums,     only : is_NaN, chksum, MOM_checksums_init
+use MOM_coms,          only : PE_here, root_PE, num_PEs
+use MOM_coms,          only : min_across_PEs, max_across_PEs, reproducing_sum
+use MOM_domains,       only : pass_vector, pass_var, pe_here
+use MOM_domains,       only : BGRID_NE, AGRID, To_All, Scalar_Pair
 use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
-use MOM_file_parser, only : log_version, param_file_type, get_param
-use MOM_grid, only : ocean_grid_type
-use MOM_hor_index, only : hor_index_type
+use MOM_file_parser,   only : log_version, param_file_type, get_param
+use MOM_grid,          only : ocean_grid_type
+use MOM_hor_index,     only : hor_index_type
+use MOM_io,            only : stdout
+use MOM_unit_scaling,  only : unit_scale_type
 
 implicit none ; private
 
@@ -837,67 +839,101 @@ end subroutine chksum_vec_A2d
 
 !> This function returns the sum over computational domain of all
 !! processors of hThick*stuff, where stuff is a 3-d array at tracer points.
-function totalStuff(HI, hThick, areaT, stuff)
+function totalStuff(HI, hThick, areaT, stuff, unscale)
   type(hor_index_type),               intent(in) :: HI     !< A horizontal index type
-  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights [m]
-  real, dimension(HI%isd:,HI%jsd:),   intent(in) :: areaT  !< The array of cell areas [m2]
-  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: stuff  !< The array of stuff to be summed in arbitrary units [a]
-  real                                         :: totalStuff !< the globally integrated amount of stuff [a m3]
+  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights
+                                                           !! [H ~> m or kg m-2] or [m] or [kg m-2]
+  real, dimension(HI%isd:,HI%jsd:),   intent(in) :: areaT  !< The array of cell areas [L2 ~> m2] or [m2]
+  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: stuff  !< The array of stuff to be summed in arbitrary
+                                                           !! units [A ~> a] or [a]
+  real,                     optional, intent(in) :: unscale !< A factor that is used to undo scaling of the array
+                                                           !! and the cell mass or volume before it is summed in
+                                                           !! [a m3 A-1 H-1 L-2 ~> 1] or [a kg A-1 H-1 L-2 ~> 1]
+  real                                           :: totalStuff !< the globally integrated amount of stuff
+                                                           !! [A H L2 ~> a m3 or a kg] or [a m3]
   ! Local variables
-  real, dimension(HI%isc:HI%iec, HI%jsc:HI%jec) :: tmp_for_sum ! The column integrated amount of stuff in a cell [a m3]
+  real  :: tmp_for_sum(HI%isc:HI%iec, HI%jsc:HI%jec)  ! The column integrated amount of stuff in a
+                                                      ! cell [A H L2 ~> a m3 or a kg] or [a m3]
   integer :: i, j, k, nz
 
   nz = size(hThick,3)
   tmp_for_sum(:,:) = 0.0
   do k=1,nz ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec
     tmp_for_sum(i,j) = tmp_for_sum(i,j) + hThick(i,j,k) * stuff(i,j,k) * areaT(i,j)
   enddo ; enddo ; enddo
-  totalStuff = reproducing_sum(tmp_for_sum)
+  totalStuff = reproducing_sum(tmp_for_sum, unscale=unscale)
 
 end function totalStuff
 
 !> This subroutine display the total thickness, temperature and salinity
 !! as well as the change since the last call.
-subroutine totalTandS(HI, hThick, areaT, temperature, salinity, mesg)
+subroutine totalTandS(HI, hThick, areaT, temperature, salinity, mesg, US, H_to_mks)
   type(hor_index_type),               intent(in) :: HI     !< A horizontal index type
-  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights [m]
-  real, dimension(HI%isd:,HI%jsd:),   intent(in) :: areaT  !< The array of cell areas [m2]
-  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: temperature !< The temperature field to sum [degC]
-  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: salinity    !< The salinity field to sum [ppt]
+  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights
+                                                           !! [H ~> m or kg m-2] or [m] or [kg m-2]
+  real, dimension(HI%isd:,HI%jsd:),   intent(in) :: areaT  !< The array of cell areas [L2 ~> m2] or [m2]
+  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: temperature !< The temperature field to sum [C ~> degC] or [degC]
+  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: salinity    !< The salinity field to sum [S ~> ppt] or [ppt]
   character(len=*),                   intent(in) :: mesg        !< An identifying message
+  type(unit_scale_type),    optional, intent(in) :: US       !< A dimensional unit scaling type
+  real,                     optional, intent(in) :: H_to_MKS !< A constant that translates thickness units to its
+                                                             !! MKS units (m or kg m-2) based on whether the model is
+                                                             !! Boussinesq [m H-1 ~> 1] or not [kg m-2 H-1 ~> 1]
   ! NOTE: This subroutine uses "save" data which is not thread safe and is purely for
   ! extreme debugging without a proper debugger.
-  real, save :: totalH = 0.   ! The total ocean volume, saved for the next call [m3]
-  real, save :: totalT = 0.   ! The total volume integrated ocean temperature, saved for the next call [degC m3]
-  real, save :: totalS = 0.   ! The total volume integrated ocean salinity, saved for the next call [ppt m3]
+  real, save :: totalH = 0.   ! The total ocean volume or mass, saved for the next
+                              ! call [H L2 ~> m3 or kg] or [m3] or [kg]
+  real, save :: totalT = 0.   ! The total volume integrated ocean temperature, saved for the next
+                              ! call [C H L2 ~> degC m3 or degC kg] or [degC m3] or [degC kg]
+  real, save :: totalS = 0.   ! The total volume integrated ocean salinity, saved for the next
+                              ! call [S H L2 ~> ppt m3 or ppt kg] or [ppt m3] or [ppt kg]
   ! Local variables
   logical, save :: firstCall = .true.
-  real, dimension(HI%isc:HI%iec, HI%jsc:HI%jec) :: tmp_for_sum ! The volume of each column [m3]
-  real :: thisH, delH  ! The total ocean volume and the change from the last call [m3]
-  real :: thisT, delT  ! The current total volume integrated temperature and the change from the last call [degC m3]
-  real :: thisS, delS  ! The current total volume integrated salinity and the change from the last call [ppt m3]
+  real :: tmp_for_sum(HI%isc:HI%iec, HI%jsc:HI%jec) ! The volume of each column [H L2 ~> m3 or kg] or [m3] or [kg]
+  real :: thisH, delH  ! The total ocean volume and the change from the last call [H L2 ~> m3 or kg] or [m3] or [kg]
+  real :: thisT, delT  ! The current total volume integrated temperature and the change from the last
+                       ! call [C H L2 ~> degC m3 or degC kg] or [degC m3] or [degC kg]
+  real :: thisS, delS  ! The current total volume integrated salinity and the change from the last
+                       ! call [S H L2 ~> ppt m3 or ppt kg] or [ppt m3] or [ppt kg]
+  real :: H_unscale    ! A constant that translates thickness units to its MKS units (m or kg m-2) based on
+                       ! whether the model is Boussinesq [m H-1 ~> 1] or non-Boussinesq [kg m-2 H-1 ~> 1]
+  real :: HL2_unscale  ! An overall unscaling factor for cell mass or volume [m3 H-1 L-2 ~> 1] or [kg H-1 L-2 ~> 1]
+  real :: T_unscale    ! An overall unscaling factor for cell-integrated temperature [degC m3 C-1 H-1 L-2 ~> 1] or
+                       ! [degC kg C-1 H-1 L-2 ~> 1]
+  real :: S_unscale    ! An overall unscaling factor for cell-integrated salinity [ppt m3 S-1 H-1 L-2 ~> 1] or
+                       ! [ppt kg S-1 H-1 L-2 ~> 1]
   integer :: i, j, k, nz
 
+  H_unscale = 1.0 ; if (present(H_to_mks)) H_unscale = H_to_mks
+  if (present(US)) then
+    HL2_unscale = US%L_to_m**2 * H_unscale
+    T_unscale = US%C_to_degC * HL2_unscale ; S_unscale = US%S_to_ppt * HL2_unscale
+  else
+    HL2_unscale = H_unscale
+    T_unscale = HL2_unscale ; S_unscale = HL2_unscale
+  endif
+
   nz = size(hThick,3)
   tmp_for_sum(:,:) = 0.0
   do k=1,nz ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec
     tmp_for_sum(i,j) = tmp_for_sum(i,j) + hThick(i,j,k) * areaT(i,j)
   enddo ; enddo ; enddo
-  thisH = reproducing_sum(tmp_for_sum)
-  thisT = totalStuff(HI, hThick, areaT, temperature)
-  thisS = totalStuff(HI, hThick, areaT, salinity)
+  thisH = reproducing_sum(tmp_for_sum, unscale=HL2_unscale)
+  thisT = totalStuff(HI, hThick, areaT, temperature, unscale=T_unscale)
+  thisS = totalStuff(HI, hThick, areaT, salinity, unscale=S_unscale)
 
   if (is_root_pe()) then
     if (firstCall) then
       totalH = thisH ; totalT = thisT ; totalS = thisS
-      write(0,*) 'Totals H,T,S:',thisH,thisT,thisS,' ',mesg
+      write(stdout,*) 'Totals H,T,S:', thisH*HL2_unscale, thisT*T_unscale, thisS*S_unscale, ' ', mesg
       firstCall = .false.
     else
       delH = thisH - totalH
       delT = thisT - totalT
       delS = thisS - totalS
       totalH = thisH ; totalT = thisT ; totalS = thisS
-      write(0,*) 'Tot/del H,T,S:',thisH,thisT,thisS,delH,delT,delS,' ',mesg
+      write(0,*) 'Tot/del H,T,S:', thisH*HL2_unscale, thisT*T_unscale, thisS*S_unscale, &
+                                   delH*HL2_unscale,  delT*T_unscale,  delS*S_unscale, ' ', mesg
     endif
   endif
 
