Skip to content

Commit

Permalink
+Refactor totalTandS to work with scaled variables
Browse files Browse the repository at this point in the history 
  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.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Feb 5, 2025
1 parent bb66d8a commit 508d45f
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 34 deletions.
9 changes: 8 additions & 1 deletion src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -1840,6 +1840,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
Expand Down
101 changes: 68 additions & 33 deletions src/diagnostics/MOM_debugging.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,17 @@ 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_unit_scaling, only : unit_scale_type

implicit none ; private

Expand Down Expand Up @@ -837,67 +838,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 not [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(0,*) '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

Expand Down

0 comments on commit 508d45f

Please sign in to comment.