Skip to content

Commit

Permalink
+Fix multiple bugs when MEKE_GM_SRC_ALT is True
Browse files Browse the repository at this point in the history 
  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.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed May 5, 2024
1 parent 46fd6e6 commit 0d826e4
Showing 1 changed file with 63 additions and 12 deletions.
75 changes: 63 additions & 12 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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<h_frac<=1
real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! Height change across layers [Z ~> 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
Expand Down Expand Up @@ -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].
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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.)
Expand All @@ -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, "//&
Expand Down

0 comments on commit 0d826e4

Please sign in to comment.