From 0d826e43595b87d56dfface4fe4924c70461a458 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 5 May 2024 09:47:59 -0400 Subject: [PATCH] +Fix multiple bugs when MEKE_GM_SRC_ALT is True Fix several bugs when MEKE_GM_SRC_ALT is True, including corrections to two dimensional rescaling factors, and optionally fix a bug that sets a limit only on positive slopes but leaving negative slopes unlimited. This bug is corrected when the new runtime parameter MEKE_GM_SRC_ALT_SLOPE_BUG is false. Additionally there is a new runtime parameter MEKE_GM_SRC_ANSWER_DATE that specifies the use of rotationally symmetric expressions for PE_release_h when it is set to 20240601 or higher, but it should be noted that rotational symmetry also requires that MEKE_GM_SRC_ALT_SLOPE_BUG is false. Four new checksum calls were also added to verify the correctness of the calculation of MEKE%GM_src when MEKE_GM_SRC_ALT and DEBUG are true. By default, all answers are bitwise identical but there are two new runtime parameters in some MOM_parameter_doc files. --- .../lateral/MOM_thickness_diffuse.F90 | 75 ++++++++++++++++--- 1 file changed, 63 insertions(+), 12 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index c510acdf0d..178e6f76e2 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -82,7 +82,14 @@ module MOM_thickness_diffuse !! used for MEKE [H ~> m or kg m-2]. When the total depth is less !! than this, the diffusivity is scaled away. logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather - !! than the streamfunction for the GM source term. + !! than the streamfunction for the GM source term for MEKE. + integer :: MEKE_src_answer_date !< The vintage of the expressions in the GM energy conversion + !! calculation when MEKE_GM_SRC_ALT is true. Values below 20240601 + !! recover the answers from the original implementation, while higher + !! values use expressions that satisfy rotational symmetry. + logical :: MEKE_src_slope_bug !< If true, use a bug that limits the positive values, but not the + !! negative values, of the slopes used when MEKE_GM_SRC_ALT is true. + !! When this is true, it breaks rotational symmetry. logical :: use_GM_work_bug !< If true, use the incorrect sign for the !! top-level work tendency on the top layer. real :: Stanley_det_coeff !< The coefficient correlating SGS temperature variance with the mean @@ -635,12 +642,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! interface of a layer that is within a layer [nondim]. 0 m] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: & - Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below) [nondim] + Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below) [Z L-1 ~> nondim] hN2_y_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency ! at v-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2], ! used for calculating the potential energy release real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: & - Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below) [nondim] + Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below) [Z L-1 ~> nondim] hN2_x_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency ! at u-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2], ! used for calculating the potential energy release @@ -998,7 +1005,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K) - Slope_x_PE(I,j,k) = MIN(Slope,CS%slope_max) + if (CS%MEKE_src_slope_bug) then + Slope_x_PE(I,j,k) = MIN(Slope, CS%slope_max) + else + Slope_x_PE(I,j,k) = Slope + if (Slope > CS%slope_max) Slope_x_PE(I,j,k) = CS%slope_max + if (Slope < -CS%slope_max) Slope_x_PE(I,j,k) = -CS%slope_max + endif if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope ! Estimate the streamfunction at each interface [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -1313,7 +1326,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K) - Slope_y_PE(i,J,k) = MIN(Slope,CS%slope_max) + if (CS%MEKE_src_slope_bug) then + Slope_y_PE(i,J,k) = MIN(Slope, CS%slope_max) + else + Slope_y_PE(i,J,k) = Slope + if (Slope > CS%slope_max) Slope_y_PE(i,J,k) = CS%slope_max + if (Slope < -CS%slope_max) Slope_y_PE(i,J,k) = -CS%slope_max + endif if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) @@ -1550,14 +1569,33 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo ; enddo ; endif if (find_work .and. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then - do j=js,je ; do i=is,ie ; do k=nz,1,-1 - PE_release_h = -0.25 * (GV%H_to_RZ*US%L_to_Z**2) * & + if (CS%MEKE_src_answer_date >= 20240601) then + do j=js,je ; do i=is,ie ; do k=nz,1,-1 + PE_release_h = -0.25 * GV%H_to_RZ * & + ( (KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + & + Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k)) + & + (Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + & + Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k)) ) + MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h + enddo ; enddo ; enddo + else + do j=js,je ; do i=is,ie ; do k=nz,1,-1 + PE_release_h = -0.25 * GV%H_to_RZ * & (KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + & Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k) + & Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + & Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k)) - MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h - enddo ; enddo ; enddo + MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h + enddo ; enddo ; enddo + endif + if (CS%debug) then + call hchksum(MEKE%GM_src, 'MEKE%GM_src', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + call uvchksum("KH_[uv]", Kh_u, Kh_v, G%HI, scale=US%L_to_m**2*US%s_to_T, & + scalar_pair=.true.) + call uvchksum("Slope_[xy]_PE", Slope_x_PE, Slope_y_PE, G%HI, scale=US%Z_to_L) + call uvchksum("hN2_[xy]_PE", hN2_x_PE, hN2_y_PE, G%HI, scale=GV%H_to_mks*US%L_to_Z**2*US%s_to_T**2, & + scalar_pair=.true.) + endif endif ; endif if (CS%id_slope_x > 0) call post_data(CS%id_slope_x, CS%diagSlopeX, CS%diag) @@ -2224,9 +2262,25 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) + call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & + "This sets the default value for the various _ANSWER_DATE parameters.", & + default=99991231, do_not_log=.true.) + call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", CS%GM_src_alt, & "If true, use the GM energy conversion form S^2*N^2*kappa rather "//& "than the streamfunction for the GM source term.", default=.false.) + call get_param(param_file, mdl, "MEKE_GM_SRC_ANSWER_DATE", CS%MEKE_src_answer_date, & + "The vintage of the expressions in the GM energy conversion calculation when "//& + "MEKE_GM_SRC_ALT is true. Values below 20240601 recover the answers from the "//& + "original implementation, while higher values use expressions that satisfy "//& + "rotational symmetry.", & + default=20240101, do_not_log=.not.CS%GM_src_alt) ! ### Change default to default_answer_date. + call get_param(param_file, mdl, "MEKE_GM_SRC_ALT_SLOPE_BUG", CS%MEKE_src_slope_bug, & + "If true, use a bug that limits the positive values, but not the negative values, "//& + "of the slopes used when MEKE_GM_SRC_ALT is true. When this is true, it breaks "//& + "all of the symmetry rules that MOM6 is supposed to obey.", & + default=.true., do_not_log=.not.CS%GM_src_alt) ! ### Change default to False. + call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & "If true, uses the GM coefficient formulation from the GEOMETRIC "//& "framework (Marshall et al., 2012).", default=.false.) @@ -2238,9 +2292,6 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "The nondimensional coefficient governing the efficiency of the GEOMETRIC "//& "thickness diffusion.", units="nondim", default=0.05) - call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & - "This sets the default value for the various _ANSWER_DATE parameters.", & - default=99991231) call get_param(param_file, mdl, "MEKE_GEOMETRIC_ANSWER_DATE", CS%MEKE_GEOM_answer_date, & "The vintage of the expressions in the MEKE_GEOMETRIC calculation. "//& "Values below 20190101 recover the answers from the original implementation, "//&