+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

@@ -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