Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into makedep_cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Aug 27, 2024
2 parents 117f16b + d234bce commit 8a2b5d7
Show file tree
Hide file tree
Showing 6 changed files with 353 additions and 234 deletions.
54 changes: 40 additions & 14 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ module MOM_barotropic
use MOM_grid, only : ocean_grid_type
use MOM_harmonic_analysis, only : HA_accum_FtSSH, harmonic_analysis_CS
use MOM_hor_index, only : hor_index_type
use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher
use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher, NORTH_FACE, EAST_FACE
use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, open_boundary_query
use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_segment_type
Expand Down Expand Up @@ -4459,6 +4459,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
! drag piston velocity.
character(len=80) :: wave_drag_var ! The wave drag piston velocity variable
! name in wave_drag_file.
character(len=80) :: wave_drag_u ! The wave drag piston velocity variable
! name in wave_drag_file.
character(len=80) :: wave_drag_v ! The wave drag piston velocity variable
! name in wave_drag_file.
real :: mean_SL ! The mean sea level that is used along with the bathymetry to estimate the
! geometry when LINEARIZED_BT_CORIOLIS is true or BT_NONLIN_STRESS is false [Z ~> m].
real :: Z_to_H ! A local unit conversion factor [H Z-1 ~> nondim or kg m-3]
Expand Down Expand Up @@ -4703,8 +4707,17 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
"piston velocities.", default="", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_VAR", wave_drag_var, &
"The name of the variable in BT_WAVE_DRAG_FILE with the "//&
"barotropic linear wave drag piston velocities at h points.", &
"barotropic linear wave drag piston velocities at h points. "//&
"It will not be used if both BT_WAVE_DRAG_U and BT_WAVE_DRAG_V are defined.", &
default="rH", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_U", wave_drag_u, &
"The name of the variable in BT_WAVE_DRAG_FILE with the "//&
"barotropic linear wave drag piston velocities at u points.", &
default="", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_V", wave_drag_v, &
"The name of the variable in BT_WAVE_DRAG_FILE with the "//&
"barotropic linear wave drag piston velocities at v points.", &
default="", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_SCALE", wave_drag_scale, &
"A scaling factor for the barotropic linear wave drag "//&
"piston velocities.", default=1.0, units="nondim", &
Expand Down Expand Up @@ -4924,19 +4937,32 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
wave_drag_file = trim(slasher(inputdir))//trim(wave_drag_file)
call log_param(param_file, mdl, "INPUTDIR/BT_WAVE_DRAG_FILE", wave_drag_file)

allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0)
if (len_trim(wave_drag_u) > 0 .and. len_trim(wave_drag_v) > 0) then
call MOM_read_data(wave_drag_file, wave_drag_u, CS%lin_drag_u, G%Domain, &
position=EAST_FACE, scale=GV%m_to_H*US%T_to_s)
call pass_var(CS%lin_drag_u, G%Domain)
CS%lin_drag_u(:,:) = wave_drag_scale * CS%lin_drag_u(:,:)

call MOM_read_data(wave_drag_file, wave_drag_var, lin_drag_h, G%Domain, scale=GV%m_to_H*US%T_to_s)
call pass_var(lin_drag_h, G%Domain)
do j=js,je ; do I=is-1,ie
CS%lin_drag_u(I,j) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j))
enddo ; enddo
do J=js-1,je ; do i=is,ie
CS%lin_drag_v(i,J) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1))
enddo ; enddo
deallocate(lin_drag_h)
endif
endif
call MOM_read_data(wave_drag_file, wave_drag_v, CS%lin_drag_v, G%Domain, &
position=NORTH_FACE, scale=GV%m_to_H*US%T_to_s)
call pass_var(CS%lin_drag_v, G%Domain)
CS%lin_drag_v(:,:) = wave_drag_scale * CS%lin_drag_v(:,:)

else
allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0)

call MOM_read_data(wave_drag_file, wave_drag_var, lin_drag_h, G%Domain, scale=GV%m_to_H*US%T_to_s)
call pass_var(lin_drag_h, G%Domain)
do j=js,je ; do I=is-1,ie
CS%lin_drag_u(I,j) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j))
enddo ; enddo
do J=js-1,je ; do i=is,ie
CS%lin_drag_v(i,J) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1))
enddo ; enddo
deallocate(lin_drag_h)
endif ! len_trim(wave_drag_u) > 0 .and. len_trim(wave_drag_v) > 0
endif ! len_trim(wave_drag_file) > 0
endif ! CS%linear_wave_drag

CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input

Expand Down
191 changes: 191 additions & 0 deletions src/diagnostics/MOM_spatial_means.F90
View file
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

0 comments on commit 8a2b5d7

Please sign in to comment.