diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index a7e8194a84..3cb1ebf399 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -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 @@ -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 diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 7bc67e2fdf..95de2fd923 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -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. @@ -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. @@ -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 @@ -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 @@ -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 @@ -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. @@ -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 @@ -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 @@ -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 diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 5de7ea7319..ba5001e427 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -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 @@ -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 @@ -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