Skip to content

Commit

Permalink
extend find_rho_bottom to get bbl thickness/index
Browse files Browse the repository at this point in the history
  • Loading branch information
Raphael Dussin authored and Raphael Dussin committed Jul 23, 2024
1 parent 6f23840 commit 871e0d2
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 6 deletions.
56 changes: 54 additions & 2 deletions src/core/MOM_interface_heights.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -394,7 +394,7 @@ end subroutine find_col_avg_SpV

!> Determine the in situ density averaged over a specified distance from the bottom,
!! calculating it as the inverse of the mass-weighted average specific volume.
subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)
subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot, h_bot, k_bot)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
Expand All @@ -409,6 +409,8 @@ subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)
!! thermodynamic fields.
integer, intent(in) :: j !< j-index of row to work on
real, dimension(SZI_(G)), intent(out) :: Rho_bot !< Near-bottom density [R ~> kg m-3].
real, dimension(SZI_(G)), optional, intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m].
integer, dimension(SZI_(G)), optional, intent(out) :: k_bot !< Bottom boundary layer top layer index [nondim].

! Local variables
real :: hb(SZI_(G)) ! Running sum of the thickness in the bottom boundary layer [H ~> m or kg m-2]
Expand Down Expand Up @@ -441,6 +443,53 @@ subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)
do i=is,ie
rho_bot(i) = GV%Rho0
enddo

! Obtain bottom boundary layer thickness and index of top layer
do i=is,ie
hb(i) = 0.0 ; h_bot(i) = 0.0 ; k_bot(i) = nz
dz_bbl_rem(i) = G%mask2dT(i,j) * max(0.0, dz_avg(i))
do_i(i) = .true.
if (G%mask2dT(i,j) <= 0.0) then
h_bbl_frac(i) = 0.0
do_i(i) = .false.
endif
enddo

do k=nz,1,-1
do_any = .false.
do i=is,ie ; if (do_i(i)) then
if (dz(i,k) < dz_bbl_rem(i)) then
! This layer is fully within the averaging depth.
dz_bbl_rem(i) = dz_bbl_rem(i) - dz(i,k)
hb(i) = hb(i) + h(i,j,k)
k_bot(i) = k
do_any = .true.
else
if (dz(i,k) > 0.0) then
frac_in = dz_bbl_rem(i) / dz(i,k)
if (frac_in >= 0.5) k_bot(i) = k ! update bbl top index if >= 50% of layer
else
frac_in = 0.0
endif
h_bbl_frac(i) = frac_in * h(i,j,k)
dz_bbl_rem(i) = 0.0
do_i(i) = .false.
endif
endif ; enddo
if (.not.do_any) exit
enddo
do i=is,ie ; if (do_i(i)) then
! The nominal bottom boundary layer is thicker than the water column, but layer 1 is
! already included in the averages. These values are set so that the call to find
! the layer-average specific volume will behave sensibly.
h_bbl_frac(i) = 0.0
endif ; enddo

do i=is,ie
if (hb(i) + h_bbl_frac(i) < GV%H_subroundoff) h_bbl_frac(i) = GV%H_subroundoff
h_bot(i) = hb(i) + h_bbl_frac(i)
enddo

else
! Check that SpV_avg has been set.
if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, &
Expand All @@ -450,7 +499,7 @@ subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)
! specified distance, with care taken to avoid having compressibility lead to an imprint
! of the layer thicknesses on this density.
do i=is,ie
hb(i) = 0.0 ; SpV_h_bot(i) = 0.0
hb(i) = 0.0 ; SpV_h_bot(i) = 0.0 ; h_bot(i) = 0.0 ; k_bot(i) = nz
dz_bbl_rem(i) = G%mask2dT(i,j) * max(0.0, dz_avg(i))
do_i(i) = .true.
if (G%mask2dT(i,j) <= 0.0) then
Expand All @@ -470,10 +519,12 @@ subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)
SpV_h_bot(i) = SpV_h_bot(i) + h(i,j,k) * tv%SpV_avg(i,j,k)
dz_bbl_rem(i) = dz_bbl_rem(i) - dz(i,k)
hb(i) = hb(i) + h(i,j,k)
k_bot(i) = k
do_any = .true.
else
if (dz(i,k) > 0.0) then
frac_in = dz_bbl_rem(i) / dz(i,k)
if (frac_in >= 0.5) k_bot(i) = k ! update bbl top index if >= 50% of layer
else
frac_in = 0.0
endif
Expand Down Expand Up @@ -516,6 +567,7 @@ subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)
do i=is,ie
if (hb(i) + h_bbl_frac(i) < GV%H_subroundoff) h_bbl_frac(i) = GV%H_subroundoff
rho_bot(i) = G%mask2dT(i,j) * (hb(i) + h_bbl_frac(i)) / (SpV_h_bot(i) + h_bbl_frac(i)*SpV_bbl(i))
h_bot(i) = hb(i) + h_bbl_frac(i)
enddo
endif

Expand Down
12 changes: 8 additions & 4 deletions src/parameterizations/vertical/MOM_set_diffusivity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,6 +264,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
! local variables
real :: N2_bot(SZI_(G)) ! Bottom squared buoyancy frequency [T-2 ~> s-2]
real :: rho_bot(SZI_(G)) ! In situ near-bottom density [T-2 ~> s-2]
real :: h_bot(SZI_(G)) ! Bottom boundary layer thickness [H ~> m].
integer :: k_bot(SZI_(G)) ! Bottom boundary layer thickness top layer index [nondim].

type(diffusivity_diags) :: dd ! structure with arrays of available diags

Expand Down Expand Up @@ -447,12 +449,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
! parameterization of Kd.

!$OMP parallel do default(shared) private(dRho_int,N2_lay,Kd_lay_2d,Kd_int_2d,Kv_bkgnd,N2_int,dz, &
!$OMP N2_bot,rho_bot,KT_extra,KS_extra,TKE_to_Kd,maxTKE,dissip,kb) &
!$OMP N2_bot,rho_bot,h_bot,k_bot,KT_extra,KS_extra,TKE_to_Kd,maxTKE,dissip,kb) &
!$OMP if(.not. CS%use_CVMix_ddiff)
do j=js,je

! Set up variables related to the stratification.
call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, N2_lay, N2_int, N2_bot, rho_bot)
call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, N2_lay, N2_int, N2_bot, rho_bot, h_bot, k_bot)

if (associated(dd%N2_3d)) then
do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo
Expand Down Expand Up @@ -1024,7 +1026,7 @@ end subroutine find_TKE_to_Kd

!> Calculate Brunt-Vaisala frequency, N^2.
subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
N2_lay, N2_int, N2_bot, Rho_bot)
N2_lay, N2_int, N2_bot, Rho_bot, h_bot, k_bot)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
Expand All @@ -1050,6 +1052,8 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
intent(out) :: N2_lay !< The squared buoyancy frequency of the layers [T-2 ~> s-2].
real, dimension(SZI_(G)), intent(out) :: N2_bot !< The near-bottom squared buoyancy frequency [T-2 ~> s-2].
real, dimension(SZI_(G)), intent(out) :: Rho_bot !< Near-bottom density [R ~> kg m-3].
real, dimension(SZI_(G)), optional, intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m].
integer, dimension(SZI_(G)), optional, intent(out) :: k_bot !< Bottom boundary layer thickness top layer index [nondim].

! Local variables
real, dimension(SZI_(G),SZK_(GV)+1) :: &
Expand Down Expand Up @@ -1194,7 +1198,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &

! Average over the larger of the envelope of the topography or a minimal distance.
do i=is,ie ; dz_BBL_avg(i) = max(h_amp(i), CS%dz_BBL_avg_min) ; enddo
call find_rho_bottom(h, dz, pres, dz_BBL_avg, tv, j, G, GV, US, Rho_bot)
call find_rho_bottom(h, dz, pres, dz_BBL_avg, tv, j, G, GV, US, Rho_bot, h_bot, k_bot)

end subroutine find_N2

Expand Down

0 comments on commit 871e0d2

Please sign in to comment.