Merge branch 'dev/gfdl' into makedep_cleanup

Author: marshallward

src/core/MOM_barotropic.F90

@@ -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
@@ -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]
@@ -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", &
@@ -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
 

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