Skip to content

Commit

Permalink
(*)Use por_face_area in zonal_face_thickness
Browse files Browse the repository at this point in the history 
  Use the por_face_area[UV] in the effective thickness calculations in
zonal_face_thickness and merid_face_thickness, so that they are more consistent
with their use elsewhere in the code for the relative weights in calculating the
barotropic accelerations.  Because these por_face_area arrays are still 1 in all
test cases, the answers are unchanged in any test cases from before a few weeks
ago, but there could be answer changes in cases that are using the very recently
added capability (in PR #3) to set fractional face areas.  This change was
discussed with Sam Ditkovsky, and agreed that there is no reason to keep the
ability to recover the previous answers in any cases that use the recently added
partial face width option.

  This commit also expanded the comments describing the h_u and h_v arguments to
btcalc(), zonal_face_thickness(), and merid_face_thickness() routines, the
diag_h[uv] elements of the accel_diag_ptrs type and the h_u and h_v elements of
the BT_cont_type.

  All answers and output are bitwise identical in the MOM6-examples test suite
and TC tests, but answer changes are possible in cases using a very recently
added code option.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Dec 29, 2021
1 parent 1028ee0 commit 2b2214d
Show file tree
Hide file tree
Showing 3 changed files with 63 additions and 36 deletions.
20 changes: 15 additions & 5 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3274,9 +3274,19 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
optional, intent(in) :: h_u !< The specified thicknesses at u-points [H ~> m or kg m-2].
optional, intent(in) :: h_u !< The specified effective thicknesses at u-points,
!! perhaps scaled down to account for viscosity and
!! fractional open areas [H ~> m or kg m-2]. These
!! are used here as non-normalized weights for each
!! layer that are converted the normalized weights
!! for determining the barotropic accelerations.
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
optional, intent(in) :: h_v !< The specified thicknesses at v-points [H ~> m or kg m-2].
optional, intent(in) :: h_v !< The specified effective thicknesses at v-points,
!! perhaps scaled down to account for viscosity and
!! fractional open areas [H ~> m or kg m-2]. These
!! are used here as non-normalized weights for each
!! layer that are converted the normalized weights
!! for determining the barotropic accelerations.
logical, optional, intent(in) :: may_use_default !< An optional logical argument
!! to indicate that the default velocity point
!! thicknesses may be used for this particular
Expand All @@ -3296,9 +3306,9 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: wt_arith ! The weight for the arithmetic mean thickness [nondim].
! The harmonic mean uses a weight of (1 - wt_arith).
real :: Rh ! A ratio of summed thicknesses, nondim.
real :: e_u(SZIB_(G),SZK_(GV)+1) ! The interface heights at u-velocity and
real :: e_v(SZI_(G),SZK_(GV)+1) ! v-velocity points [H ~> m or kg m-2].
real :: Rh ! A ratio of summed thicknesses [nondim]
real :: e_u(SZIB_(G),SZK_(GV)+1) ! The interface heights at u-velocity points [H ~> m or kg m-2]
real :: e_v(SZI_(G),SZK_(GV)+1) ! The interface heights at v-velocity points [H ~> m or kg m-2]
real :: D_shallow_u(SZI_(G)) ! The height of the shallower of the adjacent bathymetric depths
! around a u-point (positive upward) [H ~> m or kg m-2]
real :: D_shallow_v(SZIB_(G))! The height of the shallower of the adjacent bathymetric depths
Expand Down
39 changes: 24 additions & 15 deletions src/core/MOM_continuity_PPM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -604,7 +604,8 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, &
endif
end subroutine zonal_flux_layer

!> Sets the effective interface thickness at each zonal velocity point.
!> Sets the effective interface thickness at each zonal velocity point, optionally scaling
!! back these thicknesses to account for viscosity and fractional open areas.
subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, &
marginal, OBC, por_face_areaU, visc_rem_u)
type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
Expand All @@ -616,7 +617,10 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL,
!! reconstruction [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_R !< Right thickness in the
!! reconstruction [H ~> m or kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_u !< Thickness at zonal faces [H ~> m or kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_u !< Effective thickness at zonal faces,
!! scaled down to account for the effects of
!! viscoity and the fractional open area
!! [H ~> m or kg m-2].
real, intent(in) :: dt !< Time increment [T ~> s].
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
Expand Down Expand Up @@ -672,11 +676,12 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL,
else ; h_u(I,j,k) = h_avg ; endif
enddo ; enddo ; enddo
if (present(visc_rem_u)) then
!### The expression setting h_u should also be multiplied by por_face_areaU in this case,
! and in the two OBC cases below with visc_rem_u.
! Scale back the thickness to account for the effects of viscosity and the fractional open
! thickness to give an appropriate non-normalized weight for each layer in determining the
! barotropic acceleration.
!$OMP parallel do default(shared)
do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh
h_u(I,j,k) = h_u(I,j,k) * visc_rem_u(I,j,k) !### * por_face_areaU(I,j,k)
h_u(I,j,k) = h_u(I,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k))
enddo ; enddo ; enddo
endif

Expand All @@ -689,7 +694,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL,
if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
if (present(visc_rem_u)) then ; do k=1,nz
do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
h_u(I,j,k) = h(i,j,k) * visc_rem_u(I,j,k) !### * por_face_areaU(I,j,k)
h_u(I,j,k) = h(i,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k))
enddo
enddo ; else ; do k=1,nz
do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
Expand All @@ -699,7 +704,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL,
else
if (present(visc_rem_u)) then ; do k=1,nz
do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
h_u(I,j,k) = h(i+1,j,k) * visc_rem_u(I,j,k) !### * por_face_areaU(I,j,k)
h_u(I,j,k) = h(i+1,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k))
enddo
enddo ; else ; do k=1,nz
do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
Expand Down Expand Up @@ -1427,19 +1432,22 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, &
endif
end subroutine merid_flux_layer

!> Sets the effective interface thickness at each meridional velocity point.
!> Sets the effective interface thickness at each meridional velocity point, optionally scaling
!! back these thicknesses to account for viscosity and fractional open areas.
subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, &
marginal, OBC, por_face_areaV, visc_rem_v)
type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness used to calculate fluxes,
!! [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_L !< Left thickness in the reconstruction,
!! [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_R !< Right thickness in the reconstruction,
!! [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: h_v !< Thickness at meridional faces,
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: h_v !< Effective thickness at meridional faces,
!! scaled down to account for the effects of
!! viscoity and the fractional open area
!! [H ~> m or kg m-2].
real, intent(in) :: dt !< Time increment [T ~> s].
type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
Expand Down Expand Up @@ -1497,11 +1505,12 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL,
enddo ; enddo ; enddo

if (present(visc_rem_v)) then
!### This expression setting h_v should also be multiplied by por_face_areaU in this case,
! and in the two OBC cases below with visc_rem_u.
! Scale back the thickness to account for the effects of viscosity and the fractional open
! thickness to give an appropriate non-normalized weight for each layer in determining the
! barotropic acceleration.
!$OMP parallel do default(shared)
do k=1,nz ; do J=jsh-1,jeh ; do i=ish,ieh
h_v(i,J,k) = h_v(i,J,k) * visc_rem_v(i,J,k) !### * por_face_areaV(i,J,k)
h_v(i,J,k) = h_v(i,J,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k))
enddo ; enddo ; enddo
endif

Expand All @@ -1514,7 +1523,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL,
if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
if (present(visc_rem_v)) then ; do k=1,nz
do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied
h_v(i,J,k) = h(i,j,k) * visc_rem_v(i,J,k) !### * por_face_areaV(i,J,k)
h_v(i,J,k) = h(i,j,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k))
enddo
enddo ; else ; do k=1,nz
do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied
Expand All @@ -1524,7 +1533,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL,
else
if (present(visc_rem_v)) then ; do k=1,nz
do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied
h_v(i,J,k) = h(i,j+1,k) * visc_rem_v(i,J,k) !### * por_face_areaV(i,J,k)
h_v(i,J,k) = h(i,j+1,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k))
enddo
enddo ; else ; do k=1,nz
do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied
Expand Down
40 changes: 24 additions & 16 deletions src/core/MOM_variables.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -193,13 +193,15 @@ module MOM_variables
real, pointer :: rv_x_v(:,:,:) => NULL() !< rv_x_v = rv * v at u [L T-2 ~> m s-2]
real, pointer :: rv_x_u(:,:,:) => NULL() !< rv_x_u = rv * u at v [L T-2 ~> m s-2]

real, pointer :: diag_hfrac_u(:,:,:) => NULL() !< Fractional layer thickness at u points
real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points
real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points
real, pointer :: diag_hfrac_u(:,:,:) => NULL() !< Fractional layer thickness at u points [nondim]
real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points [nondim]
real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points, modulated by the viscous
!! remnant and fractional open areas [H ~> m or kg m-2]
real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points, modulated by the viscous
!! remnant and fractional open areas [H ~> m or kg m-2]

real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points [nondim]
real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points [nondim]

end type accel_diag_ptrs

Expand Down Expand Up @@ -283,10 +285,10 @@ module MOM_variables
!! drawing from nearby to the west [H L ~> m2 or kg m-1].
real, allocatable :: FA_u_WW(:,:) !< The effective open face area for zonal barotropic transport
!! drawing from locations far to the west [H L ~> m2 or kg m-1].
real, allocatable :: uBT_WW(:,:) !< uBT_WW is the barotropic velocity [L T-1 ~> m s-1], beyond which the marginal
!! open face area is FA_u_WW. uBT_WW must be non-negative.
real, allocatable :: uBT_EE(:,:) !< uBT_EE is a barotropic velocity [L T-1 ~> m s-1], beyond which the marginal
!! open face area is FA_u_EE. uBT_EE must be non-positive.
real, allocatable :: uBT_WW(:,:) !< uBT_WW is the barotropic velocity [L T-1 ~> m s-1], beyond which the
!! marginal open face area is FA_u_WW. uBT_WW must be non-negative.
real, allocatable :: uBT_EE(:,:) !< uBT_EE is a barotropic velocity [L T-1 ~> m s-1], beyond which the
!! marginal open face area is FA_u_EE. uBT_EE must be non-positive.
real, allocatable :: FA_v_NN(:,:) !< The effective open face area for meridional barotropic transport
!! drawing from locations far to the north [H L ~> m2 or kg m-1].
real, allocatable :: FA_v_N0(:,:) !< The effective open face area for meridional barotropic transport
Expand All @@ -295,12 +297,18 @@ module MOM_variables
!! drawing from nearby to the south [H L ~> m2 or kg m-1].
real, allocatable :: FA_v_SS(:,:) !< The effective open face area for meridional barotropic transport
!! drawing from locations far to the south [H L ~> m2 or kg m-1].
real, allocatable :: vBT_SS(:,:) !< vBT_SS is the barotropic velocity, [L T-1 ~> m s-1], beyond which the marginal
!! open face area is FA_v_SS. vBT_SS must be non-negative.
real, allocatable :: vBT_NN(:,:) !< vBT_NN is the barotropic velocity, [L T-1 ~> m s-1], beyond which the marginal
!! open face area is FA_v_NN. vBT_NN must be non-positive.
real, allocatable :: h_u(:,:,:) !< An effective thickness at zonal faces [H ~> m or kg m-2].
real, allocatable :: h_v(:,:,:) !< An effective thickness at meridional faces [H ~> m or kg m-2].
real, allocatable :: vBT_SS(:,:) !< vBT_SS is the barotropic velocity, [L T-1 ~> m s-1], beyond which the
!! marginal open face area is FA_v_SS. vBT_SS must be non-negative.
real, allocatable :: vBT_NN(:,:) !< vBT_NN is the barotropic velocity, [L T-1 ~> m s-1], beyond which the
!! marginal open face area is FA_v_NN. vBT_NN must be non-positive.
real, allocatable :: h_u(:,:,:) !< An effective thickness at zonal faces, taking into account the effects
!! of vertical viscosity and fractional open areas [H ~> m or kg m-2].
!! This is primarily used as a non-normalized weight in determining
!! the depth averaged accelerations for the barotropic solver.
real, allocatable :: h_v(:,:,:) !< An effective thickness at meridional faces, taking into account the effects
!! of vertical viscosity and fractional open areas [H ~> m or kg m-2].
!! This is primarily used as a non-normalized weight in determining
!! the depth averaged accelerations for the barotropic solver.
type(group_pass_type) :: pass_polarity_BT !< Structure for polarity group halo updates
type(group_pass_type) :: pass_FA_uv !< Structure for face area group halo updates
end type BT_cont_type
Expand Down

0 comments on commit 2b2214d

Please sign in to comment.