Skip to content

Commit

Permalink
Updates GME by removing dependency on MEKE
Browse files Browse the repository at this point in the history 
Following changes have been made:

* add new argument (thickness_diffuse_CS) to horizontal_viscosity;
this is needed so that the GM coeff can be called from MOM_hor_visc.

* Deletes unnecessary calls to pass_vector

* Simplifies GME by removing dependecy on MEKE. GME is now set to be
some multiple of the GM coeff. A new runtime parameter (GME_efficiency)
can be used to control the strength of GME.
  • Loading branch information
@gustavo-marques
gustavo-marques committed Mar 13, 2020
1 parent c012019 commit 1ecade5
Show file tree
Hide file tree
Showing 2 changed files with 16 additions and 63 deletions.
5 changes: 3 additions & 2 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -681,7 +681,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
call cpu_clock_begin(id_clock_horvisc)
call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, &
MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, &
OBC=CS%OBC, BT=CS%barotropic_CSp)
OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp)
call cpu_clock_end(id_clock_horvisc)
if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")

Expand Down Expand Up @@ -1149,7 +1149,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
.not. query_initialized(CS%diffv,"diffv",restart_CS)) then
call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, &
G, GV, US, CS%hor_visc_CSp, &
OBC=CS%OBC, BT=CS%barotropic_CSp)
OBC=CS%OBC, BT=CS%barotropic_CSp, &
TD=thickness_diffuse_CSp)
else
if ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. &
(US%m_to_L * US%s_to_T_restart**2 /= US%m_to_L_restart * US%s_to_T**2) ) then
Expand Down
74 changes: 13 additions & 61 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ module MOM_hor_visc
!! v[is-2:ie+2,js-2:je+2]
!! h[is-1:ie+1,js-1:je+1]
subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
CS, OBC, BT)
CS, OBC, BT, TD)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
Expand All @@ -228,7 +228,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type
type(barotropic_CS), optional, pointer :: BT !< Pointer to a structure containing
!! barotropic velocities.

type(thickness_diffuse_CS), optional, pointer :: TD !< Pointer to a structure containing
!! thickness diffusivities.
! Local variables
real, dimension(SZIB_(G),SZJ_(G)) :: &
Del2u, & ! The u-compontent of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1]
Expand Down Expand Up @@ -418,8 +419,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
G%IdxCv(i,J-1) * vbtav(i,J-1))
enddo; enddo

call pass_vector(dudx_bt, dvdy_bt, G%Domain, stagger=BGRID_NE)

do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)
enddo ; enddo
Expand All @@ -432,8 +431,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
- ubtav(I,j)*G%IdxCu(I,j))
enddo ; enddo

call pass_vector(dvdx_bt, dudy_bt, G%Domain, stagger=AGRID)

if (CS%no_slip) then
do J=js-1,Jeq ; do I=is-1,Ieq
sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
Expand Down Expand Up @@ -1063,48 +1060,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo

if (CS%use_GME) then
if (CS%answers_2018) then
do j=js,je ; do i=is,ie
grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
(0.25*((dvdx(I,J) + dvdx(I-1,J-1)) + (dvdx(I,J-1) + dvdx(I-1,J))) )**2 + &
(0.25*((dudy(I,J) + dudy(I-1,J-1)) + (dudy(I,J-1) + dudy(I-1,J))) )**2)
max_diss_rate_h(i,j,k) = 2.0 * MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j))
enddo ; enddo
else ! This form is invariant to 90-degree rotations.
do j=js,je ; do i=is,ie
grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * ((dudx(i,j)**2 + dvdy(i,j)**2) + &
((0.25*((dvdx(I,J) + dvdx(I-1,J-1)) + (dvdx(I,J-1) + dvdx(I-1,J))) )**2 + &
(0.25*((dudy(I,J) + dudy(I-1,J-1)) + (dudy(I,J-1) + dudy(I-1,J))) )**2))
max_diss_rate_h(i,j,k) = 2.0 * MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j))
enddo ; enddo
endif

if (CS%answers_2018) then
do J = G%JscB, G%JecB ; do I = G%IscB, G%IecB
grad_vel_mag_q(I,J) = boundary_mask_q(I,J) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
(0.25*((dvdx(I,J)+dvdx(I-1,J-1)) + (dvdx(I,J-1)+dvdx(I-1,J))) )**2 + &
(0.25*((dudy(I,J)+dudy(I-1,J-1)) + (dudy(I,J-1)+dudy(I-1,J))) )**2)
call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G)

max_diss_rate_q(I,J,k) = 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)+ &
MEKE%MEKE(i,j+1)+MEKE%MEKE(i+1,j+1)) * sqrt(grad_vel_mag_q(I,J))
enddo ; enddo
else ! This form is rotationally invariant
do J = G%JscB, G%JecB ; do I = G%IscB, G%IecB
grad_vel_mag_q(I,J) = boundary_mask_q(I,J) * ((dudx(i,j)**2 + dvdy(i,j)**2) + &
((0.25*((dvdx(I,J)+dvdx(I-1,J-1)) + (dvdx(I,J-1)+dvdx(I-1,J))) )**2 + &
(0.25*((dudy(I,J)+dudy(I-1,J-1)) + (dudy(I,J-1)+dudy(I-1,J))) )**2))

max_diss_rate_q(I,J,k) = 0.5*((MEKE%MEKE(i,j) + MEKE%MEKE(i+1,j+1)) + &
(MEKE%MEKE(i+1,j) + MEKE%MEKE(i,j+1))) * sqrt(grad_vel_mag_q(I,J))
enddo ; enddo
endif

do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
if ((grad_vel_mag_bt_h(i,j)>0) .and. (max_diss_rate_h(i,j,k)>0)) then
GME_coeff = (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * CS%GME_efficiency*max_diss_rate_h(i,j,k) / &
grad_vel_mag_bt_h(i,j)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
if (grad_vel_mag_bt_h(i,j)>0) then
GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_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)))
else
GME_coeff = 1.0
GME_coeff = 0.0
endif

! apply mask
Expand All @@ -1117,10 +1080,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo

do J=js-1,Jeq ; do I=is-1,Ieq

if ((grad_vel_mag_bt_q(I,J)>0) .and. (max_diss_rate_q(I,J,k)>0)) then
GME_coeff = (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * CS%GME_efficiency*max_diss_rate_q(I,J,k) / &
grad_vel_mag_bt_q(I,J)
if (grad_vel_mag_bt_q(i,j)>0) then
GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_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)))
else
GME_coeff = 0.0
endif
Expand Down Expand Up @@ -1153,7 +1115,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

if (associated(MEKE%GME_snk)) then
do j=js,je ; do i=is,ie
FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_RZ * grad_vel_mag_bt_h(i,j)
FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 * grad_vel_mag_bt_h(i,j)
enddo ; enddo
endif

Expand Down Expand Up @@ -1395,8 +1357,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
! 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, use the MEKE module for calculating eddy kinetic energy.
! If false and USE_GME = True, issue a FATAL error.
logical :: default_2018_answers

character(len=64) :: inputdir, filename
Expand Down Expand Up @@ -1667,14 +1627,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// &
"cannot be used with SPLIT=False.")

call get_param(param_file, mdl, "USE_MEKE", use_MEKE, &
"If true, turns on the MEKE scheme which calculates\n"// &
"a sub-grid mesoscale eddy kinetic energy budget.", &
default=.false.)

if (.not. use_MEKE) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// &
"cannot be used with USE_MEKE=False.")

call get_param(param_file, mdl, "GME_H0", CS%GME_h0, &
"The strength of GME tapers quadratically to zero when the bathymetric "//&
"depth is shallower than GME_H0.", units="m", scale=US%m_to_Z, &
Expand Down

0 comments on commit 1ecade5

Please sign in to comment.