Skip to content

Commit

Permalink
Added a new flag to MOM_lateral_mixing_coeffs called "USE_VISBECK", w…
Browse files Browse the repository at this point in the history 
…hich

calculates the GM coefficient using the Visbeck et al. (1997) formulation. Previously
this was being done whenever VarMix was enabled, but I wanted to separate VarMix from
this diffusivity call.
  • Loading branch information
@sdbachman
sdbachman committed Apr 16, 2019
1 parent 6542122 commit 8df6efd
Show file tree
Hide file tree
Showing 4 changed files with 61 additions and 15 deletions.
8 changes: 4 additions & 4 deletions src/parameterizations/lateral/MOM_MEKE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -491,10 +491,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
endif ! MEKE_KH>=0
!$OMP end parallel

! Ensure that MEKE is non-negative
do j=js,je ; do i=is,ie
MEKE%MEKE(i,j) = MAX(0.0, MEKE%MEKE(i,j))
enddo ; enddo
! ! Ensure that MEKE is non-negative
! do j=js,je ; do i=is,ie
! MEKE%MEKE(i,j) = MAX(0.0, MEKE%MEKE(i,j))
! enddo ; enddo


call cpu_clock_begin(CS%id_clock_pass)
Expand Down
27 changes: 22 additions & 5 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1065,14 +1065,22 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * 0.25*(KH_u_GME(I,j,k) + KH_u_GME(I-1,j,k) + &
KH_v_GME(i,J,k) + KH_v_GME(i,J-1,k)) * &
sqrt(0.5*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)),0.0) * &
(0.5*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)),0.0) * &
( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I-1,j,k)) )**2 + &
(0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i,J-1,k)) )**2 ) / &
( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
(0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
(0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
epsilon))

! GME_coeff = 2.0 * (MIN(G%bathyT(i,j)/H0,1.0)**2) * 0.25*(KH_u_GME(I,j,k) + KH_u_GME(I-1,j,k) + &
! KH_v_GME(i,J,k) + KH_v_GME(i,J-1,k)) * &
! sqrt(0.5*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)),0.0)) / &
! sqrt( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
! (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
! (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
! epsilon)

! GME_coeff = 2.0 * MAX(0.0,MEKE%MEKE(i,j)) / &
! SQRT( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
! (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
Expand All @@ -1089,7 +1097,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
! ! apply mask
! GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))

GME_coeff_limiter = 0.0
GME_coeff_limiter = 1e6 ! 1e8
! GME_coeff_limiter = 2.0 * MAX(0.0,MEKE%MEKE(i,j)) / sqrt(dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
! (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
! (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + epsilon)
Expand Down Expand Up @@ -1122,7 +1130,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * 0.25*(KH_u_GME(I,j,k) + KH_u_GME(I,j+1,k) + &
KH_v_GME(i,J,k) + KH_v_GME(i+1,J,k)) * &
sqrt( 0.25*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I,j+1,k) + &
( 0.25*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I,j+1,k) + &
VarMix%N2_v(i,J,k)+VarMix%N2_v(i+1,J,k)),0.0) * &
( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I,j+1,k)) )**2 + &
(0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i+1,J,k)) )**2 ) / &
Expand All @@ -1131,6 +1139,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
(0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
epsilon))

! GME_coeff = 2.0 * (MIN(G%bathyT(i,j)/H0,1.0)**2) * 0.25*(KH_u_GME(I,j,k) + KH_u_GME(I,j+1,k) + &
! KH_v_GME(i,J,k) + KH_v_GME(i+1,J,k)) * &
! sqrt( 0.25*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I,j+1,k) + &
! VarMix%N2_v(i,J,k)+VarMix%N2_v(i+1,J,k)),0.0)) / &
! sqrt(dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
! (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
! (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
! epsilon)

! GME_coeff = 2.0* MAX(0.0,MEKE%MEKE(i,j)) / &
! SQRT( dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
! (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
Expand All @@ -1147,7 +1164,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
! apply mask
GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))

GME_coeff_limiter = 0.0
GME_coeff_limiter = 1e6 !1e8
! GME_coeff_limiter = 2.0 * MAX(0.0,MEKE%MEKE(i,j)) / sqrt( dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
! (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
! (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + epsilon)
Expand Down Expand Up @@ -1288,7 +1305,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
enddo ; enddo
else
do j=js,je ; do i=is,ie
MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k)
MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork_diss(i,j,k)
enddo ; enddo
endif
endif ; endif
Expand Down
4 changes: 4 additions & 0 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ module MOM_lateral_mixing_coeffs
KH_v_QG !< QG Leith GM coefficient at v-points (m2 s-1)

! Parameters
logical :: use_Visbeck !< Use Visbeck formulation for thickness diffusivity
integer :: VarMix_Ktop !< Top layer to start downward integrals
real :: Visbeck_L_scale !< Fixed length scale in Visbeck formula
real :: Res_coef_khth !< A non-dimensional number that determines the function
Expand Down Expand Up @@ -928,6 +929,9 @@ subroutine VarMix_init(Time, G, GV, param_file, diag, CS)
"not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, \n"//&
"this is set to true regardless of what is in the \n"//&
"parameter file.", default=.false.)
call get_param(param_file, mdl, "USE_VISBECK", CS%use_Visbeck,&
"If true, use the Visbeck et al. (1997) formulation for \n"//&
"thickness diffusivity.", default=.false.)
call get_param(param_file, mdl, "RESOLN_SCALED_KH", CS%Resoln_scaled_Kh, &
"If true, the Laplacian lateral viscosity is scaled away \n"//&
"when the first baroclinic deformation radius is well \n"//&
Expand Down
37 changes: 31 additions & 6 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected, in H.
real, dimension(:,:), pointer :: cg1 => null() !< Wave speed (m/s)
logical :: use_VarMix, Resoln_scaled, use_stored_slopes, khth_use_ebt_struct
logical :: use_VarMix, Resoln_scaled, use_stored_slopes, khth_use_ebt_struct, use_Visbeck
logical :: use_QG_Leith
integer :: i, j, k, is, ie, js, je, nz
real :: hu(SZI_(G), SZJ_(G)) ! u-thickness (H)
Expand Down Expand Up @@ -151,6 +151,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS
Resoln_scaled = VarMix%Resoln_scaled_KhTh
use_stored_slopes = VarMix%use_stored_slopes
khth_use_ebt_struct = VarMix%khth_use_ebt_struct
use_Visbeck = VarMix%use_Visbeck
use_QG_Leith = VarMix%use_QG_Leith_GM
if (associated(VarMix%cg1)) cg1 => VarMix%cg1
else
Expand Down Expand Up @@ -183,7 +184,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS

if (use_VarMix) then
!$OMP do
if (.not. use_QG_Leith) then
if (use_Visbeck) then
do j=js,je ; do I=is-1,ie
Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + CS%KHTH_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j)
enddo ; enddo
Expand Down Expand Up @@ -255,7 +256,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS

if (use_VarMix) then
!$OMP do
if (.not. use_QG_Leith) then
if (use_Visbeck) then
do J=js-1,je ; do i=is,ie
Khth_Loc(i,j) = Khth_Loc(i,j) + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J)
enddo ; enddo
Expand Down Expand Up @@ -328,6 +329,16 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS
do K=1,nz+1 ; do J=js-1,je ; do i=is,ie ; int_slope_v(i,J,K) = 0.0 ; enddo ; enddo ; enddo
!$OMP end parallel

! if (associated(MEKE)) then
! do j=js,je ; do i=is,ie
! MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + 0.25*(Kh_u(I,j,k) + KH_u(I-1,j,k) + &
! KH_v(i,J,k) + KH_v(i,J-1,k))
!! 0.25*MAX(VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)+VarMix%N2_v(i,J,k)+VarMix%N2_v(i,J-1,k),0.0) * &
!! ( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I-1,j,k)) )**2 + &
!! (0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i,J-1,k)) )**2 )
! enddo; enddo
! endif

if (CS%detangle_interfaces) then
call add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, &
CS, int_slope_u, int_slope_v)
Expand All @@ -351,7 +362,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS
! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S
if (use_stored_slopes) then
call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, CS, &
int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y)
int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y, &
VarMix%N2_u, VarMix%N2_v)
else
call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, CS, &
int_slope_u, int_slope_v)
Expand Down Expand Up @@ -445,7 +457,7 @@ end subroutine thickness_diffuse
!! Fluxes are limited to give positive definite thicknesses.
!! Called by thickness_diffuse().
subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, &
CS, int_slope_u, int_slope_v, slope_x, slope_y)
CS, int_slope_u, int_slope_v, slope_x, slope_y, N2_x, N2_y)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness in H (m or kg/m2)
Expand All @@ -471,7 +483,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
!! density gradients.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: slope_x !< Isopycnal slope at u-points
real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: slope_y !< Isopycnal slope at v-points

real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: N2_x !< Brunt-Vaisala frequency at
!! u points (s-2)
real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: N2_y !< Brunt-Vaisala frequency at
!! v points (s-2)
! Local variables
real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
T, & ! The temperature (or density) in C, with the values in
Expand Down Expand Up @@ -556,6 +571,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics
logical :: present_int_slope_u, present_int_slope_v
logical :: present_slope_x, present_slope_y, calc_derivatives
logical :: present_N2_x, present_N2_y
integer :: is, ie, js, je, nz, IsdB
integer :: i, j, k
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke ; IsdB = G%IsdB
Expand All @@ -573,6 +589,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
present_int_slope_v = PRESENT(int_slope_v)
present_slope_x = PRESENT(slope_x)
present_slope_y = PRESENT(slope_y)
present_N2_x = PRESENT(N2_x)
present_N2_y = PRESENT(N2_y)

nk_linear = max(GV%nkml, 1)

Expand Down Expand Up @@ -1181,6 +1199,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
if (associated(CS%GMwork)) CS%GMwork(i,j) = Work_h
if (associated(MEKE)) then ; if (associated(MEKE%GM_src)) then
MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h
! if (present_N2_x .and. present_N2_y) &
! MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + 0.25*(Kh_u(I,j,k) + KH_u(I-1,j,k) + &
! KH_v(i,J,k) + KH_v(i,J-1,k)) * &
! 0.25*MAX(N2_x(I,j,k)+N2_x(I-1,j,k)+N2_y(i,J,k)+N2_y(i,J-1,k),0.0) * &
! ( (0.5*(slope_x(I,j,k)+slope_x(I-1,j,k)) )**2 + &
! (0.5*(slope_y(i,J,k)+slope_y(i,J-1,k)) )**2 )

endif ; endif
enddo ; enddo ; endif

Expand Down

0 comments on commit 8df6efd

Please sign in to comment.