From 8658edfaab36e7e4367c022546aa4ec4cddd1794 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 10 Jul 2024 18:36:08 -0400 Subject: [PATCH 1/2] +Move array_global_min_max to MOM_spatial_means Moved the publicly visible routine array_global_min_max to MOM_spatial_means from MOM_generic_tracer, along with the private supporting function ijk_loc, without any changes to the code itself. This was done because this routine now gives useful debugging information that can be useful outside of the generic tracer module. All answers are bitwise identical, but the module use statements for array_global_min_max will need to be updated. --- src/diagnostics/MOM_spatial_means.F90 | 191 ++++++++++++++++++++++++++ src/tracer/MOM_generic_tracer.F90 | 190 +------------------------ 2 files changed, 193 insertions(+), 188 deletions(-) diff --git a/src/diagnostics/MOM_spatial_means.F90 b/src/diagnostics/MOM_spatial_means.F90 index d8a8abfd99..fe33e38a80 100644 --- a/src/diagnostics/MOM_spatial_means.F90 +++ b/src/diagnostics/MOM_spatial_means.F90 @@ -7,9 +7,11 @@ module MOM_spatial_means use MOM_coms, only : EFP_to_real, real_to_EFP, EFP_sum_across_PEs use MOM_coms, only : reproducing_sum, reproducing_sum_EFP, EFP_to_real use MOM_coms, only : query_EFP_overflow_error, reset_EFP_overflow_error +use MOM_coms, only : max_across_PEs, min_across_PEs use MOM_error_handler, only : MOM_error, NOTE, WARNING, FATAL, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type +use MOM_hor_index, only : hor_index_type use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -21,6 +23,7 @@ module MOM_spatial_means public :: global_area_integral public :: global_volume_mean, global_mass_integral, global_mass_int_EFP public :: adjust_area_mean_to_zero +public :: array_global_min_max ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -701,4 +704,192 @@ subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale, unscale) end subroutine adjust_area_mean_to_zero + +!> Find the global maximum and minimum of a tracer array and return the locations of the extrema. +!! When there multiple cells with the same extreme values, the reported locations are from the +!! uppermost layer where they occur, and then from the logically northernmost and then eastermost +!! such location on the unrotated version of the grid within that layer. Only ocean points (as +!! indicated by a positive value of G%mask2dT) are evaluated, and if there are no ocean points +!! anywhere in the domain, the reported extrema and their locations are all returned as 0. +subroutine array_global_min_max(tr_array, G, nk, g_min, g_max, & + xgmin, ygmin, zgmin, xgmax, ygmax, zgmax, unscale) + integer, intent(in) :: nk !< The number of vertical levels + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZI_(G),SZJ_(G),nk), intent(in) :: tr_array !< The tracer array to search for + !! extrema in arbitrary concentration units [CU ~> conc] + real, intent(out) :: g_min !< The global minimum of tr_array, either in + !! the same units as tr_array [CU ~> conc] or in + !! unscaled units if unscale is present [conc] + real, intent(out) :: g_max !< The global maximum of tr_array, either in + !! the same units as tr_array [CU ~> conc] or in + !! unscaled units if unscale is present [conc] + real, optional, intent(out) :: xgmin !< The x-position of the global minimum in the + !! units of G%geoLonT, often [degrees_E] or [km] or [m] + real, optional, intent(out) :: ygmin !< The y-position of the global minimum in the + !! units of G%geoLatT, often [degrees_N] or [km] or [m] + real, optional, intent(out) :: zgmin !< The z-position of the global minimum [layer] + real, optional, intent(out) :: xgmax !< The x-position of the global maximum in the + !! units of G%geoLonT, often [degrees_E] or [km] or [m] + real, optional, intent(out) :: ygmax !< The y-position of the global maximum in the + !! units of G%geoLatT, often [degrees_N] or [km] or [m] + real, optional, intent(out) :: zgmax !< The z-position of the global maximum [layer] + real, optional, intent(in) :: unscale !< A factor to use to undo any scaling of + !! the input tracer array [conc CU-1 ~> 1] + + ! Local variables + real :: tmax, tmin ! Maximum and minimum tracer values, in the same units as tr_array [CU ~> conc] + integer :: ijk_min_max(2) ! Integers encoding the global grid positions of the global minimum and maximum values + real :: xyz_min_max(6) ! A single array with the x-, y- and z-positions of the minimum and + ! maximum values in units that vary between the array elements [various] + logical :: valid_PE ! True if there are any valid points on the local PE. + logical :: find_location ! If true, report the locations of the extrema + integer :: ijk_loc_max ! An integer encoding the global grid position of the maximum tracer value on this PE + integer :: ijk_loc_min ! An integer encoding the global grid position of the minimum tracer value on this PE + integer :: ijk_loc_here ! An integer encoding the global grid position of the current grid point + integer :: itmax, jtmax, ktmax, itmin, jtmin, ktmin + integer :: i, j, k, isc, iec, jsc, jec + + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec + + find_location = (present(xgmin) .or. present(ygmin) .or. present(zgmin) .or. & + present(xgmax) .or. present(ygmax) .or. present(zgmax)) + + ! The initial values set here are never used if there are any valid points. + tmax = -huge(tmax) ; tmin = huge(tmin) + + if (find_location) then + ! Find the maximum and minimum tracer values on this PE and their locations. + valid_PE = .false. + itmax = 0 ; jtmax = 0 ; ktmax = 0 ; ijk_loc_max = 0 + itmin = 0 ; jtmin = 0 ; ktmin = 0 ; ijk_loc_min = 0 + do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then + valid_PE = .true. + if (tr_array(i,j,k) > tmax) then + tmax = tr_array(i,j,k) + itmax = i ; jtmax = j ; ktmax = k + ijk_loc_max = ijk_loc(i, j, k, nk, G%HI) + elseif ((tr_array(i,j,k) == tmax) .and. (k <= ktmax)) then + ijk_loc_here = ijk_loc(i, j, k, nk, G%HI) + if (ijk_loc_here > ijk_loc_max) then + itmax = i ; jtmax = j ; ktmax = k + ijk_loc_max = ijk_loc_here + endif + endif + if (tr_array(i,j,k) < tmin) then + tmin = tr_array(i,j,k) + itmin = i ; jtmin = j ; ktmin = k + ijk_loc_min = ijk_loc(i, j, k, nk, G%HI) + elseif ((tr_array(i,j,k) == tmin) .and. (k <= ktmin)) then + ijk_loc_here = ijk_loc(i, j, k, nk, G%HI) + if (ijk_loc_here > ijk_loc_min) then + itmin = i ; jtmin = j ; ktmin = k + ijk_loc_min = ijk_loc_here + endif + endif + endif ; enddo ; enddo ; enddo + else + ! Only the maximum and minimum values are needed, and not their positions. + do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then + if (tr_array(i,j,k) > tmax) tmax = tr_array(i,j,k) + if (tr_array(i,j,k) < tmin) tmin = tr_array(i,j,k) + endif ; enddo ; enddo ; enddo + endif + + ! Find the global maximum and minimum tracer values. + g_max = tmax ; g_min = tmin + call max_across_PEs(g_max) + call min_across_PEs(g_min) + + if (find_location) then + if (g_max < g_min) then + ! This only occurs if there are no unmasked points anywhere in the domain. + xyz_min_max(:) = 0.0 + else + ! Find the global indices of the maximum and minimum locations. This can + ! occur on multiple PEs. + ijk_min_max(1:2) = 0 + if (valid_PE) then + if (g_min == tmin) ijk_min_max(1) = ijk_loc_min + if (g_max == tmax) ijk_min_max(2) = ijk_loc_max + endif + ! If MOM6 supported taking maxima on arrays of integers, these could be combined as: + ! call max_across_PEs(ijk_min_max, 2) + call max_across_PEs(ijk_min_max(1)) + call max_across_PEs(ijk_min_max(2)) + + ! Set the positions of the extrema if they occur on this PE. This will only + ! occur on a single PE. + xyz_min_max(1:6) = -huge(xyz_min_max) ! These huge negative values are never selected by max_across_PEs. + if (valid_PE) then + if (ijk_min_max(1) == ijk_loc_min) then + xyz_min_max(1) = G%geoLonT(itmin,jtmin) + xyz_min_max(2) = G%geoLatT(itmin,jtmin) + xyz_min_max(3) = real(ktmin) + endif + if (ijk_min_max(2) == ijk_loc_max) then + xyz_min_max(4) = G%geoLonT(itmax,jtmax) + xyz_min_max(5) = G%geoLatT(itmax,jtmax) + xyz_min_max(6) = real(ktmax) + endif + endif + + call max_across_PEs(xyz_min_max, 6) + endif + + if (present(xgmin)) xgmin = xyz_min_max(1) + if (present(ygmin)) ygmin = xyz_min_max(2) + if (present(zgmin)) zgmin = xyz_min_max(3) + if (present(xgmax)) xgmax = xyz_min_max(4) + if (present(ygmax)) ygmax = xyz_min_max(5) + if (present(zgmax)) zgmax = xyz_min_max(6) + endif + + if (g_max < g_min) then + ! There are no unmasked points anywhere in the domain. + g_max = 0.0 ; g_min = 0.0 + endif + + if (present(unscale)) then + ! Rescale g_min and g_max, perhaps changing their units from [CU ~> conc] to [conc] + g_max = unscale * g_max + g_min = unscale * g_min + endif + +end subroutine array_global_min_max + +! Return a positive integer encoding the rotationally invariant global position of a tracer cell +function ijk_loc(i, j, k, nk, HI) + integer, intent(in) :: i !< Local i-index + integer, intent(in) :: j !< Local j-index + integer, intent(in) :: k !< Local k-index + integer, intent(in) :: nk !< Range of k-index, used to pick out a low-k position. + type(hor_index_type), intent(in) :: HI !< Horizontal index ranges + integer :: ijk_loc ! An integer encoding the cell position in the global grid. + + ! Local variables + integer :: ig, jg ! Global index values with a global computational domain start value of 1. + integer :: ij_loc ! The encoding of the horizontal position + integer :: qturns ! The number of counter-clockwise quarter turns of the grid that have to be undone + + ! These global i-grid positions run from 1 to HI%niglobal, and analogously for jg. + ig = i + HI%idg_offset + (1 - HI%isg) + jg = j + HI%jdg_offset + (1 - HI%jsg) + + ! Compensate for the rotation of the model grid to give a rotationally invariant encoding. + qturns = modulo(HI%turns, 4) + if (qturns == 0) then + ij_loc = ig + HI%niglobal * jg + elseif (qturns == 1) then + ij_loc = jg + HI%njglobal * ((HI%niglobal+1)-ig) + elseif (qturns == 2) then + ij_loc = ((HI%niglobal+1)-ig) + HI%niglobal * ((HI%njglobal+1)-jg) + elseif (qturns == 3) then + ij_loc = ((HI%njglobal+1)-jg) + HI%njglobal * ig + endif + + ijk_loc = ij_loc + (HI%niglobal*HI%njglobal) * (nk-k) + +end function ijk_loc + + end module MOM_spatial_means diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index e998e3029c..6b10d15e2f 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -44,7 +44,7 @@ module MOM_generic_tracer use MOM_open_boundary, only : register_obgc_segments, fill_obgc_segments use MOM_open_boundary, only : set_obgc_segments_props use MOM_restart, only : register_restart_field, query_initialized, set_initialized, MOM_restart_CS - use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP + use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP, array_global_min_max use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, set_time use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut @@ -67,7 +67,7 @@ module MOM_generic_tracer public end_MOM_generic_tracer, MOM_generic_tracer_get public MOM_generic_tracer_stock public MOM_generic_flux_init - public MOM_generic_tracer_min_max, array_global_min_max + public MOM_generic_tracer_min_max public MOM_generic_tracer_fluxes_accumulate public register_MOM_generic_tracer_segments @@ -802,192 +802,6 @@ function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, G, CS, na end function MOM_generic_tracer_min_max -!> Find the global maximum and minimum of a tracer array and return the locations of the extrema. -!! When there multiple cells with the same extreme values, the reported locations are from the -!! uppermost layer where they occur, and then from the logically northernmost and then eastermost -!! such location on the unrotated version of the grid within that layer. Only ocean points (as -!! indicated by a positive value of G%mask2dT) are evaluated, and if there are no ocean points -!! anywhere in the domain, the reported extrema and their locations are all returned as 0. - subroutine array_global_min_max(tr_array, G, nk, g_min, g_max, & - xgmin, ygmin, zgmin, xgmax, ygmax, zgmax, unscale) - integer, intent(in) :: nk !< The number of vertical levels - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - real, dimension(SZI_(G),SZJ_(G),nk), intent(in) :: tr_array !< The tracer array to search for - !! extrema in arbitrary concentration units [CU ~> conc] - real, intent(out) :: g_min !< The global minimum of tr_array, either in - !! the same units as tr_array [CU ~> conc] or in - !! unscaled units if unscale is present [conc] - real, intent(out) :: g_max !< The global maximum of tr_array, either in - !! the same units as tr_array [CU ~> conc] or in - !! unscaled units if unscale is present [conc] - real, optional, intent(out) :: xgmin !< The x-position of the global minimum in the - !! units of G%geoLonT, often [degrees_E] or [km] or [m] - real, optional, intent(out) :: ygmin !< The y-position of the global minimum in the - !! units of G%geoLatT, often [degrees_N] or [km] or [m] - real, optional, intent(out) :: zgmin !< The z-position of the global minimum [layer] - real, optional, intent(out) :: xgmax !< The x-position of the global maximum in the - !! units of G%geoLonT, often [degrees_E] or [km] or [m] - real, optional, intent(out) :: ygmax !< The y-position of the global maximum in the - !! units of G%geoLatT, often [degrees_N] or [km] or [m] - real, optional, intent(out) :: zgmax !< The z-position of the global maximum [layer] - real, optional, intent(in) :: unscale !< A factor to use to undo any scaling of - !! the input tracer array [conc CU-1 ~> 1] - - ! Local variables - real :: tmax, tmin ! Maximum and minimum tracer values, in the same units as tr_array [CU ~> conc] - integer :: ijk_min_max(2) ! Integers encoding the global grid positions of the global minimum and maximum values - real :: xyz_min_max(6) ! A single array with the x-, y- and z-positions of the minimum and - ! maximum values in units that vary between the array elements [various] - logical :: valid_PE ! True if there are any valid points on the local PE. - logical :: find_location ! If true, report the locations of the extrema - integer :: ijk_loc_max ! An integer encoding the global grid position of the maximum tracer value on this PE - integer :: ijk_loc_min ! An integer encoding the global grid position of the minimum tracer value on this PE - integer :: ijk_loc_here ! An integer encoding the global grid position of the current grid point - integer :: itmax, jtmax, ktmax, itmin, jtmin, ktmin - integer :: i, j, k, isc, iec, jsc, jec - - isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec - - find_location = (present(xgmin) .or. present(ygmin) .or. present(zgmin) .or. & - present(xgmax) .or. present(ygmax) .or. present(zgmax)) - - ! The initial values set here are never used if there are any valid points. - tmax = -huge(tmax) ; tmin = huge(tmin) - - if (find_location) then - ! Find the maximum and minimum tracer values on this PE and their locations. - valid_PE = .false. - itmax = 0 ; jtmax = 0 ; ktmax = 0 ; ijk_loc_max = 0 - itmin = 0 ; jtmin = 0 ; ktmin = 0 ; ijk_loc_min = 0 - do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then - valid_PE = .true. - if (tr_array(i,j,k) > tmax) then - tmax = tr_array(i,j,k) - itmax = i ; jtmax = j ; ktmax = k - ijk_loc_max = ijk_loc(i, j, k, nk, G%HI) - elseif ((tr_array(i,j,k) == tmax) .and. (k <= ktmax)) then - ijk_loc_here = ijk_loc(i, j, k, nk, G%HI) - if (ijk_loc_here > ijk_loc_max) then - itmax = i ; jtmax = j ; ktmax = k - ijk_loc_max = ijk_loc_here - endif - endif - if (tr_array(i,j,k) < tmin) then - tmin = tr_array(i,j,k) - itmin = i ; jtmin = j ; ktmin = k - ijk_loc_min = ijk_loc(i, j, k, nk, G%HI) - elseif ((tr_array(i,j,k) == tmin) .and. (k <= ktmin)) then - ijk_loc_here = ijk_loc(i, j, k, nk, G%HI) - if (ijk_loc_here > ijk_loc_min) then - itmin = i ; jtmin = j ; ktmin = k - ijk_loc_min = ijk_loc_here - endif - endif - endif ; enddo ; enddo ; enddo - else - ! Only the maximum and minimum values are needed, and not their positions. - do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then - if (tr_array(i,j,k) > tmax) tmax = tr_array(i,j,k) - if (tr_array(i,j,k) < tmin) tmin = tr_array(i,j,k) - endif ; enddo ; enddo ; enddo - endif - - ! Find the global maximum and minimum tracer values. - g_max = tmax ; g_min = tmin - call max_across_PEs(g_max) - call min_across_PEs(g_min) - - if (find_location) then - if (g_max < g_min) then - ! This only occurs if there are no unmasked points anywhere in the domain. - xyz_min_max(:) = 0.0 - else - ! Find the global indices of the maximum and minimum locations. This can - ! occur on multiple PEs. - ijk_min_max(1:2) = 0 - if (valid_PE) then - if (g_min == tmin) ijk_min_max(1) = ijk_loc_min - if (g_max == tmax) ijk_min_max(2) = ijk_loc_max - endif - ! If MOM6 supported taking maxima on arrays of integers, these could be combined as: - ! call max_across_PEs(ijk_min_max, 2) - call max_across_PEs(ijk_min_max(1)) - call max_across_PEs(ijk_min_max(2)) - - ! Set the positions of the extrema if they occur on this PE. This will only - ! occur on a single PE. - xyz_min_max(1:6) = -huge(xyz_min_max) ! These huge negative values are never selected by max_across_PEs. - if (valid_PE) then - if (ijk_min_max(1) == ijk_loc_min) then - xyz_min_max(1) = G%geoLonT(itmin,jtmin) - xyz_min_max(2) = G%geoLatT(itmin,jtmin) - xyz_min_max(3) = real(ktmin) - endif - if (ijk_min_max(2) == ijk_loc_max) then - xyz_min_max(4) = G%geoLonT(itmax,jtmax) - xyz_min_max(5) = G%geoLatT(itmax,jtmax) - xyz_min_max(6) = real(ktmax) - endif - endif - - call max_across_PEs(xyz_min_max, 6) - endif - - if (present(xgmin)) xgmin = xyz_min_max(1) - if (present(ygmin)) ygmin = xyz_min_max(2) - if (present(zgmin)) zgmin = xyz_min_max(3) - if (present(xgmax)) xgmax = xyz_min_max(4) - if (present(ygmax)) ygmax = xyz_min_max(5) - if (present(zgmax)) zgmax = xyz_min_max(6) - endif - - if (g_max < g_min) then - ! There are no unmasked points anywhere in the domain. - g_max = 0.0 ; g_min = 0.0 - endif - - if (present(unscale)) then - ! Rescale g_min and g_max, perhaps changing their units from [CU ~> conc] to [conc] - g_max = unscale * g_max - g_min = unscale * g_min - endif - - end subroutine array_global_min_max - - ! Return a positive integer encoding the rotationally invariant global position of a tracer cell - function ijk_loc(i, j, k, nk, HI) - integer, intent(in) :: i !< Local i-index - integer, intent(in) :: j !< Local j-index - integer, intent(in) :: k !< Local k-index - integer, intent(in) :: nk !< Range of k-index, used to pick out a low-k position. - type(hor_index_type), intent(in) :: HI !< Horizontal index ranges - integer :: ijk_loc ! An integer encoding the cell position in the global grid. - - ! Local variables - integer :: ig, jg ! Global index values with a global computational domain start value of 1. - integer :: ij_loc ! The encoding of the horizontal position - integer :: qturns ! The number of counter-clockwise quarter turns of the grid that have to be undone - - ! These global i-grid positions run from 1 to HI%niglobal, and analogously for jg. - ig = i + HI%idg_offset + (1 - HI%isg) - jg = j + HI%jdg_offset + (1 - HI%jsg) - - ! Compensate for the rotation of the model grid to give a rotationally invariant encoding. - qturns = modulo(HI%turns, 4) - if (qturns == 0) then - ij_loc = ig + HI%niglobal * jg - elseif (qturns == 1) then - ij_loc = jg + HI%njglobal * ((HI%niglobal+1)-ig) - elseif (qturns == 2) then - ij_loc = ((HI%niglobal+1)-ig) + HI%niglobal * ((HI%njglobal+1)-jg) - elseif (qturns == 3) then - ij_loc = ((HI%njglobal+1)-jg) + HI%njglobal * ig - endif - - ijk_loc = ij_loc + (HI%niglobal*HI%njglobal) * (nk-k) - - end function ijk_loc - !> This subroutine calculates the surface state and sets coupler values for !! those generic tracers that have flux exchange with atmosphere. !! From 7e394219e1b227bd9bffb680e0cd1c63706daec2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 10 Jul 2024 18:36:45 -0400 Subject: [PATCH 2/2] +Add runtime parameter WRITE_TRACER_MIN_MAX Added the new runtime parameters WRITE_TRACER_MIN_MAX and WRITE_TRACER_MIN_MAX_LOC to the MOM_sum_output module to control whether the maximum and minimum values of temperature, salinity and some tracers are periodically written to stdout, perhaps with their locations, as determined by array_global_min_max. These can be expensive global reductions, so by default these are set to false. All solutions are bitwise identical, but there can be some changes to the output to stdout and there will be new entries in the MOM_parameter_doc.debugging files. --- src/diagnostics/MOM_sum_output.F90 | 108 +++++++++++++++++++++++++---- 1 file changed, 94 insertions(+), 14 deletions(-) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index edb66d225c..a97df8e94d 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -20,6 +20,7 @@ module MOM_sum_output use MOM_io, only : axis_info, set_axis_info, delete_axis_info, get_filename_appendix use MOM_io, only : attribute_info, set_attribute_info, delete_attribute_info use MOM_io, only : APPEND_FILE, SINGLE_FILE, WRITEONLY_FILE +use MOM_spatial_means, only : array_global_min_max use MOM_time_manager, only : time_type, get_time, get_date, set_time, operator(>) use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) use MOM_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<) @@ -124,6 +125,12 @@ module MOM_sum_output !! interval at which the run is stopped. logical :: write_stocks !< If true, write the integrated tracer amounts !! to stdout when the energy files are written. + logical :: write_min_max !< If true, write the maximum and minimum values of temperature, + !! salinity and some tracer concentrations to stdout when the energy + !! files are written. + logical :: write_min_max_loc !< If true, write the locations of the maximum and minimum values + !! of temperature, salinity and some tracer concentrations to stdout + !! when the energy files are written. integer :: previous_calls = 0 !< The number of times write_energy has been called. integer :: prev_n = 0 !< The value of n from the last call. type(MOM_netcdf_file) :: fileenergy_nc !< The file handle for the netCDF version of the energy file. @@ -179,6 +186,15 @@ subroutine MOM_sum_output_init(G, GV, US, param_file, directory, ntrnc, & call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", CS%use_temperature, & "If true, Temperature and salinity are used as state "//& "variables.", default=.true.) + call get_param(param_file, mdl, "WRITE_TRACER_MIN_MAX", CS%write_min_max, & + "If true, write the maximum and minimum values of temperature, salinity and "//& + "some tracer concentrations to stdout when the energy files are written.", & + default=.false., do_not_log=.not.CS%write_stocks, debuggingParam=.true.) + call get_param(param_file, mdl, "WRITE_TRACER_MIN_MAX_LOC", CS%write_min_max_loc, & + "If true, write the locations of the maximum and minimum values of "//& + "temperature, salinity and some tracer concentrations to stdout when the "//& + "energy files are written.", & + default=.false., do_not_log=.not.CS%write_min_max, debuggingParam=.true.) call get_param(param_file, mdl, "DT", CS%dt_in_T, & "The (baroclinic) dynamics time step.", & units="s", scale=US%s_to_T, fail_if_missing=.true.) @@ -404,6 +420,34 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci character(len=32) :: mesg_intro, time_units, day_str, n_str, date_str logical :: date_stamped type(time_type) :: dt_force ! A time_type version of the forcing timestep. + + real :: S_min ! The global minimum unmasked value of the salinity [ppt] + real :: S_max ! The global maximum unmasked value of the salinity [ppt] + real :: S_min_x ! The x-positions of the global salinity minima + ! in the units of G%geoLonT, often [degrees_E] or [km] + real :: S_min_y ! The y-positions of the global salinity minima + ! in the units of G%geoLatT, often [degrees_N] or [km] + real :: S_min_z ! The z-positions of the global salinity minima [layer] + real :: S_max_x ! The x-positions of the global salinity maxima + ! in the units of G%geoLonT, often [degrees_E] or [km] + real :: S_max_y ! The y-positions of the global salinity maxima + ! in the units of G%geoLatT, often [degrees_N] or [km] + real :: S_max_z ! The z-positions of the global salinity maxima [layer] + + real :: T_min ! The global minimum unmasked value of the temperature [degC] + real :: T_max ! The global maximum unmasked value of the temperature [degC] + real :: T_min_x ! The x-positions of the global temperature minima + ! in the units of G%geoLonT, often [degreeT_E] or [km] + real :: T_min_y ! The y-positions of the global temperature minima + ! in the units of G%geoLatT, often [degreeT_N] or [km] + real :: T_min_z ! The z-positions of the global temperature minima [layer] + real :: T_max_x ! The x-positions of the global temperature maxima + ! in the units of G%geoLonT, often [degreeT_E] or [km] + real :: T_max_y ! The y-positions of the global temperature maxima + ! in the units of G%geoLatT, often [degreeT_N] or [km] + real :: T_max_z ! The z-positions of the global temperature maxima [layer] + + ! The units of the tracer stock vary between tracers, with [conc] given explicitly by Tr_units. real :: Tr_stocks(MAX_FIELDS_) ! The total amounts of each of the registered tracers [kg conc] real :: Tr_min(MAX_FIELDS_) ! The global minimum unmasked value of the tracers [conc] @@ -527,11 +571,20 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci nTr_stocks = 0 Tr_minmax_avail(:) = .false. - call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, & - stock_units=Tr_units, num_stocks=nTr_stocks,& - got_min_max=Tr_minmax_avail, global_min=Tr_min, global_max=Tr_max, & - xgmin=Tr_min_x, ygmin=Tr_min_y, zgmin=Tr_min_z,& - xgmax=Tr_max_x, ygmax=Tr_max_y, zgmax=Tr_max_z) + if (CS%write_min_max .and. CS%write_min_max_loc) then + call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, & + stock_units=Tr_units, num_stocks=nTr_stocks,& + got_min_max=Tr_minmax_avail, global_min=Tr_min, global_max=Tr_max, & + xgmin=Tr_min_x, ygmin=Tr_min_y, zgmin=Tr_min_z,& + xgmax=Tr_max_x, ygmax=Tr_max_y, zgmax=Tr_max_z) + elseif (CS%write_min_max) then + call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, & + stock_units=Tr_units, num_stocks=nTr_stocks,& + got_min_max=Tr_minmax_avail, global_min=Tr_min, global_max=Tr_max) + else + call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, & + stock_units=Tr_units, num_stocks=nTr_stocks) + endif if (nTr_stocks > 0) then do m=1,nTr_stocks vars(num_nc_fields+m) = var_desc(Tr_names(m), units=Tr_units(m), & @@ -540,6 +593,13 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci num_nc_fields = num_nc_fields + nTr_stocks endif + if (CS%use_temperature .and. CS%write_stocks) then + call array_global_min_max(tv%T, G, nz, T_min, T_max, & + T_min_x, T_min_y, T_min_z, T_max_x, T_max_y, T_max_z, unscale=US%C_to_degC) + call array_global_min_max(tv%S, G, nz, S_min, S_max, & + S_min_x, S_min_y, S_min_z, S_max_x, S_max_y, S_max_z, unscale=US%S_to_ppt) + endif + if (CS%previous_calls == 0) then CS%mass_prev_EFP = mass_EFP @@ -847,6 +907,15 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci write(stdout,'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & Salt*0.001, Salt_chg*0.001, Salt_anom*0.001, Salt_anom/Salt endif + if (CS%write_min_max .and. CS%write_min_max_loc) then + write(stdout,'(16X,"Salinity Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & + S_min, S_min_x, S_min_y, S_min_z + write(stdout,'(16X,"Salinity Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & + S_max, S_max_x, S_max_y, S_max_z + elseif (CS%write_min_max) then + write(stdout,'(16X,"Salinity Global Min & Max:",ES24.16,1X,ES24.16)') S_min, S_max + endif + if (Heat == 0.) then write(stdout,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') & Heat, Heat_chg, Heat_anom @@ -854,17 +923,28 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci write(stdout,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & Heat, Heat_chg, Heat_anom, Heat_anom/Heat endif + if (CS%write_min_max .and. CS%write_min_max_loc) then + write(stdout,'(16X,"Temperature Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & + T_min, T_min_x, T_min_y, T_min_z + write(stdout,'(16X,"Temperature Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & + T_max, T_max_x, T_max_y, T_max_z + elseif (CS%write_min_max) then + write(stdout,'(16X,"Temperature Global Min & Max:",ES24.16,1X,ES24.16)') T_min, T_max + endif endif do m=1,nTr_stocks - write(stdout,'(" Total ",a,": ",ES24.16,1X,a)') & + write(stdout,'(" Total ",a,": ",ES24.16,1X,a)') & trim(Tr_names(m)), Tr_stocks(m), trim(Tr_units(m)) - if (Tr_minmax_avail(m)) then - write(stdout,'(64X,"Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & - Tr_min(m),Tr_min_x(m),Tr_min_y(m),Tr_min_z(m) - write(stdout,'(64X,"Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & - Tr_max(m),Tr_max_x(m),Tr_max_y(m),Tr_max_z(m) + if (CS%write_min_max .and. CS%write_min_max_loc .and. Tr_minmax_avail(m)) then + write(stdout,'(18X,a," Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & + trim(Tr_names(m)), Tr_min(m), Tr_min_x(m), Tr_min_y(m), Tr_min_z(m) + write(stdout,'(18X,a," Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & + trim(Tr_names(m)), Tr_max(m), Tr_max_x(m), Tr_max_y(m), Tr_max_z(m) + elseif (CS%write_min_max .and. Tr_minmax_avail(m)) then + write(stdout,'(18X,a," Global Min & Max:",ES24.16,1X,ES24.16)') & + trim(Tr_names(m)), Tr_min(m), Tr_max(m) endif enddo @@ -1269,9 +1349,9 @@ subroutine write_depth_list(G, US, DL, filename) call create_MOM_file(IO_handle, filename, vars, 3, fields, SINGLE_FILE, & extra_axes=extra_axes, global_atts=global_atts) - call MOM_write_field(IO_handle, fields(1), DL%depth, unscale=US%Z_to_m) - call MOM_write_field(IO_handle, fields(2), DL%area, unscale=US%L_to_m**2) - call MOM_write_field(IO_handle, fields(3), DL%vol_below, unscale=US%Z_to_m*US%L_to_m**2) + call MOM_write_field(IO_handle, fields(1), DL%depth, scale=US%Z_to_m) + call MOM_write_field(IO_handle, fields(2), DL%area, scale=US%L_to_m**2) + call MOM_write_field(IO_handle, fields(3), DL%vol_below, scale=US%Z_to_m*US%L_to_m**2) call delete_axis_info(extra_axes) call delete_attribute_info(global_atts)