diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 77d431442b..62337301a1 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -87,6 +87,9 @@ module MOM_hor_visc !! that gets backscattered in the Leith+E scheme. [nondim] logical :: smooth_Ah !< If true (default), then Ah and m_leithy are smoothed. !! This smoothing requires a lot of blocking communication. + logical :: taper_leithy !< If true, backscatter coeff is tapered to zero with depth + real :: leithy_depth !< If tapering leith+E, taper is applied below this depth [Z ~> m] + real :: leithy_width !< If tapering leith+E, backscatter is zero below leithy_depth+leithy_width [Z ~> m] logical :: use_QG_Leith_visc !< If true, use QG Leith nonlinear eddy viscosity. !! KH is the background value. logical :: bound_Coriolis !< If true & SMAGORINSKY_AH is used, the biharmonic @@ -471,6 +474,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, hrat_min, & ! h_min divided by the thickness at the stress point (h or q) [nondim] visc_bound_rem ! fraction of overall viscous bounds that remain to be applied (h or q) [nondim] + ! Layer depth. Used to taper Leith+E backscatter coefficient with depth + real, allocatable :: zc(:,:,:) ! depth at center of h cell [Z ~> m] + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -651,6 +657,16 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, call smooth_x9_uv(CS, G, u_smooth(:,:,k), v_smooth(:,:,k), zero_land=.false.) enddo call pass_vector(u_smooth, v_smooth, G%Domain) + ! If tapering the backscatter, calculate depths now + if (CS%taper_leithy) then + ! allocate zc + allocate(zc(SZI_(G),SZJ_(G),SZK_(GV))) ; zc(:,:,:) = 0.0 + ! Compute zc. Not actual cell centers because it starts at 0 rather than at SSH. + zc(:,:,1) = 0.5 * h(:,:,1) + do k=2,nz + zc(:,:,k) = zc(:,:,k-1) + 0.5 * (h(:,:,k-1) + h(:,:,k)) + enddo + endif endif if (CS%use_QG_Leith_visc .and. ((CS%Leith_Kh) .or. (CS%Leith_Ah))) then @@ -674,7 +690,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_bh, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & - !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, hu_cont, hv_cont, STOCH & + !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, hu_cont, hv_cont, STOCH, zc & !$OMP ) & !$OMP private( & !$OMP i, j, k, n, tmp, & @@ -1320,6 +1336,10 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, m_leithy(i,j) = CS%m_leithy_max(i,j) endif endif + if (CS%taper_leithy) then + ! Multiply m_leithy by taper function of depth + m_leithy(:,:) = m_leithy(:,:) * leithy_taper_function(CS, zc(:,:,k)) + endif enddo ; enddo if (CS%smooth_Ah) then @@ -2286,6 +2306,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, CS%dx2h, CS%dy2h, CS%dx2q, CS%dy2q) endif + if (allocated(zc)) deallocate(zc) + end subroutine horizontal_viscosity !> Allocates space for and calculates static variables used by horizontal_viscosity. @@ -2655,6 +2677,19 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "If true, Ah and m_leithy are smoothed within Leith+E. This requires "//& "lots of blocking communications, which can be expensive", & default=.true., do_not_log=.not.CS%use_Leithy) + call get_param(param_file, mdl, "TAPER_LEITHY", CS%taper_leithy, & + "If true, Leith+E c_K coefficient is tapered to zero"//& + "below a threshold depth", & + default=.false., do_not_log=.not.CS%use_Leithy) + if (CS%taper_leithy) then + call get_param(param_file, mdl, "LEITHY_DEPTH", CS%leithy_depth, & + "Leith+E backscatter starts tapering below this depth.", & + units="m", scale=US%m_to_Z, default=800.0) + call get_param(param_file, mdl, "LEITHY_WIDTH", CS%leithy_width, & + "Leith+E backscatter is zero below LEITHY_DEPTH+LEITHY_WIDTH.", & + units="m", scale=US%m_to_Z, default=400.0) + if (CS%leithy_width <= 0.0) call MOM_error(FATAL,"ERROR: LEITHY_WIDTH must be positive ") + endif if (CS%use_GME .and. .not.split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") @@ -3265,6 +3300,26 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) end subroutine hor_visc_init +!> leithy_taper_function returns 1 if zc is shallower than leithy_depth; 0 if deeper than +!! leithy_depth+leithy_width; and an interpolating cubic spline in between. +function leithy_taper_function(CS, zc) + type(hor_visc_CS), intent(in) :: CS !< Control structure for horizontal viscosity + real, intent(in) :: zc(:,:) !< depth of h-cell centersi [Z ~> m] + real :: leithy_taper_function(size(zc,dim=1),size(zc,dim=2)) ! Taper function evaluated at zc [nondim] + + ! Local variables + real :: x(size(zc,dim=1),size(zc,dim=2)) ! 0 at top of transition and 1 at bottom [nondim] + + x = (zc - CS%leithy_depth) / CS%leithy_width + where (zc <= CS%leithy_depth) + leithy_taper_function = 1.0 + elsewhere (zc >= CS%leithy_depth + CS%leithy_width) + leithy_taper_function = 0.0 + elsewhere + leithy_taper_function = (x - 1.0)**2 * (1.0 + 2 * x) + end where +end function leithy_taper_function + !> hor_visc_vel_stencil returns the horizontal viscosity input velocity stencil size function hor_visc_vel_stencil(CS) result(stencil) type(hor_visc_CS), intent(in) :: CS !< Control structure for horizontal viscosity