Skip to content

Option to taper neutral diffusion #248

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
alperaltuntas merged 5 commits into from
Aug 23, 2023
Merged

Option to taper neutral diffusion #248

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
61baca8
Option to taper neutral diffusion
gustavo-marques Jul 6, 2023
1e0b773
Merge branch 'dev/ncar' into neutral_diffusion_tapering
gustavo-marques Jul 10, 2023
8234e69
Add hbd to the control structure
gustavo-marques Jul 27, 2023
53ccbc3
Fix line length
gustavo-marques Jul 27, 2023
cf29f1b
Allocate hbl
gustavo-marques Jul 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
168 changes: 145 additions & 23 deletions src/tracer/MOM_neutral_diffusion.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,17 @@ module MOM_neutral_diffusion
!! density [R L2 T-2 ~> Pa]
logical :: interior_only !< If true, only applies neutral diffusion in the ocean interior.
!! That is, the algorithm will exclude the surface and bottom boundary layers.
logical :: tapering = .false. !< If true, neutral diffusion linearly decays towards zero within a
!! transition zone defined using boundary layer depths. Only available when
!! interior_only=true.
logical :: use_unmasked_transport_bug !< If true, use an older form for the accumulation of
!! neutral-diffusion transports that were unmasked, as used prior to Jan 2018.
real, allocatable, dimension(:,:) :: hbl !< Boundary layer depth [H ~> m or kg m-2]
! Coefficients used to apply tapering from neutral to horizontal direction
real, allocatable, dimension(:) :: coeff_l !< Non-dimensional coefficient in the left column,
!! at cell interfaces
real, allocatable, dimension(:) :: coeff_r !< Non-dimensional coefficient in the right column,
!! at cell interfaces
! Positions of neutral surfaces in both the u, v directions
real, allocatable, dimension(:,:,:) :: uPoL !< Non-dimensional position with left layer uKoL-1, u-point [nondim]
real, allocatable, dimension(:,:,:) :: uPoR !< Non-dimensional position with right layer uKoR-1, u-point [nondim]
Expand Down Expand Up @@ -172,6 +181,12 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS,
"If true, only applies neutral diffusion in the ocean interior."//&
"That is, the algorithm will exclude the surface and bottom"//&
"boundary layers.", default=.false.)
if (CS%interior_only) then
call get_param(param_file, mdl, "NDIFF_TAPERING", CS%tapering, &
"If true, neutral diffusion linearly decays to zero within "//&
"a transition zone defined using boundary layer depths. "//&
"Only applicable when NDIFF_INTERIOR_ONLY=True", default=.false.)
endif
call get_param(param_file, mdl, "NDIFF_USE_UNMASKED_TRANSPORT_BUG", CS%use_unmasked_transport_bug, &
"If true, use an older form for the accumulation of neutral-diffusion "//&
"transports that were unmasked, as used prior to Jan 2018. This is not "//&
Expand Down Expand Up @@ -252,11 +267,17 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS,
endif

if (CS%interior_only) then
allocate(CS%hbl(SZI_(G),SZJ_(G)), source=0.)
call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp)
call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp)
if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then
call MOM_error(FATAL,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found")
endif

if (CS%tapering) then
allocate(CS%coeff_l(SZK_(GV)+1), source=0.)
allocate(CS%coeff_r(SZK_(GV)+1), source=0.)
endif
endif
! Store a rescaling factor for use in diagnostic messages.
CS%R_to_kg_m3 = US%R_to_kg_m3
Expand Down Expand Up @@ -316,7 +337,6 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf)
! Variables used for reconstructions
real, dimension(SZK_(GV),2) :: ppoly_r_S ! Reconstruction slopes
real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum ! Summed effective face thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth [H ~> m or kg m-2]
integer :: iMethod
real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta [R L2 T-2 ~> Pa]
real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
Expand All @@ -335,14 +355,15 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf)

! Check if hbl needs to be extracted
if (CS%interior_only) then
if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H)
if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, &
if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, CS%hbl, G, US, m_to_BLD_units=GV%m_to_H)
if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, CS%hbl, G, US, &
m_to_MLD_units=GV%m_to_H)
call pass_var(hbl,G%Domain)
call pass_var(CS%hbl,G%Domain)
! get k-indices and zeta
do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1
if (G%mask2dT(i,j) > 0.0) then
call boundary_k_range(SURFACE, G%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j))
call boundary_k_range(SURFACE, G%ke, h(i,j,:), CS%hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), &
zeta_bot(i,j))
endif
enddo; enddo
! TODO: add similar code for BOTTOM boundary layer
Expand Down Expand Up @@ -585,7 +606,6 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS)
real, dimension(SZI_(G),SZJB_(G)) :: trans_y_2d ! depth integrated diffusive tracer y-transport diagn
real, dimension(SZK_(GV)) :: dTracer ! change in tracer concentration due to ndiffusion
! [H L2 conc ~> m3 conc or kg conc]

type(tracer_type), pointer :: Tracer => NULL() ! Pointer to the current tracer

integer :: i, j, k, m, ks, nk
Expand Down Expand Up @@ -619,24 +639,53 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS)
! x-flux
do j = G%jsc,G%jec ; do I = G%isc-1,G%iec
if (G%mask2dCu(I,j)>0.) then
call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i+1,j,:), &
tracer%t(i,j,:), tracer%t(i+1,j,:), &
CS%uPoL(I,j,:), CS%uPoR(I,j,:), &
CS%uKoL(I,j,:), CS%uKoR(I,j,:), &
CS%uhEff(I,j,:), uFlx(I,j,:), &
CS%continuous_reconstruction, h_neglect, CS%remap_CS, h_neglect_edge)
if (CS%tapering) then
! compute coeff_l and coeff_r and pass them to neutral_surface_flux
call compute_tapering_coeffs(G%ke+1, CS%hbl(I,j), CS%hbl(I+1,j), CS%coeff_l(:), CS%coeff_r(:), &
h(I,j,:), h(I+1,j,:))

call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i+1,j,:), &
tracer%t(i,j,:), tracer%t(i+1,j,:), &
CS%uPoL(I,j,:), CS%uPoR(I,j,:), &
CS%uKoL(I,j,:), CS%uKoR(I,j,:), &
CS%uhEff(I,j,:), uFlx(I,j,:), &
CS%continuous_reconstruction, h_neglect, &
CS%remap_CS, h_neglect_edge, CS%coeff_l(:), CS%coeff_r(:))
else
call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i+1,j,:), &
tracer%t(i,j,:), tracer%t(i+1,j,:), &
CS%uPoL(I,j,:), CS%uPoR(I,j,:), &
CS%uKoL(I,j,:), CS%uKoR(I,j,:), &
CS%uhEff(I,j,:), uFlx(I,j,:), &
CS%continuous_reconstruction, h_neglect, CS%remap_CS, h_neglect_edge)
endif
endif
enddo ; enddo

! y-flux
do J = G%jsc-1,G%jec ; do i = G%isc,G%iec
if (G%mask2dCv(i,J)>0.) then
call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i,j+1,:), &
tracer%t(i,j,:), tracer%t(i,j+1,:), &
CS%vPoL(i,J,:), CS%vPoR(i,J,:), &
CS%vKoL(i,J,:), CS%vKoR(i,J,:), &
CS%vhEff(i,J,:), vFlx(i,J,:), &
CS%continuous_reconstruction, h_neglect, CS%remap_CS, h_neglect_edge)
if (CS%tapering) then
! compute coeff_l and coeff_r and pass them to neutral_surface_flux
call compute_tapering_coeffs(G%ke+1, CS%hbl(i,J), CS%hbl(i,J+1), CS%coeff_l(:), CS%coeff_r(:), &
h(i,J,:), h(i,J+1,:))

call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i,j+1,:), &
tracer%t(i,j,:), tracer%t(i,j+1,:), &
CS%vPoL(i,J,:), CS%vPoR(i,J,:), &
CS%vKoL(i,J,:), CS%vKoR(i,J,:), &
CS%vhEff(i,J,:), vFlx(i,J,:), &
CS%continuous_reconstruction, h_neglect, &
CS%remap_CS, h_neglect_edge, CS%coeff_l(:), CS%coeff_r(:))
else

call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i,j+1,:), &
tracer%t(i,j,:), tracer%t(i,j+1,:), &
CS%vPoL(i,J,:), CS%vPoR(i,J,:), &
CS%vKoL(i,J,:), CS%vKoR(i,J,:), &
CS%vhEff(i,J,:), vFlx(i,J,:), &
CS%continuous_reconstruction, h_neglect, CS%remap_CS, h_neglect_edge)
endif
endif
enddo ; enddo

Expand Down Expand Up @@ -736,6 +785,62 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS)

end subroutine neutral_diffusion

!> Computes linear tapering coefficients at interfaces of the left and right columns
!! within a region defined by the boundary layer depths in the two columns.
subroutine compute_tapering_coeffs(ne, bld_l, bld_r, coeff_l, coeff_r, h_l, h_r)
integer, intent(in) :: ne !< Number of interfaces
real, intent(in) :: bld_l !< Boundary layer depth, left column [H ~> m or kg m-2]
real, intent(in) :: bld_r !< Boundary layer depth, right column [H ~> m or kg m-2]
real, dimension(ne-1), intent(in) :: h_l !< Layer thickness, left column [H ~> m or kg m-2]
real, dimension(ne-1), intent(in) :: h_r !< Layer thickness, right column [H ~> m or kg m-2]
real, dimension(ne), intent(inout) :: coeff_l !< Tapering coefficient, left column [nondim]
real, dimension(ne), intent(inout) :: coeff_r !< Tapering coefficient, right column [nondim]

! Local variables
real :: min_bld, max_bld ! Min/Max boundary layer depth in two adjacent columns
integer :: dummy1 ! dummy integer
real :: dummy2 ! dummy real
integer :: k_min_l, k_min_r, k_max_l, k_max_r ! Min/max vertical indices in two adjacent columns
real :: zeta_l, zeta_r ! dummy variables
integer :: k ! vertical index

! initialize coeffs
coeff_l(:) = 1.0
coeff_r(:) = 1.0

! Calculate vertical indices containing the boundary layer depths
max_bld = MAX(bld_l, bld_r)
min_bld = MIN(bld_l, bld_r)

! k_min
call boundary_k_range(SURFACE, ne-1, h_l, min_bld, dummy1, dummy2, k_min_l, &
zeta_l)
call boundary_k_range(SURFACE, ne-1, h_r, min_bld, dummy1, dummy2, k_min_r, &
zeta_r)

! k_max
call boundary_k_range(SURFACE, ne-1, h_l, max_bld, dummy1, dummy2, k_max_l, &
zeta_l)
call boundary_k_range(SURFACE, ne-1, h_r, max_bld, dummy1, dummy2, k_max_r, &
zeta_r)
! left
do k=1,k_min_l
coeff_l(k) = 0.0
enddo
do k=k_min_l+1,k_max_l+1
coeff_l(k) = (real(k - k_min_l) + 1.0)/(real(k_max_l - k_min_l) + 2.0)
enddo

! right
do k=1,k_min_r
coeff_r(k) = 0.0
enddo
do k=k_min_r+1,k_max_r+1
coeff_r(k) = (real(k - k_min_r) + 1.0)/(real(k_max_r - k_min_r) + 2.0)
enddo

end subroutine compute_tapering_coeffs

!> Returns interface scalar, Si, for a column of layer values, S.
subroutine interface_scalar(nk, h, S, Si, i_method, h_neglect)
integer, intent(in) :: nk !< Number of levels
Expand Down Expand Up @@ -1921,7 +2026,8 @@ end function absolute_positions

!> Returns a single column of neutral diffusion fluxes of a tracer.
subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, &
hEff, Flx, continuous, h_neglect, remap_CS, h_neglect_edge)
hEff, Flx, continuous, h_neglect, remap_CS, h_neglect_edge, &
coeff_l, coeff_r)
integer, intent(in) :: nk !< Number of levels
integer, intent(in) :: nsurf !< Number of neutral surfaces
integer, intent(in) :: deg !< Degree of polynomial reconstructions
Expand All @@ -1945,11 +2051,14 @@ subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, K
!! to create sublayers
real, optional, intent(in) :: h_neglect_edge !< A negligibly small width used for
!! edge value calculations if continuous is false [H ~> m or kg m-2]
real, dimension(nk+1), optional, intent(in) :: coeff_l !< Left-column diffusivity [L2 T-1 ~> m2 s-1]
real, dimension(nk+1), optional, intent(in) :: coeff_r !< Right-column diffusivity [L2 T-1 ~> m2 s-1]

! Local variables
integer :: k_sublayer, klb, klt, krb, krt
real :: T_right_top, T_right_bottom, T_right_layer, T_right_sub, T_right_top_int, T_right_bot_int
real :: T_left_top, T_left_bottom, T_left_layer, T_left_sub, T_left_top_int, T_left_bot_int
real :: dT_top, dT_bottom, dT_layer, dT_ave, dT_sublayer, dT_top_int, dT_bot_int
real :: dT_top, dT_bottom, dT_layer, dT_ave, dT_sublayer, dT_top_int, dT_bot_int, khtr_ave
real, dimension(nk+1) :: Til !< Left-column interface tracer (conc, e.g. degC)
real, dimension(nk+1) :: Tir !< Right-column interface tracer (conc, e.g. degC)
real, dimension(nk) :: aL_l !< Left-column left edge value of tracer (conc, e.g. degC)
Expand All @@ -1964,7 +2073,12 @@ subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, K
real, dimension(nk,deg+1) :: ppoly_r_coeffs_r
real, dimension(nk,deg+1) :: ppoly_r_S_l
real, dimension(nk,deg+1) :: ppoly_r_S_r
logical :: down_flux
logical :: down_flux, tapering

tapering = .false.
if (present(coeff_l) .and. present(coeff_r)) tapering = .true.
khtr_ave = 1.0

! Setup reconstruction edge values
if (continuous) then
call interface_scalar(nk, hl, Tl, Til, 2, h_neglect)
Expand All @@ -1987,6 +2101,14 @@ subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, K
if (hEff(k_sublayer) == 0.) then
Flx(k_sublayer) = 0.
else
if (tapering) then
klb = KoL(k_sublayer+1)
klt = KoL(k_sublayer)
krb = KoR(k_sublayer+1)
krt = KoR(k_sublayer)
! these are added in this order to preserve vertically-uniform diffusivity answers
khtr_ave = 0.25 * ((coeff_l(klb) + coeff_l(klt)) + (coeff_r(krb) + coeff_r(krt)))
endif
if (continuous) then
klb = KoL(k_sublayer+1)
T_left_bottom = ( 1. - PiL(k_sublayer+1) ) * Til(klb) + PiL(k_sublayer+1) * Til(klb+1)
Expand All @@ -2010,7 +2132,7 @@ subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, K
else
dT_ave = dT_layer
endif
Flx(k_sublayer) = dT_ave * hEff(k_sublayer)
Flx(k_sublayer) = dT_ave * hEff(k_sublayer) * khtr_ave
else ! Discontinuous reconstruction
! Calculate tracer values on left and right side of the neutral surface
call neutral_surface_T_eval(nk, nsurf, k_sublayer, KoL, PiL, Tl, Tid_l, deg, iMethod, &
Expand All @@ -2036,7 +2158,7 @@ subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, K
dT_sublayer >= 0. .and. dT_top_int >= 0. .and. &
dT_bot_int >= 0.)
if (down_flux) then
Flx(k_sublayer) = dT_sublayer * hEff(k_sublayer)
Flx(k_sublayer) = dT_sublayer * hEff(k_sublayer) * khtr_ave
else
Flx(k_sublayer) = 0.
endif
Expand Down