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/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) 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. !!