From 15d9b37eea34bc5d00470144ab308aa58bef099b Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 5 Feb 2025 05:37:24 -0500 Subject: [PATCH] +Refactor totalTandS to work with scaled variables 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. --- src/core/MOM.F90 | 9 ++- src/diagnostics/MOM_debugging.F90 | 102 ++++++++++++++++++++---------- 2 files changed, 77 insertions(+), 34 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 937cda08f7..89221a2b3f 100644 --- 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 index f454ac8d4a..85af39e377 100644 --- 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,14 +839,21 @@ 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) @@ -852,52 +861,79 @@ function totalStuff(HI, hThick, areaT, stuff) 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