From 79bea687698f60b06ae8d7f1ad14a81b170216be Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Thu, 12 Sep 2019 10:22:00 -0600 Subject: [PATCH] New unit tests for surface boundary fluxes - Add two unit tests for cases where the surface boundary layer intersects partly through a cell. 1. Right column same BLT, same thicknesses, flux from right to left, constant in the vertical 2. Right column same BLT, same thicknesses, flux from right to left, linear profile on right TODO: 1. Uncomment out previous unit tests 2. Update API in those test cases 3. Need to add similar unit tests for the bottom boundary --- src/tracer/MOM_boundary_lateral_mixing.F90 | 145 ++++++++++++--------- 1 file changed, 85 insertions(+), 60 deletions(-) diff --git a/src/tracer/MOM_boundary_lateral_mixing.F90 b/src/tracer/MOM_boundary_lateral_mixing.F90 index 52c3ac4823..301535dae3 100644 --- a/src/tracer/MOM_boundary_lateral_mixing.F90 +++ b/src/tracer/MOM_boundary_lateral_mixing.F90 @@ -49,11 +49,13 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe integer :: method !< Remapping scheme to use integer :: k_top !< Index of the first layer within the boundary - real :: 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] + 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 !< Distance of the lower layer to the boundary layer depth - !! (0 at top, 1 at bottom) [nondim] + 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 @@ -70,7 +72,8 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe enddo elseif (boundary == BOTTOM) then htot = (h(k_top) * zeta_top) - bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, 1.) * htot + ! (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) @@ -84,7 +87,6 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe else bulk_average = 0. endif - write(*,*)'bulk_average:', bulk_average end function bulk_average @@ -197,17 +199,23 @@ subroutine layer_fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, p 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) - khtr_avg = (h_u(k_bot_u) * zeta_bot_u) * khtr_u(k_bot_u) - do k=k_bot_u,1,-1 - khtr_avg = khtr_avg + h_u(k) * khtr_u(k) - enddo + 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. @@ -236,7 +244,6 @@ subroutine layer_fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, p enddo endif - if (boundary == BOTTOM) then k_max = MAX(k_top_L, k_top_R) @@ -298,46 +305,46 @@ logical function near_boundary_unit_tests( verbose ) integer :: method ! Polynomial method 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) +!! ! 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)' @@ -405,19 +412,37 @@ logical function near_boundary_unit_tests( verbose ) ! ! Cases where hbl < column thickness (polynomial coefficients specified for pseudo-linear reconstruction) - test_name = 'hbl < column thickness' +! 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./) +! 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. +! method = 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 = (/1.,1./) + 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) = 1.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(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) = 1; ppoly0_E_R(1,2) = 1 - ppoly0_E_R(2,1) = 1; ppoly0_E_R(2,2) = 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. method = 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)