Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

+Fix multiple bugs when MEKE_GM_SRC_ALT is True #620

Merged
merged 2 commits into from
May 8, 2024
Merged

+Fix multiple bugs when MEKE_GM_SRC_ALT is True #620

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0d826e4
+Fix multiple bugs when MEKE_GM_SRC_ALT is True
Hallberg-NOAA May 5, 2024
01ca494
Merge branch 'dev/gfdl' into MEKE_GM_src_alt_bugfix
marshallward May 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
75 changes: 63 additions & 12 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Open in desktop
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
Loading