From 374acc55e5046a99aa385d230178c2135a2f12a4 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 5 Oct 2021 14:48:49 -0400 Subject: [PATCH] +Add BAROTROPIC_TIDAL_SAL_BUG to fix a tide bug Added a runtime flag, BAROTROPIC_TIDAL_SAL_BUG, to fix a sign error in the tidal self-attraction and loading anomalies in the barotropic solver when tides are enabled. The default is to keep the previous bug so that answers do not change, but this default will be changed after solutions have been corrected. This commit partly addresses MOM6 issue #1496, but it should only be considered to be properly handled once the default has been changed to avoid this bug. This commit will change the MOM_parameter_doc files in cases where this bug matters, but by default all answers are bitwise identical. --- src/core/MOM_barotropic.F90 | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index abc2e228f6..f912ff3275 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -225,6 +225,9 @@ module MOM_barotropic !! pressure [nondim]. Stable values are < ~1.0. !! The default is 0.9. logical :: tides !< If true, apply tidal momentum forcing. + logical :: tidal_sal_bug !< If true, the tidal self-attraction and loading anomaly in the + !! barotropic solver has the wrong sign, replicating a long-standing + !! bug. real :: G_extra !< A nondimensional factor by which gtot is enhanced. integer :: hvel_scheme !< An integer indicating how the thicknesses at !! velocity points are calculated. Valid values are @@ -1048,7 +1051,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (CS%tides) then call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de) - dgeo_de = 1.0 + det_de + CS%G_extra + if (CS%tidal_sal_bug) then + dgeo_de = 1.0 + det_de + CS%G_extra + else + dgeo_de = (1.0 - det_de) + CS%G_extra + endif else dgeo_de = 1.0 + CS%G_extra endif @@ -2792,7 +2799,11 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) det_de = 0.0 if (CS%tides) call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de) - dgeo_de = 1.0 + max(0.0, det_de + CS%G_extra) + if (CS%tidal_sal_bug) then + dgeo_de = 1.0 + max(0.0, det_de + CS%G_extra) + else + dgeo_de = 1.0 + max(0.0, CS%G_extra - det_de) + endif if (present(pbce)) then do j=js,je ; do i=is,ie gtot_E(i,j) = 0.0 ; gtot_W(i,j) = 0.0 @@ -4329,6 +4340,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, ! a restart file to the internal representation in this run. real :: mean_SL ! The mean sea level that is used along with the bathymetry to estimate the ! geometry when LINEARIZED_BT_CORIOLIS is true or BT_NONLIN_STRESS is false [Z ~> m]. + real :: det_de ! The partial derivative due to self-attraction and loading of the reference + ! geopotential with the sea surface height when tides are enabled. + ! This is typically ~0.09 or less. real, allocatable, dimension(:,:) :: lin_drag_h type(memory_size_type) :: MS type(group_pass_type) :: pass_static_data, pass_q_D_Cor @@ -4483,6 +4497,14 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, call get_param(param_file, mdl, "TIDES", CS%tides, & "If true, apply tidal momentum forcing.", default=.false.) + det_de = 0.0 + if (CS%tides .and. associated(CS%tides_CSp)) & + call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de) + call get_param(param_file, mdl, "BAROTROPIC_TIDAL_SAL_BUG", CS%tidal_sal_bug, & + "If true, the tidal self-attraction and loading anomaly in the barotropic "//& + "solver has the wrong sign, replicating a long-standing bug with a scalar "//& + "self-attraction and loading term or the SAL term from a previous simulation.", & + default=.true., do_not_log=(det_de==0.0)) call get_param(param_file, mdl, "SADOURNY", CS%Sadourny, & "If true, the Coriolis terms are discretized with the "//& "Sadourny (1975) energy conserving scheme, otherwise "//&