Skip to content

Commit

Permalink
+(*)Fix bug when RES_SCALE_MEKE_VISC = True (#1512)
Browse files Browse the repository at this point in the history 
* +(*)Fix bug when RES_SCALE_MEKE_VISC = True

  Fix a bug that can lead to segmentation faults when RES_SCALE_MEKE_VISC is
true, as noted in MOM6 issue #1464, and add better error handling to detect
related problems.  The refactoring that is a part of this may also avoid some
problems with optimized code even when RES_SCALE_MEKE_VISC it false, as
described in MOM6 issue #1463.  In addition, logic was added so that the value
of RES_SCALE_MEKE_VISC is only logged if USE_MEKE is true.  All answers are
bitwise identical in cases that worked, including the MOM6-examples test suite,
but there are some changes in the MOM_parameter_doc files, due to an irrelevant
parameter no longer being logged.

* (*)Initialize CS%res_scale_MEKE

  Initialize the logical element res_scale_MEKE in the hor_visc_CS even if
there is no Laplacian viscosity, so that a self-consistency test does not use an
initialized value.  Also improved some comments.  All answers are bitwise
identical for all cases that successfully ran before, including in all cases in
the MOM6-examples test suite.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA authored Oct 15, 2021
1 parent 5bc0be3 commit 85afd24
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 28 deletions.
61 changes: 40 additions & 21 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -433,21 +433,29 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
backscat_subround = (1.0e-16/MEKE%backscatter_Ro_c)**(1.0/MEKE%backscatter_Ro_Pow)
endif

! Toggle whether to use a Laplacian viscosity derived from MEKE
if (associated(MEKE)) then
use_MEKE_Ku = associated(MEKE%Ku)
use_MEKE_Au = associated(MEKE%Au)
else
use_MEKE_Ku = .false. ; use_MEKE_Au = .false.
endif

rescale_Kh = .false.
if (associated(VarMix)) then
rescale_Kh = VarMix%Resoln_scaled_Kh
if (rescale_Kh .and. &
if ((rescale_Kh .or. CS%res_scale_MEKE) .and. &
(.not.associated(VarMix%Res_fn_h) .or. .not.associated(VarMix%Res_fn_q))) &
call MOM_error(FATAL, "MOM_hor_visc: VarMix%Res_fn_h and " //&
"VarMix%Res_fn_q both need to be associated with Resoln_scaled_Kh.")
call MOM_error(FATAL, "MOM_hor_visc: VarMix%Res_fn_h and VarMix%Res_fn_q "//&
"both need to be associated with Resoln_scaled_Kh or RES_SCALE_MEKE_VISC.")
elseif (CS%res_scale_MEKE) then
call MOM_error(FATAL, "MOM_hor_visc: VarMix needs to be associated if "//&
"RES_SCALE_MEKE_VISC is True.")
endif

legacy_bound = (CS%Smagorinsky_Kh .or. CS%Leith_Kh) .and. &
(CS%bound_Kh .and. .not.CS%better_bound_Kh)

! Toggle whether to use a Laplacian viscosity derived from MEKE
use_MEKE_Ku = associated(MEKE%Ku)
use_MEKE_Au = associated(MEKE%Au)

if (CS%use_GME) then
do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
boundary_mask_h(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))
Expand Down Expand Up @@ -892,8 +900,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

endif ! CS%Leith_Kh

meke_res_fn = 1.

if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
sh_xx_sq = sh_xx(i,j) * sh_xx(i,j)
Expand Down Expand Up @@ -977,13 +983,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min)
enddo ; enddo

do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j)

if (use_MEKE_Ku) &
! *Add* the MEKE contribution (might be negative)
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * meke_res_fn
enddo ; enddo
if (use_MEKE_Ku) then
! *Add* the MEKE contribution (which might be negative)
if (CS%res_scale_MEKE) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j)
enddo ; enddo
else
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j)
enddo ; enddo
endif
endif

if (CS%anisotropic) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand Down Expand Up @@ -1803,6 +1814,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
! valid parameters.
logical :: split ! If true, use the split time stepping scheme.
! If false and USE_GME = True, issue a FATAL error.
logical :: use_MEKE ! If true, the MEKE parameterization is in use.
logical :: default_2018_answers
character(len=64) :: inputdir, filename
real :: deg2rad ! Converts degrees to radians
Expand All @@ -1828,16 +1840,19 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
CS%diag => diag
! Read parameters and write them to the model log.
call log_version(param_file, mdl, version, "")
! It is not clear whether these initialization lines are needed for the
! It is not clear whether all of these initialization lines are needed for the
! cases where the corresponding parameters are not read.
CS%bound_Kh = .false. ; CS%better_bound_Kh = .false. ; CS%Smagorinsky_Kh = .false. ; CS%Leith_Kh = .false.
CS%bound_Ah = .false. ; CS%better_bound_Ah = .false. ; CS%Smagorinsky_Ah = .false. ; CS%Leith_Ah = .false.
CS%use_QG_Leith_visc = .false.
CS%bound_Coriolis = .false.
CS%Modified_Leith = .false.
CS%anisotropic = .false.
CS%dynamic_aniso = .false.
Kh = 0.0 ; Ah = 0.0
! These initialization lines are needed because they are used even in cases where they are not read.
CS%anisotropic = .false.
CS%res_scale_MEKE = .false.

! If GET_ALL_PARAMS is true, all parameters are read in all cases to enable
! parameter spelling checks.
call get_param(param_file, mdl, "GET_ALL_PARAMS", get_all, default=.false.)
Expand Down Expand Up @@ -1885,13 +1900,17 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, &
"If true, use a Leith nonlinear eddy viscosity.", &
default=.false.)
! This call duplicates one that occurs 26 lines later, and is probably unneccessary.
call get_param(param_file, mdl, "MODIFIED_LEITH", CS%Modified_Leith, &
"If true, add a term to Leith viscosity which is "//&
"proportional to the gradient of divergence.", &
default=.false.)
call get_param(param_file, mdl, "USE_MEKE", use_MEKE, &
default=.false., do_not_log=.true.)
call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", CS%res_scale_MEKE, &
"If true, the viscosity contribution from MEKE is scaled by "//&
"the resolution function.", default=.false.)
"the resolution function.", default=.false., do_not_log=.not.use_MEKE)
if (.not.use_MEKE) CS%res_scale_MEKE = .false.
if (CS%Leith_Kh .or. get_all) then
call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
"The nondimensional Laplacian Leith constant, "//&
Expand All @@ -1901,15 +1920,15 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
"If true, use QG Leith nonlinear eddy viscosity.", &
default=.false.)
if (CS%use_QG_Leith_visc .and. .not. CS%Leith_Kh) call MOM_error(FATAL, &
"MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
"MOM_hor_visc.F90, hor_visc_init:"//&
"LEITH_KH must be True when USE_QG_LEITH_VISC=True.")
endif
if (CS%Leith_Kh .or. CS%Leith_Ah .or. get_all) then
call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, &
"If true, include the beta term in the Leith nonlinear eddy viscosity.", &
default=CS%Leith_Kh)
call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, &
"If true, add a term to Leith viscosity which is \n"//&
"If true, add a term to Leith viscosity which is "//&
"proportional to the gradient of divergence.", &
default=.false.)
endif
Expand Down
23 changes: 16 additions & 7 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@ module MOM_lateral_mixing_coeffs
!> Variable mixing coefficients
type, public :: VarMix_CS
logical :: use_variable_mixing !< If true, use the variable mixing.
logical :: Resoln_scaling_used !< If true, a resolution function is used somewhere to scale
!! away one of the viscosities or diffusivities when the
!! deformation radius is well resolved.
logical :: Resoln_scaled_Kh !< If true, scale away the Laplacian viscosity
!! when the deformation radius is well resolved.
logical :: Resoln_scaled_KhTh !< If true, scale away the thickness diffusivity
Expand Down Expand Up @@ -1162,6 +1165,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
real :: grid_sp_u3, grid_sp_v3 ! Intermediate quantities for Leith metrics [L3 ~> m3]
real :: wave_speed_min ! A floor in the first mode speed below which 0 is returned [L T-1 ~> m s-1]
real :: wave_speed_tol ! The fractional tolerance for finding the wave speeds [nondim]
logical :: Resoln_scaled_MEKE_visc ! If true, the viscosity contribution from MEKE is
! scaled by the resolution function.
logical :: better_speed_est ! If true, use a more robust estimate of the first
! mode wave speed as the starting point for iterations.
! This include declares and sets the variable "version".
Expand Down Expand Up @@ -1217,6 +1222,12 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"If true, the epipycnal tracer diffusivity is scaled "//&
"away when the first baroclinic deformation radius is "//&
"well resolved.", default=.false.)
call get_param(param_file, mdl, "USE_MEKE", use_MEKE, &
default=.false., do_not_log=.true.)
call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", Resoln_scaled_MEKE_visc, &
"If true, the viscosity contribution from MEKE is scaled by "//&
"the resolution function.", default=.false., do_not_log=.true.) ! Logged elsewhere.
if (.not.use_MEKE) Resoln_scaled_MEKE_visc = .false.
call get_param(param_file, mdl, "RESOLN_USE_EBT", CS%Resoln_use_ebt, &
"If true, uses the equivalent barotropic wave speed instead "//&
"of first baroclinic wave for calculating the resolution fn.",&
Expand Down Expand Up @@ -1245,13 +1256,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_FGNV_streamfn, &
default=.false., do_not_log=.true.)
CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn
call get_param(param_file, mdl, "USE_MEKE", use_MEKE, &
default=.false., do_not_log=.true.)
CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE
! Indicate whether to calculate the Eady growth rate
CS%calculate_Eady_growth_rate = use_MEKE &
.or. (KhTr_Slope_Cff>0.) &
.or. (KhTh_Slope_Cff>0.)
CS%calculate_Eady_growth_rate = use_MEKE .or. (KhTr_Slope_Cff>0.) .or. (KhTh_Slope_Cff>0.)
call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", KhTr_passivity_coeff, &
default=0., do_not_log=.true.)
CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (KhTr_passivity_coeff>0.)
Expand Down Expand Up @@ -1383,7 +1390,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
endif

oneOrTwo = 1.0
if (CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. CS%Resoln_scaled_KhTr) then
CS%Resoln_scaling_used = CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. &
CS%Resoln_scaled_KhTr .or. Resoln_scaled_MEKE_visc
if (CS%Resoln_scaling_used) then
CS%calculate_Rd_dx = .true.
CS%calculate_res_fns = .true.
allocate(CS%Res_fn_h(isd:ied,jsd:jed), source=0.0)
Expand Down Expand Up @@ -1615,7 +1624,7 @@ subroutine VarMix_end(CS)
if (associated(CS%L2u)) deallocate(CS%L2u)
if (associated(CS%L2v)) deallocate(CS%L2v)

if (CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. CS%Resoln_scaled_KhTr) then
if (CS%Resoln_scaling_used) then
deallocate(CS%Res_fn_h)
deallocate(CS%Res_fn_q)
deallocate(CS%Res_fn_u)
Expand Down

0 comments on commit 85afd24

Please sign in to comment.