Merge pull request #1319 from Hallberg-NOAA/generic_tracer_minmax

Included array_global_min_man in MOM_generic_tracer.F90

src/tracer/MOM_generic_tracer.F90

@@ -28,6 +28,8 @@ module MOM_generic_tracer
   use g_tracer_utils,   only: g_tracer_send_diag,g_tracer_get_values
   use g_tracer_utils,   only: g_tracer_get_pointer,g_tracer_get_alias,g_tracer_set_csdiag
 
+  use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
+  use MOM_coms, only : max_across_PEs, min_across_PEs, PE_here
   use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
   use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end
   use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe
@@ -36,18 +38,17 @@ module MOM_generic_tracer
   use MOM_grid, only : ocean_grid_type
   use MOM_hor_index, only : hor_index_type
   use MOM_io, only : file_exists, MOM_read_data, slasher
+  use MOM_open_boundary, only : ocean_OBC_type
   use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS
   use MOM_spatial_means, only : global_area_mean
   use MOM_sponge, only : set_up_sponge_field, sponge_CS
-  use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
   use MOM_time_manager, only : time_type, set_time
   use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
   use MOM_tracer_registry, only : register_tracer, tracer_registry_type
   use MOM_tracer_Z_init, only : tracer_Z_init
   use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z
   use MOM_unit_scaling, only : unit_scale_type
   use MOM_variables, only : surface, thermo_var_ptrs
-  use MOM_open_boundary, only : ocean_OBC_type
   use MOM_verticalGrid, only : verticalGrid_type
 
 
@@ -633,7 +634,6 @@ end function MOM_generic_tracer_stock
   !! requested specifically, returning the number of tracers it has gone through.
   function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, xgmin, ygmin, zgmin, &
                                       xgmax, ygmax, zgmax , G, CS, names, units)
-    use mpp_utilities_mod, only: mpp_array_global_min_max
     integer,                        intent(in)    :: ind_start !< The index of the tracer to start with
     logical, dimension(:),          intent(out)   :: got_minmax !< Indicates whether the global min and
                                                             !! max are found for each tracer
@@ -693,8 +693,8 @@ function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, xgmin, yg
 
       tr_ptr => tr_field(:,:,:,1)
 
-      call mpp_array_global_min_max(tr_ptr, grid_tmask,isd,jsd,isc,iec,jsc,jec,nk , gmin(m), gmax(m), &
-                                    G%geoLonT,G%geoLatT,geo_z,xgmin(m), ygmin(m), zgmin(m), &
+      call array_global_min_max(tr_ptr, grid_tmask, isd, jsd, isc, iec, jsc, jec, nk, gmin(m), gmax(m), &
+                                    G%geoLonT, G%geoLatT, geo_z, xgmin(m), ygmin(m), zgmin(m), &
                                     xgmax(m), ygmax(m), zgmax(m))
 
       got_minmax(m) = .true.
@@ -710,6 +710,133 @@ function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, xgmin, yg
 
   end function MOM_generic_tracer_min_max
 
+  !> Find the global maximum and minimum of a tracer array and return the locations of the extrema.
+  subroutine array_global_min_max(tr_array, tmask, isd, jsd, isc, iec, jsc, jec, nk, g_min, g_max, &
+                                  geo_x, geo_y, geo_z, xgmin, ygmin, zgmin, xgmax, ygmax, zgmax)
+    integer,                      intent(in)  :: isd   !< The starting data domain i-index
+    integer,                      intent(in)  :: jsd   !< The starting data domain j-index
+    real, dimension(isd:,jsd:,:), intent(in)  :: tr_array !< The tracer array to search for extrema
+    real, dimension(isd:,jsd:,:), intent(in)  :: tmask !< A mask that is 0 for points to exclude
+    integer,                      intent(in)  :: isc   !< The starting compute domain i-index
+    integer,                      intent(in)  :: iec   !< The ending compute domain i-index
+    integer,                      intent(in)  :: jsc   !< The starting compute domain j-index
+    integer,                      intent(in)  :: jec   !< The ending compute domain j-index
+    integer,                      intent(in)  :: nk    !< The number of vertical levels
+    real,                         intent(out) :: g_min !< The global minimum of tr_array
+    real,                         intent(out) :: g_max !< The global maximum of tr_array
+    real, dimension(isd:,jsd:),   intent(in)  :: geo_x !< The geographic x-positions of points
+    real, dimension(isd:,jsd:),   intent(in)  :: geo_y !< The geographic y-positions of points
+    real, dimension(:),           intent(in)  :: geo_z !< The vertical pseudo-positions of points
+    real,                         intent(out) :: xgmin !< The x-position of the global minimum
+    real,                         intent(out) :: ygmin !< The y-position of the global minimum
+    real,                         intent(out) :: zgmin !< The z-position of the global minimum
+    real,                         intent(out) :: xgmax !< The x-position of the global maximum
+    real,                         intent(out) :: ygmax !< The y-position of the global maximum
+    real,                         intent(out) :: zgmax !< The z-position of the global maximum
+
+    ! This subroutine is an exact transcription (bugs and all) of mpp_array_global_min_max()
+    ! from the version in FMS/mpp/mpp_utilities.F90, but with some whitespace changes to match
+    ! MOM6 code styles and to use infrastructure routines via the MOM6 framework code, and with
+    ! added comments to document its arguments.i
+
+    !### The obvious problems with this routine as currently written include:
+    !  1. It does not return exactly the maximum and minimum values.
+    !  2. The reported maximum and minimum are dependent on PE count and layout.
+    !  3. For all-zero arrays, the reported maxima scale with the PE_count
+    !  4. For arrays with a large enough offset or scaling, so that the magnitude of values exceed
+    !     1e10, the values it returns are simply wrong.
+    !  5. The results do not scale appropriately if the argument is rescaled.
+    !  6. The extrema and locations are not rotationally invariant.
+    !  7. It is inefficient because it uses 8 blocking global reduction calls when it could use just 2 or 3.
+
+    ! Local variables
+    real    :: tmax, tmin   ! Maximum and minimum tracer values, in the same units as tr_array
+    real    :: tmax0, tmin0 ! First-guest values of tmax and tmin.
+    integer :: itmax, jtmax, ktmax, itmin, jtmin, ktmin
+    integer :: igmax, jgmax, kgmax, igmin, jgmin, kgmin
+    real    :: fudge ! A factor that is close to 1 that is used to find the location of the extrema.
+
+     ! arrays to enable vectorization
+    integer :: iminarr(3), imaxarr(3)
+
+    !### These dimensional constant values mean that the results can not be guaranteed to be rescalable.
+    g_min = -88888888888.0 ; g_max = -999999999.0
+    tmax = -1.e10 ; tmin = 1.e10
+    itmax = 0 ; jtmax = 0 ; ktmax = 0
+    itmin = 0 ; jtmin = 0 ; ktmin = 0
+
+    if (ANY(tmask(isc:iec,jsc:jec,:) > 0.)) then
+      ! Vectorized using maxloc() and minloc() intrinsic functions by Russell.Fiedler@csiro.au.
+      iminarr = minloc(tr_array(isc:iec,jsc:jec,:), (tmask(isc:iec,jsc:jec,:) > 0.))
+      imaxarr = maxloc(tr_array(isc:iec,jsc:jec,:), (tmask(isc:iec,jsc:jec,:) > 0.))
+      itmin = iminarr(1)+isc-1
+      jtmin = iminarr(2)+jsc-1
+      ktmin = iminarr(3)
+      itmax = imaxarr(1)+isc-1
+      jtmax = imaxarr(2)+jsc-1
+      ktmax = imaxarr(3)
+      tmin = tr_array(itmin,jtmin,ktmin)
+      tmax = tr_array(itmax,jtmax,ktmax)
+    end if
+
+    ! use "fudge" to distinguish processors when tracer extreme is independent of processor
+    !### This fudge factor is not independent of PE layout, and while it mostly works for finding
+    !    a positive maximum or a negative minimum, it could miss the true extrema in the opposite
+    !    cases, for which the fudge factor should be slightly reduced.  The fudge factor should
+    !    be based on global index-space conventions, which are decomposition invariant, and
+    !    not the PE-number!
+    fudge = 1.0 + 1.e-12*real(PE_here() )
+    tmax = tmax*fudge
+    tmin = tmin*fudge
+    if (tmax == 0.0) then
+      tmax = tmax + 1.e-12*real(PE_here() )
+    endif
+    if (tmin == 0.0) then
+      tmin = tmin + 1.e-12*real(PE_here() )
+    endif
+
+    tmax0 = tmax ; tmin0 = tmin
+
+    call max_across_PEs(tmax)
+    call min_across_PEs(tmin)
+
+    g_max = tmax
+    g_min = tmin
+
+    ! Now find the location of the global extrema.
+    !
+    ! Note that the fudge factor above guarantees that the location of max (min) is uinque,
+    ! since tmax0 (tmin0) has slightly different values on each processor.
+    ! Otherwise, the function tr_array(i,j,k) could be equal to global max (min) at more
+    ! than one point in space and this would be a much more difficult problem to solve.
+    !
+    !-999 on all current PE's
+    xgmax = -999. ; ygmax = -999. ; zgmax = -999.
+    xgmin = -999. ; ygmin = -999. ; zgmin = -999.
+
+    if (tmax0 == tmax) then !This happens ONLY on ONE processor because of fudge factor above.
+      xgmax = geo_x(itmax,jtmax)
+      ygmax = geo_y(itmax,jtmax)
+      zgmax = geo_z(ktmax)
+    endif
+
+    !### These three calls and the three calls that follow in about 10 lines should be combined
+    !    into a single call for efficiency.
+    call max_across_PEs(xgmax)
+    call max_across_PEs(ygmax)
+    call max_across_PEs(zgmax)
+
+    if (tmin0 == tmin) then !This happens ONLY on ONE processor because of fudge factor above.
+      xgmin = geo_x(itmin,jtmin)
+      ygmin = geo_y(itmin,jtmin)
+      zgmin = geo_z(ktmin)
+    endif
+
+    call max_across_PEs(xgmin)
+    call max_across_PEs(ygmin)
+    call max_across_PEs(zgmin)
+
+  end subroutine array_global_min_max
 
   !> This subroutine calculates the surface state and sets coupler values for
   !! those generic tracers that have flux exchange with atmosphere.