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

Included array_global_min_man in MOM_generic_tracer.F90 #1319

Merged
merged 2 commits into from
Feb 10, 2021
Merged

Included array_global_min_man in MOM_generic_tracer.F90 #1319

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9a315f2
Included array_global_min_man in MOM_generic_tracer.F90
Hallberg-NOAA Feb 8, 2021
d62ca60
Merge branch 'dev/gfdl' into generic_tracer_minmax
adcroft Feb 10, 2021
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
137 changes: 132 additions & 5 deletions src/tracer/MOM_generic_tracer.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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


Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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.
Expand Down