Skip to content

Commit

Permalink
Add flux limiter for bulk layer fluxes
Browse files Browse the repository at this point in the history
  • Loading branch information
@gustavo-marques
gustavo-marques committed Sep 25, 2019
1 parent 9f15e4f commit 5aaf34b
Showing 1 changed file with 47 additions and 13 deletions.
60 changes: 47 additions & 13 deletions src/tracer/MOM_lateral_boundary_mixing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ subroutine lateral_boundary_mixing(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS)
if (G%mask2dCu(I,j)>0.) then
call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(i,j,:), h(i+1,j,:), hbl(i,j), hbl(i+1,j), &
tracer%t(i,j,:), tracer%t(i+1,j,:), ppoly0_coefs(i,j,:,:), ppoly0_coefs(i+1,j,:,:), ppoly0_E(i,j,:,:), &
ppoly0_E(i+1,j,:,:), remap_method, Coef_x(I,j), uFlx_bulk(I,j), uFlx(I,j,:))
ppoly0_E(i+1,j,:,:), remap_method, Coef_x(I,j), G%areaT(I,j), G%areaT(I+1,j), uFlx_bulk(I,j), uFlx(I,j,:))
endif
enddo
enddo
Expand All @@ -179,7 +179,7 @@ subroutine lateral_boundary_mixing(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS)
if (G%mask2dCv(i,J)>0.) then
call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), &
tracer%t(i,J,:), tracer%t(i,J+1,:), ppoly0_coefs(i,J,:,:), ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), &
ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx_bulk(i,J), vFlx(i,J,:))
ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), G%areaT(i,J), G%areaT(i,J+1), vFlx_bulk(i,J), vFlx(i,J,:))
endif
enddo
enddo
Expand Down Expand Up @@ -458,7 +458,7 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L,
end subroutine fluxes_layer_method

!> Calculate the near-boundary diffusive fluxes calculated from a 'bulk model'
subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, ppoly0_coefs_L, &
subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, ppoly0_coefs_L, &
ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim]
integer, intent(in ) :: nk !< Number of layers [nondim]
Expand All @@ -469,6 +469,8 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L,
!! layer (left) [m]
real, intent(in ) :: hbl_R !< Thickness of the boundary boundary
!! layer (left) [m]
real, intent(in ) :: area_L !< Area of the horizontal grid (left) [m^2]
real, intent(in ) :: area_R !< Area of the horizontal grid (right) [m^2]
real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [ nondim m^-3 ]
real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [ nondim m^-3 ]
real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [ nondim m^-3 ]
Expand All @@ -479,6 +481,8 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L,
real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2]
real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux [m^2 trunit]
real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^2 trunit]
real, optional, dimension(nk), intent( out) :: F_limit !< The amount of flux not applied due to limiter
!! F_layer(k) - F_max [m^2 trunit]
! Local variables
real, dimension(nk) :: h_means ! Calculate the layer-wise harmonic means [m]
real, dimension(nk) :: h_u ! Thickness at the u-point [m]
Expand All @@ -495,6 +499,8 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L,
real :: zeta_top_L, zeta_top_R, zeta_top_u
real :: zeta_bot_L, zeta_bot_R, zeta_bot_u
real :: h_work_L, h_work_R ! dummy variables
real :: F_max !< The maximum amount of flux that can leave a cell
logical :: limited !< True if the flux limiter was applied
if (hbl_L == 0. .or. hbl_R == 0.) then
F_bulk = 0.
F_layer(:) = 0.
Expand Down Expand Up @@ -573,8 +579,33 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L,
inv_heff = 1./SUM(h_means)
! Decompose the bulk flux onto the individual layers
do k=1,nk
! Limit the tracer flux so that the donor cell with positive concentration can't go negative
! If a tracer can go negative, it is unclear what the limiter should be. BOB ALISTAIR?!
if ( SIGN(1.,F_bulk) == SIGN(1., -(phi_R(k)-phi_L(k))) ) then
if (F_bulk < 0. .and. phi_R(k) > 0.) then
F_max = 0.25 * (area_R*(phi_R(k)*h_R(k)))
elseif (F_bulk > 0. .and. phi_L(k) > 0.) then
F_max = 0.25 * (area_L*(phi_L(k)*h_L(k)))
else ! The above quantities are always positive, so we can use F_max < -1 to see if we don't need to limit
F_max = -1.
endif
! Distribute bulk flux onto layers
F_layer(k) = F_bulk * (h_means(k)*inv_heff)
! Apply flux limiter calculated above
if (F_max > 0.) then
if (F_layer(k) > 0.) then
F_layer(k) = MIN(F_layer(k),F_max)
elseif (F_layer(k) < 0.) then
F_layer(k) = MAX(F_layer(k),F_max)
endif
endif
if (PRESENT(F_limiter)) then
if (limited) then
F_limiter(k) = F_layer(k) - F_max
else
F_limiter(k) = 0.
endif
endif
else
F_layer(k) = 0.
endif
Expand Down Expand Up @@ -611,6 +642,9 @@ logical function near_boundary_unit_tests( verbose )
real :: zeta_top ! Nondimension position
integer :: k_bot ! Index of cell containing bottom of boundary
real :: zeta_bot ! Nondimension position
real :: area_L,area_R ! Area of grid cell [m^2]
area_L = 1.; area_R = 1. ! Set to unity for all unit tests

near_boundary_unit_tests = .false.

! Unit tests for boundary_k_range
Expand Down Expand Up @@ -668,7 +702,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R, &
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-5.0,-5.0/) )

Expand All @@ -685,7 +719,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0.
ppoly0_E_R(2,1) = 0.; ppoly0_E_R(2,2) = 0.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/5.0,5.0/) )

Expand All @@ -702,7 +736,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) )

Expand All @@ -719,7 +753,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-8.0,-8.0/) )

Expand All @@ -736,7 +770,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) )

Expand All @@ -753,7 +787,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) )

Expand All @@ -770,7 +804,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) )

Expand All @@ -785,7 +819,7 @@ logical function near_boundary_unit_tests( verbose )
phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0.
phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0.
khtr_u = 1.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) )

Expand All @@ -802,7 +836,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0.
ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) )

Expand All @@ -819,7 +853,7 @@ logical function near_boundary_unit_tests( verbose )
ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0.
ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1.
ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3.
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,&
ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer)
near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-2.,-2./) )
! unit tests for layer by layer method
Expand Down

0 comments on commit 5aaf34b

Please sign in to comment.