diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 301969ed50..fe170563a4 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2362,7 +2362,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif call tracer_advect_init(Time, G, param_file, diag, CS%tracer_adv_CSp) - call tracer_hor_diff_init(Time, G, param_file, diag, CS%tv%eqn_of_state, & + call tracer_hor_diff_init(Time, G, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, & CS%tracer_diff_CSp) call lock_tracer_registry(CS%tracer_Reg) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 25d4eadb7d..1c2e23c9d8 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2397,19 +2397,23 @@ end subroutine legacy_diabatic !> Returns pointers or values of members within the diabatic_CS type. For extensibility, !! each returned argument is an optional argument subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, & - evap_CFL_limit, minimum_forcing_depth) - type(diabatic_CS), intent(in ) :: CS !< module control structure + evap_CFL_limit, minimum_forcing_depth, KPP_CSp, energetic_PBL_CSp) + type(diabatic_CS), intent(in ) :: CS !< module control structure ! All output arguments are optional type(opacity_CS), optional, pointer :: opacity_CSp !< A pointer to be set to the opacity control structure type(optics_type), optional, pointer :: optics_CSp !< A pointer to be set to the optics control structure + type(KPP_CS), optional, pointer :: KPP_CSp !< A pointer to be set to the KPP CS + type(energetic_PBL_CS), optional, pointer :: energetic_PBL_CSp !< A pointer to be set to the ePBL CS real, optional, intent( out) :: evap_CFL_limit ! CS%opacity_CSp - if (present(optics_CSp)) optics_CSp => CS%optics + if (present(opacity_CSp)) opacity_CSp => CS%opacity_CSp + if (present(optics_CSp)) optics_CSp => CS%optics + if (present(KPP_CSp)) KPP_CSp => CS%KPP_CSp + if (present(energetic_PBL_CSp)) energetic_PBL_CSp => CS%energetic_PBL_CSp ! Constants within diabatic_CS if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit diff --git a/src/tracer/MOM_lateral_boundary_mixing.F90 b/src/tracer/MOM_lateral_boundary_mixing.F90 new file mode 100644 index 0000000000..4e4cc9f455 --- /dev/null +++ b/src/tracer/MOM_lateral_boundary_mixing.F90 @@ -0,0 +1,667 @@ +!> Calculate and apply diffusive fluxes as a parameterization of lateral mixing (non-neutral) by +!! mesoscale eddies near the top and bottom boundary layers of the ocean. +module MOM_lateral_boundary_mixing + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_diag_mediator, only : diag_ctrl, time_type +use MOM_diag_mediator, only : post_data, register_diag_field +use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe +use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_grid, only : ocean_grid_type +use MOM_remapping, only : remapping_CS, initialize_remapping +use MOM_remapping, only : extract_member_remapping_CS, build_reconstructions_1d +use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme +use MOM_tracer_registry, only : tracer_registry_type, tracer_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member + +implicit none ; private + +public near_boundary_unit_tests, lateral_boundary_mixing, lateral_boundary_mixing_init + +! Private parameters to avoid doing string comparisons for bottom or top boundary layer +integer, parameter :: SURFACE = -1 !< Set a value that corresponds to the surface bopundary +integer, parameter :: BOTTOM = 1 !< Set a value that corresponds to the bottom boundary +#include + +type, public :: lateral_boundary_mixing_CS ; private + integer :: method !< Determine which of the three methods calculate + !! 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 + type(remapping_CS) :: remap_CS !< Control structure to hold remapping configuration + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get MLD + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. +end type lateral_boundary_mixing_CS + +! This include declares and sets the variable "version". +#include "version_variable.h" +character(len=40) :: mdl = "MOM_lateral_boundary_mixing" + +contains + +!> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be +!! needed for lateral boundary mixing +logical function lateral_boundary_mixing_init(Time, G, param_file, diag, diabatic_CSp, CS) + type(time_type), target, intent(in) :: Time !< Time structure + type(ocean_grid_type), intent(in) :: G !< Grid structure + 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_mixing_CS), pointer :: CS !< Lateral boundary mixing control structure + + character(len=80) :: string ! Temporary strings + logical :: boundary_extrap + + if (ASSOCIATED(CS)) then + call MOM_error(FATAL, "lateral_boundary_mixing_init called with associated control structure.") + return + endif + + ! Log this module and master switch for turning it on/off + call log_version(param_file, mdl, version, & + "This module implements lateral boundary mixing of tracers") + call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_MIXING", lateral_boundary_mixing_init, & + "If true, enables the lateral boundary mixing module.", & + default=.false.) + + if (.not. lateral_boundary_mixing_init) then + return + endif + + allocate(CS) + CS%diag => diag + call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) + call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) + + CS%surface_boundary_scheme = -1 + if ( ASSOCIATED(CS%energetic_PBL_CSp) ) CS%surface_boundary_scheme = 1 + if ( ASSOCIATED(CS%KPP_CSp) ) CS%surface_boundary_scheme = 2 + if (CS%surface_boundary_scheme < 0) then + call MOM_error(FATAL,"Lateral boundary mixing is true, but no valid boundary layer scheme was found") + endif + + ! Read all relevant parameters and write them to the model log. + call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", CS%method, & + "Determine how to apply near-boundary lateral mixing of tracers"//& + "1. Bulk layer approach"//& + "2. Along layer approach"//& + "3. Decomposition on to pressure levels", default=1) + call get_param(param_file, mdl, "LBM_BOUNDARY_EXTRAP", boundary_extrap, & + "Use boundary extrapolation in LBM code", & + default=.false.) + call get_param(param_file, mdl, "LBM_REMAPPING_SCHEME", string, & + "This sets the reconstruction scheme used "//& + "for vertical remapping for all variables. "//& + "It can be one of the following schemes: "//& + trim(remappingSchemesDoc), default=remappingDefaultScheme) + call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ) + call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) + +end function lateral_boundary_mixing_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 +subroutine lateral_boundary_mixing(G, GV, 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 + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [m2] + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [m2] + real, intent(in) :: dt !< Tracer time step * I_numitts + !! (I_numitts in tracer_hordiff) + type(tracer_registry_type), pointer :: Reg !< Tracer registry + type(lateral_boundary_mixing_CS), intent(in) :: CS !< Control structure for this module + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [m] + + + + + + + +end subroutine lateral_boundary_mixing + +!< Calculate bulk layer value of a scalar quantity as the thickness weighted average +real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, zeta_bot) + integer :: boundary !< SURFACE or BOTTOM [nondim] + integer :: nk !< Number of layers [nondim] + integer :: deg !< Degree of polynomial [nondim] + real, dimension(nk) :: h !< Layer thicknesses [m] + real :: hBLT !< Depth of the mixing layer [m] + real, dimension(nk) :: phi !< Scalar quantity + real, dimension(nk,2) :: ppoly0_E(:,:) !< Edge value of polynomial + real, dimension(nk,deg+1) :: ppoly0_coefs(:,:) !< Coefficients of polynomial + integer :: method !< Remapping scheme to use + + integer :: k_top !< Index of the first layer within the boundary + real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer + !! (0 if none, 1. if all). For the surface, this is always 0. because + !! integration starts at the surface [nondim] + integer :: k_bot !< Index of the last layer within the boundary + real :: zeta_bot !< Fraction of the layer encompassed by the surface boundary layer + !! (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 + + + htot = 0. + bulk_average = 0. + if (boundary == SURFACE) then + htot = (h(k_bot) * zeta_bot) + bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_bot, 0., zeta_bot) * htot + do k = k_bot-1,1,-1 + bulk_average = bulk_average + phi(k)*h(k) + htot = htot + h(k) + enddo + elseif (boundary == BOTTOM) then + htot = (h(k_top) * zeta_top) + ! (note 1-zeta_top because zeta_top is the fraction of the layer) + bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_top, (1.-zeta_top), 1.) * htot + do k = k_top+1,nk + bulk_average = bulk_average + phi(k)*h(k) + htot = htot + h(k) + enddo + else + call MOM_error(FATAL, "bulk_average: a valid boundary type must be provided.") + endif + + if (htot > 0.) then + bulk_average = bulk_average / hBLT + else + bulk_average = 0. + endif + +end function bulk_average + +!> Calculate the harmonic mean of two quantities +real function harmonic_mean(h1,h2) + real :: h1 !< Scalar quantity + real :: h2 !< Scalar quantity + + harmonic_mean = 2.*(h1*h2)/(h1+h2) +end function harmonic_mean + +!> Find the k-index range corresponding to the layers that are within the boundary-layer region +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, 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) + integer, intent( out) :: k_top !< Index of the first layer within the boundary + real, intent( out) :: zeta_top !< Distance from the top of a layer to the intersection of the + !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim] + integer, intent( out) :: k_bot !< Index of the last layer within the boundary + real, intent( out) :: zeta_bot !< Distance of the lower layer to the boundary layer depth + !! (0 at top, 1 at bottom) [nondim] + ! Local variables + real :: htot + integer :: k + ! Surface boundary layer + if ( boundary == SURFACE ) then + k_top = 1 + zeta_top = 0. + htot = 0. + do k=1,nk + htot = htot + h(k) + if ( htot >= hbl) then + k_bot = k + zeta_bot = 1 - (htot - hbl)/h(k) + return + endif + enddo + ! Bottom boundary layer + elseif ( boundary == BOTTOM ) then + k_bot = nk + zeta_bot = 1. + htot = 0. + do k=nk,1,-1 + htot = htot + h(k) + if (htot >= hbl) then + k_top = k + zeta_top = 1 - (htot - hbl)/h(k) + return + endif + enddo + else + call MOM_error(FATAL,"Houston, we've had a problem in boundary_k_range") + endif + +end subroutine boundary_k_range + +!> Calculate the near-boundary diffusive fluxes calculated from a 'bulk model' +subroutine layer_fluxes_bulk_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) + integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + integer, intent(in ) :: nk !< Number of layers [nondim] + integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [m] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [m] + 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] + 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 ) :: khtr_u !< Horizontal diffusivities at U-point [m^2 s^-1] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [trunit s^-1] + ! Local variables + real :: F_bulk ! Total diffusive flux across the U point [trunit s^-1] + 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] + 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) + ! [trunit 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 + + ! Calculate vertical indices containing the boundary layer + call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) + call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) + ! Calculate bulk averages of various quantities + phi_L_avg = bulk_average(boundary, nk, deg, h_L, hbl_L, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_top_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) + if ( boundary == SURFACE ) then + khtr_avg = (h_u(k_bot_u) * zeta_bot_u) * khtr_u(k_bot_u) + do k=k_bot_u-1,1,-1 + khtr_avg = khtr_avg + h_u(k) * khtr_u(k) + enddo + elseif ( boundary == BOTTOM ) then + khtr_avg = (h_u(k_top_u) * (1.-zeta_top_u)) * khtr_u(k_top_u) + do k=k_top_u+1,nk + khtr_avg = khtr_avg + h_u(k) * khtr_u(k) + enddo + endif + + khtr_avg = khtr_avg / hbl_u + + ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities + heff = harmonic_mean(hbl_L, hbl_R) + F_bulk = -(khtr_avg * heff) * (phi_R_avg - phi_L_avg) + ! Calculate the layerwise sum of the vertical effective thickness. This is different than the heff calculated + ! above, but is used as a way to decompose decompose the fluxes onto the individual layers + h_means(:) = 0. + + if (boundary == SURFACE) then + k_min = MIN(k_bot_L, k_bot_R) + + ! left hand side + if (k_bot_L == k_min) then + h_work_L = h_L(k_min) * zeta_bot_L + else + h_work_L = h_L(k_min) + endif + + ! right hand side + if (k_bot_R == k_min) then + h_work_R = h_R(k_min) * zeta_bot_R + else + h_work_R = h_R(k_min) + endif + + h_means(k_min) = harmonic_mean(h_work_L,h_work_R) + + do k=1,k_min-1 + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + endif + + if (boundary == BOTTOM) then + k_max = MAX(k_top_L, k_top_R) + + ! left hand side + if (k_top_L == k_max) then + h_work_L = h_L(k_max) * zeta_top_L + else + h_work_L = h_L(k_max) + endif + + ! right hand side + if (k_top_R == k_max) then + h_work_R = h_R(k_max) * zeta_top_R + else + h_work_R = h_R(k_max) + endif + + h_means(k_max) = harmonic_mean(h_work_L,h_work_R) + + do k=nk,k_max+1,-1 + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + endif + + inv_heff = 1./SUM(h_means) + ! Decompose the bulk flux onto the individual layers + do k=1,nk + if ( SIGN(1.,F_bulk) == SIGN(1., -(phi_R(k)-phi_L(k))) ) then + F_layer(k) = F_bulk * (h_means(k)*inv_heff) + else + F_layer(k) = 0. + endif + enddo + +end subroutine layer_fluxes_bulk_method + +!> 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 + + ! Local variables + integer, parameter :: nk = 2 ! Number of layers + integer, parameter :: deg = 1 ! Degree of reconstruction (linear here) + integer, parameter :: method = 1 ! Method used for integrating polynomials + real, dimension(nk) :: phi_L, phi_R ! Tracer values (left and right column) [ nondim m^-3 ] + real, dimension(nk) :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) + real, dimension(nk,deg+1) :: phi_pp_L, phi_pp_R ! Coefficients for the linear pseudo-reconstructions + ! [ nondim m^-3 ] + + real, dimension(nk,2) :: ppoly0_E_L, ppoly0_E_R! Polynomial edge values (left and right) [concentration] + real, dimension(nk) :: h_L, h_R ! Layer thickness (left and right) [m] + real, dimension(nk) :: khtr_u ! Horizontal diffusivities at U-point [m^2 s^-1] + real :: hbl_L, hbl_R ! Depth of the boundary layer (left and right) [m] + real :: F_bulk ! Total diffusive flux across the U point [nondim s^-1] + real, dimension(nk) :: F_layer ! Diffusive flux within each layer at U-point [nondim s^-1] + real :: h_u, hblt_u ! Thickness at the u-point [m] + real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] + real :: heff ! Harmonic mean of layer thicknesses [m] + real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] + character(len=120) :: test_name ! Title of the unit test + integer :: k_top ! Index of cell containing top of boundary + real :: zeta_top ! Nondimension position + integer :: k_bot ! Index of cell containing bottom of boundary + real :: zeta_bot ! Nondimension position + near_boundary_unit_tests = .false. + + ! Unit tests for boundary_k_range + test_name = 'Surface boundary spans the entire top cell' + h_L = (/5.,5./) + call boundary_k_range(SURFACE, nk, h_L, 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, 1., test_name, 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) + + test_name = 'Bottom boundary spans the entire bottom cell' + h_L = (/5.,5./) + call boundary_k_range(BOTTOM, nk, h_L, 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, 2, 0., 2, 1., test_name, verbose) + + test_name = 'Bottom boundary spans the entire column' + h_L = (/5.,5./) + call boundary_k_range(BOTTOM, 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) + + test_name = 'Surface boundary intersects second layer' + h_L = (/10.,10./) + call boundary_k_range(SURFACE, nk, h_L, 17.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., 2, 0.75, test_name, verbose) + + test_name = 'Surface boundary intersects first layer' + h_L = (/10.,10./) + 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 = '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) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 1., test_name, verbose) + + test_name = 'Bottom boundary intersects second layer' + h_L = (/10.,10./) + call boundary_k_range(BOTTOM, 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, 2, 0.25, 2, 1., test_name, verbose) + + ! All cases in this section have hbl which are equal to the column thicknesses + test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)' + hbl_L = 10; hbl_R = 10 + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-5.0,-5.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)' + hbl_L = 10.; hbl_R = 10. + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,1./) ; phi_R = (/0.,0./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 0.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. + ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 1. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/5.0,5.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (no gradient)' + hbl_L = 10; hbl_R = 10 + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,1./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 0. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + + test_name = 'Equal hbl and different layer thicknesses (gradient right to left)' + hbl_L = 16.; hbl_R = 16. + h_L = (/10.,6./) ; h_R = (/6.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-8.0,-8.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (diagonal tracer values)' + hbl_L = 10.; hbl_R = 10. + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,0./) ; phi_R = (/0.,1./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + + test_name = 'Different hbl and different column thicknesses (gradient from right to left)' + hbl_L = 12; hbl_R = 20 + h_L = (/6.,6./) ; h_R = (/10.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) ) + + test_name = 'Different hbl and different layer thicknesses (gradient from right to left)' + hbl_L = 12; hbl_R = 20 + h_L = (/6.,6./) ; h_R = (/10.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) ) + + ! Cases where hbl < column thickness (polynomial coefficients specified for pseudo-linear reconstruction) + + test_name = 'hbl < column thickness, hbl same, constant concentration each column' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + 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.,1./) + call layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) + + test_name = 'hbl < column thickness, hbl same, linear profile right' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/0.5,2./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 1. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 2. + khtr_u = (/1.,1./) + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + 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 layer_fluxes_bulk_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) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) +end function near_boundary_unit_tests + +!> Returns true if output of near-boundary unit tests does not match correct computed values +!! and conditionally writes results to stream +logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans) + logical, intent(in) :: verbose !< If true, write results to stdout + character(len=80), intent(in) :: test_name !< Brief description of the unit test + integer, intent(in) :: nk !< Number of layers + real, dimension(nk), intent(in) :: F_calc !< Fluxes of the unitless tracer from the algorithm [s^-1] + real, dimension(nk), intent(in) :: F_ans !< Fluxes of the unitless tracer calculated by hand [s^-1] + ! Local variables + integer :: k + integer, parameter :: stdunit = 6 + + test_layer_fluxes = .false. + do k=1,nk + if ( F_calc(k) /= F_ans(k) ) then + test_layer_fluxes = .true. + write(stdunit,*) "UNIT TEST FAILED: ", test_name + write(stdunit,10) k, F_calc(k), F_ans(k) + elseif (verbose) then + write(stdunit,10) k, F_calc(k), F_ans(k) + endif + enddo + +10 format("Layer=",i3," F_calc=",f20.16," F_ans",f20.16) +end function test_layer_fluxes + +!> Return true if output of unit tests for boundary_k_range does not match answers +logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_ans, zeta_top_ans,& + k_bot_ans, zeta_bot_ans, test_name, verbose) + integer :: k_top !< Index of cell containing top of boundary + real :: zeta_top !< Nondimension position + integer :: k_bot !< Index of cell containing bottom of boundary + real :: zeta_bot !< Nondimension position + integer :: k_top_ans !< Index of cell containing top of boundary + real :: zeta_top_ans !< Nondimension position + integer :: k_bot_ans !< Index of cell containing bottom of boundary + real :: zeta_bot_ans !< Nondimension position + character(len=80) :: test_name !< Name of the unit test + logical :: verbose !< If true always print output + + integer, parameter :: stdunit = 6 + + test_boundary_k_range = k_top .ne. k_top_ans + test_boundary_k_range = test_boundary_k_range .or. (zeta_top .ne. zeta_top_ans) + test_boundary_k_range = test_boundary_k_range .or. (k_bot .ne. k_bot_ans) + test_boundary_k_range = test_boundary_k_range .or. (zeta_bot .ne. zeta_bot_ans) + + if (test_boundary_k_range) write(stdunit,*) "UNIT TEST FAILED: ", test_name + if (test_boundary_k_range .or. verbose) then + write(stdunit,20) "k_top", k_top, "k_top_ans", k_top_ans + write(stdunit,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans + write(stdunit,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans + write(stdunit,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans + endif + + 20 format(A,"=",i3,X,A,"=",i3) + 30 format(A,"=",f20.16,X,A,"=",f20.16) + + +end function test_boundary_k_range +end module MOM_lateral_boundary_mixing diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 4eb986bacd..eb62a1e07d 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -11,6 +11,7 @@ module MOM_tracer_hor_diff use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : pass_vector use MOM_debugging, only : hchksum, uvchksum +use MOM_diabatic_driver, only : diabatic_CS use MOM_EOS, only : calculate_density, EOS_type use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery @@ -58,7 +59,11 @@ module MOM_tracer_hor_diff !! the CFL limit is not violated. logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within !! tracer_hor_diff. + logical :: use_lateral_boundary_mixing !< If true, use the lateral_boundary_mixing module from within + !! tracer_hor_diff. type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. + type(lateral_boundary_mixing_CS), pointer :: lateral_boundary_mixing_CSp => NULL() !< Control structure for lateral + !! boundary mixing. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. logical :: debug !< If true, write verbose checksums for debugging purposes. @@ -1377,11 +1382,12 @@ end subroutine tracer_epipycnal_ML_diff !> Initialize lateral tracer diffusion module -subroutine tracer_hor_diff_init(Time, G, param_file, diag, EOS, CS) +subroutine tracer_hor_diff_init(Time, G, param_file, diag, EOS, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< current model time type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control type(EOS_type), target, intent(in) :: EOS !< Equation of state CS + type(diabatic_CS), pointer, intent(in) :: diabatic_CSp !< Equation of state CS type(param_file_type), intent(in) :: param_file !< parameter file type(tracer_hor_diff_CS), pointer :: CS !< horz diffusion control structure @@ -1448,9 +1454,13 @@ subroutine tracer_hor_diff_init(Time, G, param_file, diag, EOS, CS) units="nondim", default=1.0) endif - CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, param_file, diag, EOS, CS%neutral_diffusion_CSp) + CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, param_file, diag, EOS, CS%neutral_diffusion_CSp ) if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") + CS%use_lateral_boundary_mixing = lateral_boundary_mixing_init(Time, G, param_file, diag, diabatic_CSp, & + CS%lateral_boundary_mixing_CSp) + if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & + "USE_LATERAL_BOUNDARY_MIXING and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)