Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

+Move array_global_min_max() and use it for Temperature & Salinity #683

Merged
merged 2 commits into from
Aug 10, 2024
Merged

+Move array_global_min_max() and use it for Temperature & Salinity #683

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8658edf
+Move array_global_min_max to MOM_spatial_means
Hallberg-NOAA Jul 10, 2024
7e39421
+Add runtime parameter WRITE_TRACER_MIN_MAX
Hallberg-NOAA Jul 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
191 changes: 191 additions & 0 deletions src/diagnostics/MOM_spatial_means.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Loading
Loading