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