From 42c1a3294d567b110f7a8a9efb3f85ecff59114a Mon Sep 17 00:00:00 2001 From: "Alan J. Wallcraft" Date: Mon, 29 Jul 2024 15:43:35 +0000 Subject: [PATCH] Fix USE_CONT_THICKNESS bug Optionally correct a bug due to a missing halo exchange. This likely isn't needed when compiled for nonsymmetric memory, but the added halo exchange does no harm in such cases. The pass_vector call could probably be replaced with fill_symmetric_edges, except there is no subroutine fill_vector_symmetric_edges_3d. USE_CONT_THICKNESS is not yet widely used, so rather than preserving the old (incorrect) solutions by default, this bug is corrected by default. However, the previous answers can be recovered by setting the new runtime parameter USE_CONT_THICKNESS_BUG to true. This parameter is only used (and logged) when USE_CONT_THICKNESS set to true. By default, this commit does change answers in symmetric memory cases with USE_CONT_THICKNESS = True, and there is a new runtime parameter in such cases. --- .../lateral/MOM_hor_visc.F90 | 32 +++++++++++-------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 94afd0c858..a5f354beba 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -123,6 +123,7 @@ module MOM_hor_visc real :: min_grid_Ah !< Minimun horizontal biharmonic viscosity used to !! limit grid Reynolds number [L4 T-1 ~> m4 s-1] logical :: use_cont_thick !< If true, thickness at velocity points adopts h[uv] in BT_cont from continuity solver. + logical :: use_cont_thick_bug !< If true, retain an answer-changing bug for thickness at velocity points. type(ZB2020_CS) :: ZB2020 !< Zanna-Bolton 2020 control structure. logical :: use_ZB2020 !< If true, use Zanna-Bolton 2020 parameterization. @@ -282,9 +283,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, type(thickness_diffuse_CS), optional, intent(in) :: TD !< Thickness diffusion control structure type(accel_diag_ptrs), optional, intent(in) :: ADp !< Acceleration diagnostics real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - optional, intent(in) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2]. + optional, intent(inout) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - optional, intent(in) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2]. + optional, intent(inout) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2]. ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & @@ -477,6 +478,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, use_MEKE_Au = allocated(MEKE%Au) use_cont_huv = CS%use_cont_thick .and. present(hu_cont) .and. present(hv_cont) + if (use_cont_huv .and. .not.CS%use_cont_thick_bug) then + call pass_vector(hu_cont, hv_cont, G%domain, To_All+Scalar_Pair, halo=2) + endif rescale_Kh = .false. if (VarMix%use_variable_mixing) then @@ -709,7 +713,14 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, ! in OBCs, which are not ordinarily be necessary, and might not be necessary ! even with OBCs if the accelerations are zeroed at OBC points, in which ! case the j-loop for h_u could collapse to j=js=1,je+1. -RWH - if (CS%use_land_mask) then + if (use_cont_huv) then + do j=js-2,je+2 ; do I=Isq-1,Ieq+1 + h_u(I,j) = hu_cont(I,j,k) + enddo ; enddo + do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2 + h_v(i,J) = hv_cont(i,J,k) + enddo ; enddo + elseif (CS%use_land_mask) then do j=js-2,je+2 ; do I=is-2,Ieq+1 h_u(I,j) = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i+1,j)*h(i+1,j,k)) enddo ; enddo @@ -725,16 +736,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - ! The following should obviously be combined with the previous block if adopted. - if (use_cont_huv) then - do j=js-2,je+2 ; do I=Isq-1,Ieq+1 - h_u(I,j) = hu_cont(I,j,k) - enddo ; enddo - do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2 - h_v(i,J) = hv_cont(i,J,k) - enddo ; enddo - endif - ! Adjust contributions to shearing strain and interpolated values of ! thicknesses on open boundaries. if (apply_OBC) then ; do n=1,OBC%number_of_segments @@ -2140,6 +2141,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) call get_param(param_file, mdl, "USE_CONT_THICKNESS", CS%use_cont_thick, & "If true, use thickness at velocity points from continuity solver. This option "//& "currently only works with split mode.", default=.false.) + call get_param(param_file, mdl, "USE_CONT_THICKNESS_BUG", CS%use_cont_thick_bug, & + "If true, retain an answer-changing halo update bug when "//& + "USE_CONT_THICKNESS=True. This is not recommended.", & + default=.false., do_not_log=.not.CS%use_cont_thick) + call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, & "If true, use a Laplacian horizontal viscosity.", & default=.false.)