From 5aaf34b7f4791a37150dde01eac785d987030248 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 25 Sep 2019 10:29:08 -0600 Subject: [PATCH] Add flux limiter for bulk layer fluxes --- src/tracer/MOM_lateral_boundary_mixing.F90 | 60 +++++++++++++++++----- 1 file changed, 47 insertions(+), 13 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_mixing.F90 b/src/tracer/MOM_lateral_boundary_mixing.F90 index 5da3d3924d..e6d68af3a2 100644 --- a/src/tracer/MOM_lateral_boundary_mixing.F90 +++ b/src/tracer/MOM_lateral_boundary_mixing.F90 @@ -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 @@ -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 @@ -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] @@ -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 ] @@ -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] @@ -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. @@ -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 @@ -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 @@ -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/) ) @@ -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/) ) @@ -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/) ) @@ -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/) ) @@ -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/) ) @@ -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/) ) @@ -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/) ) @@ -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./) ) @@ -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./) ) @@ -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