Skip to content

Commit

Permalink
Add GL90 parameterization for stacked shallow water (#268)
Browse files Browse the repository at this point in the history 
This adds a new vertical viscosity parameterization as in Greatbatch and Lamb (1990),
Ferreira & Marshall (2006) and Zhao & Vallis (2008), hereafter referred to as the GL90
vertical viscosity parameterization. This vertical viscosity scheme redistributes momentum
in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization,
but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient
nu is computed from kappa_GM via thermal wind balance, and the following relation:

nu = kappa_GM * f^2 / N^2.

The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free
boundary conditions at the top and bottom.

In the current implementation, kappa_GM is assumed either (a) constant or as (b) having an
EBT structure. A third possible formulation of nu is depth-independent:

nu = f^2 * alpha

The latter formulation would be equivalent to a kappa_GM that varies as N^2 with depth.
Currently, the GL90 parameterization is only implemented in stacked shallow water (SSW) mode,
in which case we have 1/N^2 = h/g'.

More specifically, this commit adds a new subroutine that computes the coupling coefficient
associated with GL90 via

a_cpl_gl90 = nu / h = kappa_GM * f^2 / g'

or

a_cpl_gl90 = nu / h = f^2 * alpha / h.

Further, a_cpl_gl90 is multiplied by a function (botfn), which is 0 within the GL90 bottom boundary
layer, whose depth is set by Hbbl_gl90, and 1 otherwise. This modification is necessary to avoid fluxing
momentum into vanished layers that ride over steep topography.

Finally, a_cpl_gl90 is added to a_cpl, where the latter is the coupling coefficient associated with the
remaining vertical stresses, used in the vertical viscosity solver.

More information can be found in Loose et al. (https://www.essoar.org/doi/abs/10.1002/essoar.10512867.1),
Appendix B.

* Introduce logical variable KD_GL90_USE_EBT_STRUCT

This variable is analogous to KHTH_USE_EBT_STRUCT, but is specifically for the GL90 scheme.
If the user sets KD_GL90_USE_EBT_STRUCT = True, an EBT structure will be applied to KD_GL90.
  • Loading branch information
@NoraLoose
NoraLoose authored Dec 23, 2022
1 parent 053752d commit 7869030
Show file tree
Hide file tree
Showing 5 changed files with 289 additions and 25 deletions.
6 changes: 3 additions & 3 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -567,7 +567,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
if (CS%debug) then
call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
endif
call vertvisc_coef(up, vp, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(up, vp, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
call cpu_clock_end(id_clock_vertvisc)
if (showCallTree) call callTree_wayPoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")
Expand Down Expand Up @@ -660,7 +660,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
endif
call vertvisc_coef(up, vp, h, forces, visc, dt_pred, G, GV, US, CS%vertvisc_CSp, &
CS%OBC)
CS%OBC, VarMix)
call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, &
GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)")
Expand Down Expand Up @@ -880,7 +880,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
! u <- u + dt d/dz visc d/dz u
! u_av <- u_av + dt d/dz visc d/dz u_av
call cpu_clock_begin(id_clock_vertvisc)
call vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves)
if (G%nonblocking_updates) then
Expand Down
6 changes: 3 additions & 3 deletions src/core/MOM_dynamics_unsplit.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -345,7 +345,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
call disable_averaging(CS%diag)

dt_visc = 0.5*dt ; if (CS%use_correct_dt_visc) dt_visc = dt_pred
call vertvisc_coef(up, vp, h_av, forces, visc, dt_visc, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(up, vp, h_av, forces, visc, dt_visc, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(up, vp, h_av, forces, visc, dt_visc, CS%OBC, CS%ADp, CS%CDp, &
G, GV, US, CS%vertvisc_CSp, Waves=Waves)
call cpu_clock_end(id_clock_vertvisc)
Expand Down Expand Up @@ -405,7 +405,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &

! upp <- upp + dt/2 d/dz visc d/dz upp
call cpu_clock_begin(id_clock_vertvisc)
call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(upp, vpp, hp, forces, visc, dt*0.5, CS%OBC, CS%ADp, CS%CDp, &
G, GV, US, CS%vertvisc_CSp, Waves=Waves)
call cpu_clock_end(id_clock_vertvisc)
Expand Down Expand Up @@ -489,7 +489,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &

! u <- u + dt d/dz visc d/dz u
call cpu_clock_begin(id_clock_vertvisc)
call vertvisc_coef(u, v, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(u, v, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(u, v, h_av, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, &
G, GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, Waves=Waves)
call cpu_clock_end(id_clock_vertvisc)
Expand Down
6 changes: 3 additions & 3 deletions src/core/MOM_dynamics_unsplit_RK2.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -341,7 +341,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
call set_viscous_ML(u_in, v_in, h_av, tv, forces, visc, dt_visc, G, GV, US, CS%set_visc_CSp)
call disable_averaging(CS%diag)

call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(up, vp, h_av, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, &
G, GV, US, CS%vertvisc_CSp)
call cpu_clock_end(id_clock_vertvisc)
Expand Down Expand Up @@ -392,10 +392,10 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
! up[n] <- up* + dt d/dz visc d/dz up
! u[n] <- u*[n] + dt d/dz visc d/dz u[n]
call cpu_clock_begin(id_clock_vertvisc)
call vertvisc_coef(up, vp, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(up, vp, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(up, vp, h_av, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, &
G, GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot)
call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(u_in, v_in, h_av, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp,&
G, GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot)
call cpu_clock_end(id_clock_vertvisc)
Expand Down
14 changes: 10 additions & 4 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,8 @@ module MOM_lateral_mixing_coeffs
!! of first baroclinic wave for calculating the resolution fn.
logical :: khth_use_ebt_struct !< If true, uses the equivalent barotropic structure
!! as the vertical structure of thickness diffusivity.
logical :: kdgl90_use_ebt_struct !< If true, uses the equivalent barotropic structure
!! as the vertical structure of diffusivity in the GL90 scheme.
logical :: calculate_cg1 !< If true, calls wave_speed() to calculate the first
!! baroclinic wave speed and populate CS%cg1.
!! This parameter is set depending on other parameters.
Expand Down Expand Up @@ -229,7 +231,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
if (CS%calculate_cg1) then
if (.not. allocated(CS%cg1)) call MOM_error(FATAL, &
"calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
if (CS%khth_use_ebt_struct) then
if (CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) then
if (.not. allocated(CS%ebt_struct)) call MOM_error(FATAL, &
"calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.")
if (CS%Resoln_use_ebt) then
Expand Down Expand Up @@ -1177,6 +1179,10 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"If true, uses the equivalent barotropic structure "//&
"as the vertical structure of thickness diffusivity.",&
default=.false.)
call get_param(param_file, mdl, "KD_GL90_USE_EBT_STRUCT", CS%kdgl90_use_ebt_struct, &
"If true, uses the equivalent barotropic structure "//&
"as the vertical structure of diffusivity in the GL90 scheme.",&
default=.false.)
call get_param(param_file, mdl, "KHTH_SLOPE_CFF", KhTh_Slope_Cff, &
"The nondimensional coefficient in the Visbeck formula "//&
"for the interface depth diffusivity", units="nondim", default=0.0)
Expand All @@ -1194,7 +1200,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
default=1.0e-17, units="s-1", scale=US%T_to_s)
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 .or. CS%khth_use_ebt_struct
CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct
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.)
Expand All @@ -1218,7 +1224,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.")
endif

if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct) then
if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) then
in_use = .true.
call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, &
"The depth below which N2 is monotonized to avoid stratification "//&
Expand Down Expand Up @@ -1568,7 +1574,7 @@ end subroutine VarMix_init
subroutine VarMix_end(CS)
type(VarMix_CS), intent(inout) :: CS

if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct) &
if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) &
deallocate(CS%ebt_struct)

if (CS%use_stored_slopes) then
Expand Down
Loading

0 comments on commit 7869030

Please sign in to comment.