From 4fe019156dec13b9f652b8a8732ed37a221e1371 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 9 Jan 2020 16:55:35 -0700 Subject: [PATCH 1/8] Compute tracer tendency due to lateral diffusion This commit adds the option to compute the tracer tendency due to lateral boundary diffusion. Three new diagnostics have been added: 1) Lateral diffusion tracer content tendency (*_lbd_xycont_tendency) 2) Depth integrated lateral diffusion tracer content (*_lbdxy_cont_tendency_2d) 3) Lateral diffusion tracer concentration tendency (*_lbdxy_conc_tendency) --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 40 ++++++++++++++++++ src/tracer/MOM_tracer_registry.F90 | 42 +++++++++++++++---- 2 files changed, 73 insertions(+), 9 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 07f062d1d2..6485d16c59 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -142,6 +142,8 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) real, dimension(SZI_(G),SZJB_(G)) :: vFlx_bulk !< Total calculated bulk-layer v-flux for the tracer real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport + real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency ! tendency array for diagn + real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn type(tracer_type), pointer :: Tracer => NULL() !< Pointer to the current tracer integer :: remap_method !< Reconstruction method integer :: i,j,k,m !< indices to loop over @@ -156,6 +158,12 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) do m = 1,Reg%ntr tracer => Reg%tr(m) + + ! for diagnostics + if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then + tendency(:,:,:) = 0.0 + endif + do j = G%jsc-1, G%jec+1 ! Interpolate state to interface do i = G%isc-1, G%iec+1 @@ -224,6 +232,11 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (G%mask2dT(i,j)>0.) then tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & (G%IareaT(i,j)/( h(i,j,k) + GV%H_subroundoff)) + + if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then + tendency(i,j,k) = (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) )) * G%IareaT(i,j) * Idt + endif + endif enddo ; enddo ; enddo @@ -245,6 +258,33 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) enddo; enddo; enddo call post_data(tracer%id_lbd_dfy_2d, vwork_2d, CS%diag) endif + + ! post tendency of tracer content + if (tracer%id_lbdxy_cont > 0) then + call post_data(tracer%id_lbdxy_cont, tendency(:,:,:), CS%diag) + endif + + ! post depth summed tendency for tracer content + if (tracer%id_lbdxy_cont_2d > 0) then + tendency_2d(:,:) = 0. + do j = G%jsc,G%jec ; do i = G%isc,G%iec + do k = 1, GV%ke + tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k) + enddo + enddo ; enddo + call post_data(tracer%id_lbdxy_cont_2d, tendency_2d(:,:), CS%diag) + endif + + ! post tendency of tracer concentration; this step must be + ! done after posting tracer content tendency, since we alter + ! the tendency array. + if (tracer%id_lbdxy_conc > 0) then + do k = 1, GV%ke ; do j = G%jsc,G%jec ; do i = G%isc,G%iec + tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff ) + enddo ; enddo ; enddo + call post_data(tracer%id_lbdxy_conc, tendency, CS%diag) + endif + enddo end subroutine lateral_boundary_diffusion diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index cb3e7d13af..9229074099 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -127,6 +127,7 @@ module MOM_tracer_registry integer :: id_adx_2d = -1, id_ady_2d = -1, id_dfx_2d = -1, id_dfy_2d = -1 integer :: id_adv_xy = -1, id_adv_xy_2d = -1 integer :: id_dfxy_cont = -1, id_dfxy_cont_2d = -1, id_dfxy_conc = -1 + integer :: id_lbdxy_cont = -1, id_lbdxy_cont_2d = -1, id_lbdxy_conc = -1 integer :: id_remap_conc = -1, id_remap_cont = -1, id_remap_cont_2d = -1 integer :: id_tendency = -1, id_trxh_tendency = -1, id_trxh_tendency_2d = -1 integer :: id_tr_vardec = -1 @@ -532,37 +533,60 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) enddo ; enddo ; enddo endif - ! Lateral diffusion convergence tendencies + ! Neutral/Lateral diffusion convergence tendencies if (Tr%diag_form == 1) then Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & - diag%axesTL, Time, "Lateral or neutral diffusion tracer content tendency for "//trim(shortnm), & + diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & conv_units, conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer concentration "//& + diag%axesT1, Time, "Depth integrated neutral diffusion tracer content "//& + "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & + x_cell_method='sum', y_cell_method= 'sum') + + Tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', & + diag%axesTL, Time, "Lateral diffusion tracer content tendency for "//trim(shortnm), & + conv_units, conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + + Tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', & + diag%axesT1, Time, "Depth integrated lateral diffusion tracer content "//& "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method= 'sum') else cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//& ' expressed as '//trim(lowercase(flux_longname))//& - ' content due to parameterized mesoscale diffusion' + ' content due to parameterized mesoscale neutral diffusion' Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & - diag%axesTL, Time, "Lateral or neutral diffusion tracer concentration tendency for "//trim(shortnm), & + diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & conv_units, conversion=Tr%conv_scale*US%s_to_T, cmor_field_name = trim(Tr%cmor_tendprefix)//'pmdiff', & cmor_long_name = trim(cmor_var_lname), cmor_standard_name = trim(cmor_long_std(cmor_var_lname)), & x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.) cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& - trim(lowercase(flux_longname))//' content due to parameterized mesoscale diffusion' + trim(lowercase(flux_longname))//' content due to parameterized mesoscale neutral diffusion' Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer "//& - "concentration tendency for "//trim(shortnm), conv_units, & + diag%axesT1, Time, "Depth integrated neutral diffusion tracer "//& + "content tendency for "//trim(shortnm), conv_units, & conversion=Tr%conv_scale*US%s_to_T, cmor_field_name=trim(Tr%cmor_tendprefix)//'pmdiff_2d', & cmor_long_name=trim(cmor_var_lname), cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & x_cell_method='sum', y_cell_method='sum') + + Tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', & + diag%axesTL, Time, "Lateral diffusion tracer content tendency for "//trim(shortnm), & + conv_units, conversion=Tr%conv_scale*US%s_to_T, & + x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.) + + Tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', & + diag%axesT1, Time, "Depth integrated lateral diffusion tracer "//& + "content tendency for "//trim(shortnm), conv_units, & + conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum') endif Tr%id_dfxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_conc_tendency', & - diag%axesTL, Time, "Lateral (neutral) tracer concentration tendency for "//trim(shortnm), & + diag%axesTL, Time, "Neutral diffusion tracer concentration tendency for "//trim(shortnm), & + trim(units)//' s-1', conversion=US%s_to_T) + + Tr%id_lbdxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_conc_tendency', & + diag%axesTL, Time, "Lateral diffusion tracer concentration tendency for "//trim(shortnm), & trim(units)//' s-1', conversion=US%s_to_T) var_lname = "Net time tendency for "//lowercase(flux_longname) From 31d2941d0bc23973d06f45f0a596cb48dc3cb79f Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Fri, 10 Jan 2020 00:59:53 +0000 Subject: [PATCH 2/8] Fix bug in boundary_k_range if hbl > htot In cases where the boundary layer depth is larger than the column thickness, the returned indices of the layer would point to the top of the column. This behavior is fixed such that if hbl > htot, k_bot and z_bot point to the bottom of the column. --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 07f062d1d2..83b981f04b 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -338,6 +338,11 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b k_bot = 1 zeta_bot = 0. if (hbl == 0.) return + if ( hbl >= htot ); then + k_bot = nk + zeta_bot = 0. + return + endif do k=1,nk htot = htot + h(k) if ( htot >= hbl) then @@ -354,10 +359,15 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b zeta_bot = 1. htot = 0. if (hbl == 0.) return + if (hbl >= htot) then + k_top = 1 + zeta_top = 0. + return + endif do k=nk,1,-1 htot = htot + h(k) if (htot >= hbl) then - k_top = k + k_top = k zeta_top = 1 - (htot - hbl)/h(k) return endif From 5c8b32fb1b0971d7820b3d237a3a2681d2ecb679 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 10 Jan 2020 13:52:33 -0700 Subject: [PATCH 3/8] Fix a bug when checking if hbl > htot --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 1df84ad130..d16ea81291 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -155,7 +155,6 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US) call pass_var(hbl,G%Domain) - do m = 1,Reg%ntr tracer => Reg%tr(m) @@ -378,7 +377,7 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b k_bot = 1 zeta_bot = 0. if (hbl == 0.) return - if ( hbl >= htot ); then + if (hbl >= SUM(h(:))) then k_bot = nk zeta_bot = 0. return @@ -399,7 +398,7 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b zeta_bot = 1. htot = 0. if (hbl == 0.) return - if (hbl >= htot) then + if (hbl >= SUM(h(:))) then k_top = 1 zeta_top = 0. return From dce59f430dfc1cc9bc2d9524af0993f843c42c90 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 15 Jan 2020 15:19:18 -0700 Subject: [PATCH 4/8] Fix a bug in the LBD method 2 khtr_u was missing in the F_layer calculations for the surface in method 2. --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index d16ea81291..edd7cf597c 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -209,9 +209,9 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) do j=G%jsc,G%jec do i=G%isc-1,G%iec if (G%mask2dCu(I,j)>0.) then - call fluxes_layer_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(I,j,:)) + call fluxes_layer_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(I,j,:)) endif enddo enddo @@ -356,7 +356,7 @@ end function harmonic_mean subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot) integer, intent(in ) :: boundary !< SURFACE or BOTTOM [nondim] integer, intent(in ) :: nk !< Number of layers [nondim] - real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the coluymn [m] + real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the column [m] real, intent(in ) :: hbl !< Thickness of the boundary layer [m] !! If surface, with respect to zbl_ref = 0. !! If bottom, with respect to zbl_ref = SUM(h) @@ -431,7 +431,7 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, real, intent(in ) :: hbl_L !< Thickness of the boundary boundary !! layer (left) [m] real, intent(in ) :: hbl_R !< Thickness of the boundary boundary - !! layer (left) [m] + !! layer (right) [m] 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 ] @@ -487,10 +487,10 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_bot_R, 0., zeta_bot_R) heff = harmonic_mean(h_work_L, h_work_R) ! tracer flux where the minimum BLD intersets layer - F_layer(k_bot_min) = -heff * (phi_R_avg - phi_L_avg) + F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) do k = k_bot_min-1,1,-1 heff = harmonic_mean(h_L(k), h_R(k)) - F_layer(k) = -heff * (phi_R(k) - phi_L(k)) + F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) enddo endif From c8361e799817da00481f0410bcb7ee7a0fcc401d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 15 Jan 2020 15:39:30 -0700 Subject: [PATCH 5/8] Add a note saying khtr_avg should be computed once khtr is 3D --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index edd7cf597c..b4c0b9b9ac 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -446,6 +446,8 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, real, dimension(nk) :: h_u ! Thickness at the u-point [m] real :: hbl_u ! Boundary layer Thickness at the u-point [m] real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] + ! This is just to remind developers that khtr_avg should be + ! computed once khtr is 3D. real :: heff ! Harmonic mean of layer thicknesses [m] real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] real :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) @@ -487,6 +489,7 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_bot_R, 0., zeta_bot_R) heff = harmonic_mean(h_work_L, h_work_R) ! tracer flux where the minimum BLD intersets layer + ! GMM, khtr_avg should be computed once khtr is 3D F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) do k = k_bot_min-1,1,-1 heff = harmonic_mean(h_L(k), h_R(k)) @@ -556,6 +559,8 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, real, dimension(nk) :: h_u ! Thickness at the u-point [m] real :: hbl_u ! Boundary layer Thickness at the u-point [m] real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] + ! This is just to remind developers that khtr_avg should be + ! computed once khtr is 3D. real :: heff ! Harmonic mean of layer thicknesses [m] real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] real :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) @@ -593,6 +598,7 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, call boundary_k_range(boundary, nk, h_u, hbl_u, k_top_u, zeta_top_u, k_bot_u, zeta_bot_u) ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities + ! GMM, khtr_avg should be computed once khtr is 3D heff = harmonic_mean(hbl_L, hbl_R) F_bulk = -(khtr_u * heff) * (phi_R_avg - phi_L_avg) F_bulk_remain = F_bulk From 248a87ca467f124526cc0c6680250af456e58275 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 22 Jan 2020 16:15:49 -0700 Subject: [PATCH 6/8] Fix units description and delete placeholder for the pressure reconstruction --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 219 ++---------------- 1 file changed, 17 insertions(+), 202 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index b4c0b9b9ac..22a82ecad5 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -136,9 +136,9 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) real, dimension(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1) :: ppoly0_coefs !< Coefficients of polynomial real, dimension(SZI_(G),SZJ_(G),SZK_(G),2) :: ppoly0_E !< Edge values from reconstructions real, dimension(SZK_(G),CS%deg+1) :: ppoly_S !< Slopes from reconstruction (placeholder) - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uFlx !< Zonal flux of tracer [H conc ~> m conc or conc kg m-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uFlx !< Zonal flux of tracer [conc m^3] real, dimension(SZIB_(G),SZJ_(G)) :: uFLx_bulk !< Total calculated bulk-layer u-flux for the tracer - real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vFlx !< Meridional flux of tracer + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vFlx !< Meridional flux of tracer [conc m^3] real, dimension(SZI_(G),SZJB_(G)) :: vFlx_bulk !< Total calculated bulk-layer v-flux for the tracer real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport @@ -432,15 +432,15 @@ subroutine fluxes_layer_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 (right) [m] - 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 ] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [ nondim m^-3 ] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [ nondim ] real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [ nondim ] integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] - real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^2 conc] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point [m^3 conc] ! 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] @@ -542,18 +542,18 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, !! 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 ] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [ nondim m^-3 ] - real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [ nondim ] - real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [ nondim ] - integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] + real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] + real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] + integer, intent(in ) :: method !< Method of polynomial integration [nondim] 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 conc] - real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^2 conc] + real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux [m^3 conc] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^3 conc] real, optional, dimension(nk), intent( out) :: F_limit !< The amount of flux not applied due to limiter - !! F_layer(k) - F_max [m^2 conc] + !! F_layer(k) - F_max [m^3 conc] ! 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] @@ -709,191 +709,6 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, end subroutine fluxes_bulk_method -! TODO: GMM, this is a placeholder for the pressure reconstruction. -! get rid of all the T/S related variables below. We need to use the -! continuous version since pressure will be continuous. However, -! for tracer we will need to use a discontinuous reconstruction. -! Mimic the neutral diffusion driver to calculate and apply sub-layer -! fluxes. - -!> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S -!subroutine find_neutral_surface_positions_continuous(nk, Pl, Pr, PoL, PoR, KoL, KoR, hEff) -! integer, intent(in) :: nk !< Number of levels -! real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure [Pa] -! real, dimension(2*nk+2), intent(inout) :: PoL !< Fractional position of neutral surface within -! !! layer KoL of left column -! real, dimension(2*nk+2), intent(inout) :: PoR !< Fractional position of neutral surface within -! !! layer KoR of right column -! integer, dimension(2*nk+2), intent(inout) :: KoL !< Index of first left interface above neutral surface -! integer, dimension(2*nk+2), intent(inout) :: KoR !< Index of first right interface above neutral surface -! real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa] -! -! ! Local variables -! integer :: ns ! Number of neutral surfaces -! integer :: k_surface ! Index of neutral surface -! integer :: kl ! Index of left interface -! integer :: kr ! Index of right interface -! real :: dRdT, dRdS ! dRho/dT and dRho/dS for the neutral surface -! logical :: searching_left_column ! True if searching for the position of a right interface in the left column -! logical :: searching_right_column ! True if searching for the position of a left interface in the right column -! logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target -! integer :: krm1, klm1 -! real :: dRho, dRhoTop, dRhoBot, hL, hR -! integer :: lastK_left, lastK_right -! real :: lastP_left, lastP_right -! -! ns = 2*nk+2 -! ! Initialize variables for the search -! kr = 1 ; lastK_right = 1 ; lastP_right = 0. -! kl = 1 ; lastK_left = 1 ; lastP_left = 0. -! reached_bottom = .false. -! -! ! Loop over each neutral surface, working from top to bottom -! neutral_surfaces: do k_surface = 1, ns -! klm1 = max(kl-1, 1) -! if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!' -! krm1 = max(kr-1, 1) -! if (krm1>nk) stop 'find_neutral_surface_positions(): krm1 went out of bounds!' -! -! ! TODO: GMM, instead of dRho we need dP (pressure at right - pressure at left) -! -! ! Potential density difference, rho(kr) - rho(kl) -! dRho = 0.5 * ( ( dRdTr(kr) + dRdTl(kl) ) * ( Tr(kr) - Tl(kl) ) & -! + ( dRdSr(kr) + dRdSl(kl) ) * ( Sr(kr) - Sl(kl) ) ) -! ! Which column has the lighter surface for the current indexes, kr and kl -! if (.not. reached_bottom) then -! if (dRho < 0.) then -! searching_left_column = .true. -! searching_right_column = .false. -! elseif (dRho > 0.) then -! searching_right_column = .true. -! searching_left_column = .false. -! else ! dRho == 0. -! if (kl + kr == 2) then ! Still at surface -! searching_left_column = .true. -! searching_right_column = .false. -! else ! Not the surface so we simply change direction -! searching_left_column = .not. searching_left_column -! searching_right_column = .not. searching_right_column -! endif -! endif -! endif -! -! if (searching_left_column) then -! ! Interpolate for the neutral surface position within the left column, layer klm1 -! ! Potential density difference, rho(kl-1) - rho(kr) (should be negative) -! dRhoTop = 0.5 * ( ( dRdTl(klm1) + dRdTr(kr) ) * ( Tl(klm1) - Tr(kr) ) & -! + ( dRdSl(klm1) + dRdSr(kr) ) * ( Sl(klm1) - Sr(kr) ) ) -! ! Potential density difference, rho(kl) - rho(kr) (will be positive) -! dRhoBot = 0.5 * ( ( dRdTl(klm1+1) + dRdTr(kr) ) * ( Tl(klm1+1) - Tr(kr) ) & -! + ( dRdSl(klm1+1) + dRdSr(kr) ) * ( Sl(klm1+1) - Sr(kr) ) ) -! -! ! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1 -! ! unless we are still at the top of the left column (kl=1) -! if (dRhoTop > 0. .or. kr+kl==2) then -! PoL(k_surface) = 0. ! The right surface is lighter than anything in layer klm1 -! elseif (dRhoTop >= dRhoBot) then ! Left layer is unstratified -! PoL(k_surface) = 1. -! else -! ! Linearly interpolate for the position between Pl(kl-1) and Pl(kl) where the density difference -! ! between right and left is zero. -! -! ! TODO: GMM, write the linear solution instead of using interpolate_for_nondim_position -! PoL(k_surface) = interpolate_for_nondim_position( dRhoTop, Pl(klm1), dRhoBot, Pl(klm1+1) ) -! endif -! if (PoL(k_surface)>=1. .and. klm1= is really ==, when PoL==1 we point to the bottom of the cell -! klm1 = klm1 + 1 -! PoL(k_surface) = PoL(k_surface) - 1. -! endif -! if (real(klm1-lastK_left)+(PoL(k_surface)-lastP_left)<0.) then -! PoL(k_surface) = lastP_left -! klm1 = lastK_left -! endif -! KoL(k_surface) = klm1 -! if (kr <= nk) then -! PoR(k_surface) = 0. -! KoR(k_surface) = kr -! else -! PoR(k_surface) = 1. -! KoR(k_surface) = nk -! endif -! if (kr <= nk) then -! kr = kr + 1 -! else -! reached_bottom = .true. -! searching_right_column = .true. -! searching_left_column = .false. -! endif -! elseif (searching_right_column) then -! ! Interpolate for the neutral surface position within the right column, layer krm1 -! ! Potential density difference, rho(kr-1) - rho(kl) (should be negative) -! dRhoTop = 0.5 * ( ( dRdTr(krm1) + dRdTl(kl) ) * ( Tr(krm1) - Tl(kl) ) & -! + ( dRdSr(krm1) + dRdSl(kl) ) * ( Sr(krm1) - Sl(kl) ) ) -! ! Potential density difference, rho(kr) - rho(kl) (will be positive) -! dRhoBot = 0.5 * ( ( dRdTr(krm1+1) + dRdTl(kl) ) * ( Tr(krm1+1) - Tl(kl) ) & -! + ( dRdSr(krm1+1) + dRdSl(kl) ) * ( Sr(krm1+1) - Sl(kl) ) ) -! -! ! Because we are looking right, the left surface, kl, is lighter than krm1+1 and should be denser than krm1 -! ! unless we are still at the top of the right column (kr=1) -! if (dRhoTop >= 0. .or. kr+kl==2) then -! PoR(k_surface) = 0. ! The left surface is lighter than anything in layer krm1 -! elseif (dRhoTop >= dRhoBot) then ! Right layer is unstratified -! PoR(k_surface) = 1. -! else -! ! Linearly interpolate for the position between Pr(kr-1) and Pr(kr) where the density difference -! ! between right and left is zero. -! PoR(k_surface) = interpolate_for_nondim_position( dRhoTop, Pr(krm1), dRhoBot, Pr(krm1+1) ) -! endif -! if (PoR(k_surface)>=1. .and. krm1= is really ==, when PoR==1 we point to the bottom of the cell -! krm1 = krm1 + 1 -! PoR(k_surface) = PoR(k_surface) - 1. -! endif -! if (real(krm1-lastK_right)+(PoR(k_surface)-lastP_right)<0.) then -! PoR(k_surface) = lastP_right -! krm1 = lastK_right -! endif -! KoR(k_surface) = krm1 -! if (kl <= nk) then -! PoL(k_surface) = 0. -! KoL(k_surface) = kl -! else -! PoL(k_surface) = 1. -! KoL(k_surface) = nk -! endif -! if (kl <= nk) then -! kl = kl + 1 -! else -! reached_bottom = .true. -! searching_right_column = .false. -! searching_left_column = .true. -! endif -! else -! stop 'Else what?' -! endif -! -! lastK_left = KoL(k_surface) ; lastP_left = PoL(k_surface) -! lastK_right = KoR(k_surface) ; lastP_right = PoR(k_surface) -! -! ! Effective thickness -! ! NOTE: This would be better expressed in terms of the layers thicknesses rather -! ! than as differences of position - AJA -! -! ! TODO: GMM, we need to import absolute_position from neutral diffusion. This gives us -! !! the depth of the interface on the left and right side. -! -! if (k_surface>1) then -! hL = absolute_position(nk,ns,Pl,KoL,PoL,k_surface) - absolute_position(nk,ns,Pl,KoL,PoL,k_surface-1) -! hR = absolute_position(nk,ns,Pr,KoR,PoR,k_surface) - absolute_position(nk,ns,Pr,KoR,PoR,k_surface-1) -! if ( hL + hR > 0.) then -! hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) ! Harmonic mean of layer thicknesses -! else -! hEff(k_surface-1) = 0. -! endif -! endif -! -! enddo neutral_surfaces -!end subroutine find_neutral_surface_positions_continuous - !> Unit tests for near-boundary horizontal mixing logical function near_boundary_unit_tests( verbose ) logical, intent(in) :: verbose !< If true, output additional information for debugging unit tests From 0525a4c106df3d32503610a9e6413594c96a469e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 22 Jan 2020 17:42:11 -0700 Subject: [PATCH 7/8] Change flux limiting calculation Previously, F_max was calculated based on the sign of F_bulk, F_layer and phi_*, as follows: F_max = 0.25 * (area_R*(phi_R(k)*h_R(k))) or F_max = 0.25 * (area_L*(phi_L(k)*h_L(k))), This is only based on the concentration at the donor cell and can be problematic (i.e., create new extrema). In addition, this limitor was not being applied in the layer by layer method. This commit adds the following limitor to both methods: F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) In this case, F_max is based on the tracer *gradient* and, therefore, should not create new extrema. The 0.2 comes from the following: Imagine you have a tracer extrema in the center of the domain at time = 0: t=0 0 0 1 0 0 If diffusion acts on this tracer in all directions (EWNS), the final result should look like the following: t=inf .2 .2.2.2 .2 --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 110 ++++++++++++------ 1 file changed, 75 insertions(+), 35 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 22a82ecad5..48d813faa3 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -202,16 +202,14 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uFlx_bulk*Idt, CS%diag) if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vFlx_bulk*Idt, CS%diag) - ! TODO: this is where we would filter vFlx and uFlux to get rid of checkerboard noise - ! Method #2 elseif (CS%method == 2) then do j=G%jsc,G%jec do i=G%isc-1,G%iec if (G%mask2dCu(I,j)>0.) then call fluxes_layer_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(I,j,:)) + G%areaT(I,j), G%areaT(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(I,j,:)) endif enddo enddo @@ -219,8 +217,8 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) do i=G%isc,G%iec if (G%mask2dCv(i,J)>0.) then call fluxes_layer_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(i,J,:)) + G%areaT(i,J), G%areaT(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(i,J,:)) endif enddo enddo @@ -420,8 +418,8 @@ end subroutine boundary_k_range !> Calculate the lateral boundary diffusive fluxes using the layer by layer method. !! See \ref section_method2 -subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, ppoly0_coefs_L, & - ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) +subroutine fluxes_layer_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_layer) integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] integer, intent(in ) :: nk !< Number of layers [nondim] @@ -432,6 +430,8 @@ subroutine fluxes_layer_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 (right) [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) [conc] real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] @@ -452,7 +452,8 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] real :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) ! [conc m^-3 ] - real :: htot ! Total column thickness [m] + real :: htot ! Total column thickness [m] + real :: F_max ! The maximum amount of flux that can leave a cell integer :: k, k_bot_min, k_top_max integer :: k_top_L, k_bot_L, k_top_u integer :: k_top_R, k_bot_R, k_bot_u @@ -491,9 +492,32 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, ! tracer flux where the minimum BLD intersets layer ! GMM, khtr_avg should be computed once khtr is 3D F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) + + ! limit the flux to 0.2 of the tracer *gradient* + ! Why 0.2? + ! t=0 t=inf + ! 0 .2 + ! 0 1 0 .2.2.2 + ! 0 .2 + + F_max = -0.2 * ((area_R*(phi_R_avg*h_work_R))-(area_L*(phi_L_avg*h_work_L))) + ! Apply flux limiter calculated above + if (F_max >= 0.) then + F_layer(k_bot_min) = MIN(F_layer(k_bot_min),F_max) + else + F_layer(k_bot_min) = MAX(F_layer(k_bot_min),F_max) + endif + do k = k_bot_min-1,1,-1 heff = harmonic_mean(h_L(k), h_R(k)) F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) + ! Apply flux limiter calculated above + if (F_max >= 0.) then + F_layer(k) = MIN(F_layer(k),F_max) + else + F_layer(k) = MAX(F_layer(k),F_max) + endif enddo endif @@ -518,9 +542,25 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, ! tracer flux where the minimum BLD intersets layer F_layer(k_top_max) = (-heff * khtr_u) * (phi_R_avg - phi_L_avg) + + F_max = -0.2 * ((area_R*(phi_R_avg*h_work_R))-(area_L*(phi_L_avg*h_work_L))) + ! Apply flux limiter calculated above + if (F_max >= 0.) then + F_layer(k_top_max) = MIN(F_layer(k_top_max),F_max) + else + F_layer(k_top_max) = MAX(F_layer(k_top_max),F_max) + endif + do k = k_top_max+1,nk heff = harmonic_mean(h_L(k), h_R(k)) F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) + ! Apply flux limiter calculated above + if (F_max >= 0.) then + F_layer(k) = MIN(F_layer(k),F_max) + else + F_layer(k) = MAX(F_layer(k),F_max) + endif enddo endif @@ -654,42 +694,41 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, if ( SUM(h_means) == 0. ) then return + ! Decompose the bulk flux onto the individual layers else ! Initialize remaining thickness inv_heff = 1./SUM(h_means) - ! Decompose the bulk flux onto the individual layers do k=1,nk if (h_means(k) > 0.) then - ! 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?! hfrac = h_means(k)*inv_heff F_layer(k) = F_bulk * hfrac - if ( SIGN(1.,F_bulk) == SIGN(1., F_layer(k))) then - ! limit the flux to 0.25 of the total tracer in the cell - 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 + ! limit the flux to 0.2 of the tracer *gradient* + ! Why 0.2? + ! t=0 t=inf + ! 0 .2 + ! 0 1 0 .2.2.2 + ! 0 .2 + ! + F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) + + ! check if bulk flux (or F_layer) and F_max have same direction + if ( SIGN(1.,F_bulk) == SIGN(1., F_max)) then ! Distribute bulk flux onto layers if ( ((boundary == SURFACE) .and. (k == k_min)) .or. ((boundary == BOTTOM) .and. (k == nk)) ) then - F_layer(k) = F_bulk_remain + F_layer(k) = F_bulk_remain ! GMM, are not using F_bulk_remain for now. Should we keep it? endif F_bulk_remain = F_bulk_remain - F_layer(k) ! Apply flux limiter calculated above if (F_max >= 0.) then - if (F_layer(k) > 0.) then - limited = F_layer(k) > F_max - F_layer(k) = MIN(F_layer(k),F_max) - elseif (F_layer(k) < 0.) then - limited = F_layer(k) < -F_max - F_layer(k) = MAX(F_layer(k),-F_max) ! Note negative to make the sign of flux consistent - endif + limited = F_layer(k) > F_max + F_layer(k) = MIN(F_layer(k),F_max) + else + limited = F_layer(k) < F_max + F_layer(k) = MAX(F_layer(k),F_max) endif + ! GMM, again we are not using F_limit. Should we delete it? if (PRESENT(F_limit)) then if (limited) then F_limit(k) = F_layer(k) - F_max @@ -698,6 +737,7 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, endif endif else + ! do not apply a flux on this layer F_bulk_remain = F_bulk_remain - F_layer(k) F_layer(k) = 0. endif @@ -931,8 +971,8 @@ 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, 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) + call fluxes_layer_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_layer) near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-2.,-2./) ) ! unit tests for layer by layer method @@ -949,8 +989,8 @@ 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_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + call fluxes_layer_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_layer) near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,0.0/) ) test_name = 'Different hbl and different column thicknesses (linear profile right)' @@ -967,8 +1007,8 @@ logical function near_boundary_unit_tests( verbose ) ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 2. ppoly0_E_R(2,1) = 2.; ppoly0_E_R(2,2) = 4. khtr_u = 1. - call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + call fluxes_layer_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_layer) near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-3.75,0.0/) ) end function near_boundary_unit_tests From 9db5ba14a0c4b2f5ec68655b04fc3e47bc0d4f83 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 7 Feb 2020 10:21:56 -0700 Subject: [PATCH 8/8] Improve documentation and unit tests * Fix broken unit tests * Add new unit tests (Surface boundary is deeper than column thickness) * Update documentation * Delete unneeded code --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 218 +++++++++--------- 1 file changed, 108 insertions(+), 110 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 48d813faa3..4fda621abc 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -40,7 +40,6 @@ module MOM_lateral_boundary_diffusion !! and apply near boundary layer fluxes !! 1. Bulk-layer approach !! 2. Along layer - !! 3. Decomposition onto pressure levels integer :: deg !< Degree of polynomial reconstruction integer :: surface_boundary_scheme !< Which boundary layer scheme to use !! 1. ePBL; 2. KPP @@ -65,7 +64,7 @@ logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diab type(param_file_type), intent(in) :: param_file !< Parameter file structure type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(diabatic_CS), pointer :: diabatic_CSp !< KPP control structure needed to get BLD - type(lateral_boundary_diffusion_CS), pointer :: CS !< Lateral boundary mixing control structure + type(lateral_boundary_diffusion_CS), pointer :: CS !< Lateral boundary mixing control structure ! local variables character(len=80) :: string ! Temporary strings @@ -101,8 +100,7 @@ logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diab call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", CS%method, & "Determine how to apply boundary lateral diffusion of tracers: \n"//& "1. Bulk layer approach \n"//& - "2. Along layer approach \n"//& - "3. Decomposition on to pressure levels", default=1) + "2. Along layer approach", default=1) call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, & "Use boundary extrapolation in LBD code", & default=.false.) @@ -116,8 +114,11 @@ logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diab end function lateral_boundary_diffusion_init -!> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries. Two different methods -!! Method 1: Calculate fluxes from bulk layer integrated quantities +!> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries. +!! Two different methods are available: +!! Method 1: lower order representation, calculate fluxes from bulk layer integrated quantities. +!! Method 2: more straight forward, diffusion is applied layer by layer using only information +!! from neighboring cells. subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) type(ocean_grid_type), intent(inout) :: G !< Grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -142,8 +143,8 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) real, dimension(SZI_(G),SZJB_(G)) :: vFlx_bulk !< Total calculated bulk-layer v-flux for the tracer real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport - real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency ! tendency array for diagn - real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn + real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency !< tendency array for diagn + real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d !< depth integrated content tendency for diagn type(tracer_type), pointer :: Tracer => NULL() !< Pointer to the current tracer integer :: remap_method !< Reconstruction method integer :: i,j,k,m !< indices to loop over @@ -308,8 +309,8 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1. !! because integration starts at the bottom [nondim] ! Local variables - real :: htot ! Running sum of the thicknesses (top to bottom) - integer :: k ! k indice + real :: htot !< Running sum of the thicknesses (top to bottom) + integer :: k !< k indice htot = 0. @@ -404,7 +405,7 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b do k=nk,1,-1 htot = htot + h(k) if (htot >= hbl) then - k_top = k + k_top = k zeta_top = 1 - (htot - hbl)/h(k) return endif @@ -441,26 +442,26 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point [m^3 conc] + ! 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] - real :: hbl_u ! Boundary layer Thickness at the u-point [m] - real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] - ! This is just to remind developers that khtr_avg should be - ! computed once khtr is 3D. - real :: heff ! Harmonic mean of layer thicknesses [m] - real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] - real :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) - ! [conc m^-3 ] - real :: htot ! Total column thickness [m] - real :: F_max ! The maximum amount of flux that can leave a cell - integer :: k, k_bot_min, k_top_max - integer :: k_top_L, k_bot_L, k_top_u - integer :: k_top_R, k_bot_R, k_bot_u - 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 :: hbl_min ! minimum BLD (left and right) + real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [m] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + !! This is just to remind developers that khtr_avg should be + !! computed once khtr is 3D. + real :: heff !< Harmonic mean of layer thicknesses [m] + real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) + !! [conc m^-3 ] + real :: htot !< Total column thickness [m] + integer :: k, k_bot_min, k_top_max !< k-indices, min and max for top and bottom, respectively + integer :: k_top_L, k_bot_L !< k-indices left + integer :: k_top_R, k_bot_R !< k-indices right + real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary + !! layer depth [nondim] + real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the boundary + !!layer depth [nondim] + real :: h_work_L, h_work_R !< dummy variables + real :: hbl_min !< minimum BLD (left and right) [m] F_layer(:) = 0.0 if (hbl_L == 0. .or. hbl_R == 0.) then @@ -493,31 +494,9 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L ! GMM, khtr_avg should be computed once khtr is 3D F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) - ! limit the flux to 0.2 of the tracer *gradient* - ! Why 0.2? - ! t=0 t=inf - ! 0 .2 - ! 0 1 0 .2.2.2 - ! 0 .2 - - F_max = -0.2 * ((area_R*(phi_R_avg*h_work_R))-(area_L*(phi_L_avg*h_work_L))) - ! Apply flux limiter calculated above - if (F_max >= 0.) then - F_layer(k_bot_min) = MIN(F_layer(k_bot_min),F_max) - else - F_layer(k_bot_min) = MAX(F_layer(k_bot_min),F_max) - endif - do k = k_bot_min-1,1,-1 heff = harmonic_mean(h_L(k), h_R(k)) F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) - F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) - ! Apply flux limiter calculated above - if (F_max >= 0.) then - F_layer(k) = MIN(F_layer(k),F_max) - else - F_layer(k) = MAX(F_layer(k),F_max) - endif enddo endif @@ -543,24 +522,9 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L ! tracer flux where the minimum BLD intersets layer F_layer(k_top_max) = (-heff * khtr_u) * (phi_R_avg - phi_L_avg) - F_max = -0.2 * ((area_R*(phi_R_avg*h_work_R))-(area_L*(phi_L_avg*h_work_L))) - ! Apply flux limiter calculated above - if (F_max >= 0.) then - F_layer(k_top_max) = MIN(F_layer(k_top_max),F_max) - else - F_layer(k_top_max) = MAX(F_layer(k_top_max),F_max) - endif - do k = k_top_max+1,nk heff = harmonic_mean(h_L(k), h_R(k)) F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) - F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) - ! Apply flux limiter calculated above - if (F_max >= 0.) then - F_layer(k) = MIN(F_layer(k),F_max) - else - F_layer(k) = MAX(F_layer(k),F_max) - endif enddo endif @@ -595,25 +559,26 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, real, optional, dimension(nk), intent( out) :: F_limit !< The amount of flux not applied due to limiter !! F_layer(k) - F_max [m^3 conc] ! 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] - real :: hbl_u ! Boundary layer Thickness at the u-point [m] - real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] - ! This is just to remind developers that khtr_avg should be - ! computed once khtr is 3D. - real :: heff ! Harmonic mean of layer thicknesses [m] - real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] - real :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) - ! [conc m^-3 ] - real :: htot ! Total column thickness [m] - integer :: k, k_min, k_max - integer :: k_top_L, k_bot_L, k_top_u - integer :: k_top_R, k_bot_R, k_bot_u - 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 + real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [m] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + !! This is just to remind developers that khtr_avg should be + !! computed once khtr is 3D. + real :: heff !< Harmonic mean of layer thicknesses [m] + real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) + !! [conc m^-3 ] + real :: htot ! Total column thickness [m] + integer :: k, k_min, k_max !< k-indices, min and max for top and bottom, respectively + integer :: k_top_L, k_bot_L !< k-indices left + integer :: k_top_R, k_bot_R !< k-indices right + real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the + !! boundary layer [nondim] + real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the + !! boundary layer [nondim] + real :: h_work_L, h_work_R !< dummy variables + real :: F_max !< The maximum amount of flux that can leave a + !! cell [m^3 conc] + logical :: limited !< True if the flux limiter was applied real :: hfrac, F_bulk_remain if (hbl_L == 0. .or. hbl_R == 0.) then @@ -631,11 +596,6 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, zeta_top_L, k_bot_L, zeta_bot_L) phi_R_avg = bulk_average(boundary, nk, deg, h_R, hbl_R, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, & zeta_top_R, k_bot_R, zeta_bot_R) - do k=1,nk - h_u(k) = 0.5 * (h_L(k) + h_R(k)) - enddo - hbl_u = 0.5*(hbl_L + hbl_R) - call boundary_k_range(boundary, nk, h_u, hbl_u, k_top_u, zeta_top_u, k_bot_u, zeta_bot_u) ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities ! GMM, khtr_avg should be computed once khtr is 3D @@ -791,7 +751,7 @@ logical function near_boundary_unit_tests( verbose ) test_name = 'Surface boundary spans the entire column' h_L = (/5.,5./) call boundary_k_range(SURFACE, nk, h_L, 10., k_top, zeta_top, k_bot, zeta_bot) - near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0., test_name, verbose) test_name = 'Bottom boundary spans the entire bottom cell' h_L = (/5.,5./) @@ -813,6 +773,11 @@ logical function near_boundary_unit_tests( verbose ) call boundary_k_range(SURFACE, nk, h_L, 2.5, k_top, zeta_top, k_bot, zeta_bot) near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 0.25, test_name, verbose) + test_name = 'Surface boundary is deeper than column thickness' + h_L = (/10.,10./) + call boundary_k_range(SURFACE, nk, h_L, 21.0, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0., test_name, verbose) + test_name = 'Bottom boundary intersects first layer' h_L = (/10.,10./) call boundary_k_range(BOTTOM, nk, h_L, 17.5, k_top, zeta_top, k_bot, zeta_bot) @@ -1088,41 +1053,74 @@ end function test_boundary_k_range !! Boundary lateral diffusion can be applied using one of the three methods: !! !! * [Method #1: Bulk layer](@ref section_method1) (default); -!! * [Method #2: Along layer](ref section_method2); -!! * [Method #3: Decomposition on to pressure levels](@ref section_method3). +!! * [Method #2: Along layer](@ref section_method2); !! !! A brief summary of these methods is provided below. !! !! \subsection section_method1 Bulk layer approach (Method #1) !! -!! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' +!! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'.This +!! is a lower order representation (Kraus-Turner like approach) which assumes that +!! eddies are acting along well mixed layers (i.e., eddies do not know care about +!! vertical tracer gradients within the boundary layer). +!! +!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). +!! For the TOP boundary layer, these are: !! -!! Step #1: get vertical indices containing the boundary layer depth. These are !! k_top, k_bot, zeta_top, zeta_bot !! -!! Step #2: compute bulk averages (thickness weighted). phi_L and phi_R +!! Step #2: compute bulk averages (thickness weighted) tracer averages (phi_L and phi_R), +!! then calculate the bulk diffusive flux (F_{bulk}): +!! +!! \f[ F_{bulk} = -KHTR \times h_{eff} \times (\phi_R - \phi_L), \f] +!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the boundary layer depth +!! in the left and right columns (\f[ HBL_L \f] and \f[ HBL_R \f], respectively). +!! +!! Step #3: decompose F_bulk onto individual layers: +!! +!! \f[ F_{layer}(k) = F_{bulk} \times h_{frac}(k) , \f] +!! +!! where h_{frac} is !! -!! Step #3: compute a diffusive bulk flux -!! \f[ F_{bulk} = -(KHTR \times heff) \times (\phi_R - \phi_L), \f] -!! where heff is the harmonic mean of the boundary layer depth in the left and -!! right columns (\f[ HBL_L \f] and \f[ HBL_R \f], respectively). +!! \f[ h_{frac}(k) = h_u(k) \times \frac{1}{\sum(h_u)}. \f] !! -!! Step #4: limit the tracer flux so that the donor cell, with positive -!! concentration, cannot go negative. If a tracer can go negative (e.g., -!! temperature at high latitudes) it is unclear what limiter should be used. -!! (TODO: ask Bob and Alistair). +!! h_u is the [harmonic mean](@ref section_harmonic_mean) of thicknesses at each layer. +!! Special care (layer reconstruction) must be taken at k_min = min(k_botL, k_bot_R). !! -!! Step #5: decompose the bulk flux into individual layers and keep track of -!! the remaining flux. The limiter described above is also applied during -!! this step. +!! Step #4: limit the tracer flux so that 1) only down-gradient fluxes are applied, +!! and 2) the flux cannot be larger than F_max, which is defined using the tracer +!! gradient: +!! +!! \f[ F_{max} = -0.2 \times [(V_R(k) \times \phi_R(k)) - (V_L(k) \times \phi_L(k))], \f] +!! where V is the cell volume. Why 0.2? +!! t=0 t=inf +!! 0 .2 +!! 0 1 0 .2.2.2 +!! 0 .2 !! !! \subsection section_method2 Along layer approach (Method #2) !! -!! \subsection section_method3 Decomposition on to pressure levels (Method #3) +!! This is a more straight forward method where diffusion is applied layer by layer using +!! only information from neighboring cells. +!! +!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). +!! For the TOP boundary layer, these are: +!! +!! k_top, k_bot, zeta_top, zeta_bot +!! +!! Step #2: calculate the diffusive flux at each layer: !! -!! To be implemented +!! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \f] +!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness +!! in the left and right columns. Special care (layer reconstruction) must be taken at +!! k_min = min(k_botL, k_bot_R). This method does not require a limiter since KHTR +!! is already limted based on a diffusive CFL condition prior to the call of this +!! module. !! !! \subsection section_harmonic_mean Harmonic Mean !! +!! The harmonic mean (HM) betwen h1 and h2 is defined as: +!! +!! \f[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \f] !! end module MOM_lateral_boundary_diffusion