Skip to content

Commit

Permalink
Revert "Added biharmonic option for MEKE viscosity."
Browse files Browse the repository at this point in the history 
This reverts commit 6b32174.

Conflicts:
	src/parameterizations/lateral/MOM_MEKE.F90
	src/parameterizations/lateral/MOM_hor_visc.F90
  • Loading branch information
@sdbachman
sdbachman committed May 3, 2019
1 parent a6cd3fe commit 4b6ef52
Show file tree
Hide file tree
Showing 3 changed files with 29 additions and 80 deletions.
30 changes: 5 additions & 25 deletions src/parameterizations/lateral/MOM_MEKE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,6 @@ module MOM_MEKE
type(group_pass_type) :: pass_MEKE !< Type for group halo pass calls
type(group_pass_type) :: pass_Kh !< Type for group halo pass calls
type(group_pass_type) :: pass_Ku !< Type for group halo pass calls
type(group_pass_type) :: pass_Au !< Type for group halo pass calls
type(group_pass_type) :: pass_del2MEKE !< Type for group halo pass calls
end type MEKE_CS

Expand Down Expand Up @@ -566,11 +565,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
if (CS%viscosity_coeff/=0.) then
do j=js,je ; do i=is,ie
MEKE%Ku(i,j) = CS%viscosity_coeff*sqrt(2.*max(0.,MEKE%MEKE(i,j)))*LmixScale(i,j)
MEKE%Au(i,j) = CS%viscosity_coeff*sqrt(2.*max(0.,MEKE%MEKE(i,j)))*LmixScale(i,j)**3
enddo ; enddo
call cpu_clock_begin(CS%id_clock_pass)
call do_group_pass(CS%pass_Ku, G%Domain)
call do_group_pass(CS%pass_Au, G%Domain)
call cpu_clock_end(CS%id_clock_pass)
endif

Expand All @@ -581,7 +578,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
if (CS%id_Ut>0) call post_data(CS%id_Ut, sqrt(max(0.,2.0*MEKE%MEKE*barotrFac2)), CS%diag)
if (CS%id_Kh>0) call post_data(CS%id_Kh, MEKE%Kh, CS%diag)
if (CS%id_Ku>0) call post_data(CS%id_Ku, MEKE%Ku, CS%diag)
if (CS%id_Au>0) call post_data(CS%id_Au, MEKE%Au, CS%diag)
if (CS%id_KhMEKE_u>0) call post_data(CS%id_KhMEKE_u, Kh_u, CS%diag)
if (CS%id_KhMEKE_v>0) call post_data(CS%id_KhMEKE_v, Kh_v, CS%diag)
if (CS%id_src>0) call post_data(CS%id_src, src, CS%diag)
Expand Down Expand Up @@ -723,8 +719,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, SN_u, SN_v, drag_rate_visc, I_mass)
else
EKE = 0.
endif
! MEKE%MEKE(i,j) = EKE
MEKE%MEKE(i,j) = (G%Zd_to_m*G%bathyT(i,j)*SN / (8*CS%cdrag))**2
MEKE%MEKE(i,j) = EKE
enddo ; enddo

end subroutine MEKE_equilibrium
Expand Down Expand Up @@ -844,7 +839,7 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
type(MOM_restart_CS), pointer :: restart_CS !< Restart control structure for MOM_MEKE.
! Local variables
integer :: is, ie, js, je, isd, ied, jsd, jed, nz
logical :: laplacian, biharmonic, useVarMix, coldStart
logical :: laplacian, useVarMix, coldStart
! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mdl = "MOM_MEKE" ! This module's name.
Expand Down Expand Up @@ -1009,10 +1004,8 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
"the velocity field to the bottom stress.", units="nondim", &
default=0.003)
call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)
call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.)

if (CS%viscosity_coeff/=0. .and. .not. laplacian .and. .not. biharmonic) call MOM_error(FATAL, &
"Either LAPLACIAN or BIHARMONIC must be true if MEKE_VISCOSITY_COEFF is true.")
if (CS%viscosity_coeff/=0. .and. .not. laplacian) call MOM_error(FATAL, &
"LAPLACIAN must be true if MEKE_VISCOSITY_COEFF is true.")

call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)

Expand All @@ -1034,10 +1027,6 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
call create_group_pass(CS%pass_Ku, MEKE%Ku, G%Domain)
call do_group_pass(CS%pass_Ku, G%Domain)
endif
if (associated(MEKE%Au)) then
call create_group_pass(CS%pass_Au, MEKE%Au, G%Domain)
call do_group_pass(CS%pass_Au, G%Domain)
endif
if (allocated(CS%del2MEKE)) then
call create_group_pass(CS%pass_del2MEKE, CS%del2MEKE, G%Domain)
call do_group_pass(CS%pass_del2MEKE, G%Domain)
Expand All @@ -1054,9 +1043,6 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
CS%id_Ku = register_diag_field('ocean_model', 'MEKE_KU', diag%axesT1, Time, &
'MEKE derived lateral viscosity', 'm2 s-1')
if (.not. associated(MEKE%Ku)) CS%id_Ku = -1
CS%id_Au = register_diag_field('ocean_model', 'MEKE_AU', diag%axesT1, Time, &
'MEKE derived lateral biharmonic viscosity', 'm4 s-1')
if (.not. associated(MEKE%Au)) CS%id_Au = -1
CS%id_Ue = register_diag_field('ocean_model', 'MEKE_Ue', diag%axesT1, Time, &
'MEKE derived eddy-velocity scale', 'm s-1')
if (.not. associated(MEKE%MEKE)) CS%id_Ue = -1
Expand Down Expand Up @@ -1163,14 +1149,9 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
allocate(MEKE%Rd_dx_h(isd:ied,jsd:jed)) ; MEKE%Rd_dx_h(:,:) = 0.0
if (MEKE_viscCoeff/=0.) then
allocate(MEKE%Ku(isd:ied,jsd:jed)) ; MEKE%Ku(:,:) = 0.0
vd = var_desc("MEKE_Ku", "m2 s-1", hor_grid='h', z_grid='1', &
vd = var_desc("MEKE_Ah", "m2 s-1", hor_grid='h', z_grid='1', &
longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy")
call register_restart_field(MEKE%Ku, vd, .false., restart_CS)

allocate(MEKE%Au(isd:ied,jsd:jed)) ; MEKE%Au(:,:) = 0.0
vd = var_desc("MEKE_Au", "m4 s-1", hor_grid='h', z_grid='1', &
longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy")
call register_restart_field(MEKE%Au, vd, .false., restart_CS)
endif

end subroutine MEKE_alloc_register_restart
Expand All @@ -1191,7 +1172,6 @@ subroutine MEKE_end(MEKE, CS)
if (associated(MEKE%GME_snk)) deallocate(MEKE%GME_snk)
if (associated(MEKE%Kh)) deallocate(MEKE%Kh)
if (associated(MEKE%Ku)) deallocate(MEKE%Ku)
if (associated(MEKE%Au)) deallocate(MEKE%Au)
if (allocated(CS%del2MEKE)) deallocate(CS%del2MEKE)
deallocate(MEKE)

Expand Down
1 change: 0 additions & 1 deletion src/parameterizations/lateral/MOM_MEKE_types.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@ module MOM_MEKE_types
real, dimension(:,:), pointer :: Ku => NULL() !< The MEKE-derived lateral viscosity coefficient in m2 s-1.
!! This viscosity can be negative when representing backscatter
!! from unresolved eddies (see Jansen and Held, 2014).
real, dimension(:,:), pointer :: Au => NULL() !< The MEKE-derived lateral biharmonic viscosity coefficient in m4 s-1.
! Parameters
real :: KhTh_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTh, nondim
real :: KhTr_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTr, nondim.
Expand Down
78 changes: 24 additions & 54 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,10 +245,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points (m-1 s-1)
grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points (m-1 s-1)
grad_div_mag_h, & ! Magnitude of divergence gradient at h-points (m-1 s-1)
dudx, dvdy, & ! components in the horizontal tension (s-1)
grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points (s-2)
grad_vel_mag_bt_h, & ! Magnitude of the barotropic velocity gradient tensor squared at h-points (s-2)
max_diss_rate_bt ! maximum possible energy dissipated by barotropic lateral friction (m2 s-3)
dudx, dvdy ! components in the horizontal tension (s-1)

real, dimension(SZIB_(G),SZJB_(G)) :: &
dvdx, dudy, & ! components in the shearing strain (s-1)
Expand All @@ -265,9 +262,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points (m-1 s-1)
grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points (m-1 s-1)
grad_div_mag_q, & ! Magnitude of divergence gradient at q-points (m-1 s-1)
grad_vel_mag_q, & ! Magnitude of the velocity gradient tensor squared at q-points (s-2)
hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses, in H; This form guarantees that hq/hu < 4.
grad_vel_mag_bt_q ! Magnitude of the barotropic velocity gradient tensor squared at q-points (s-2)
hq ! harmonic mean of the harmonic means of the u- & v point thicknesses, in H; This form guarantees that hq/hu < 4.

real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
Ah_q, & ! biharmonic viscosity at corner points (m4/s)
Expand All @@ -282,14 +277,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
Ah_h, & ! biharmonic viscosity at thickness points (m4/s)
Kh_h, & ! Laplacian viscosity at thickness points (m2/s)
diss_rate, & ! MKE dissipated by parameterized shear production (m2 s-3)
max_diss_rate, & ! maximum possible energy dissipated by lateral friction (m2 s-3)
target_diss_rate_GME, & ! target amount of energy to add via GME (m2 s-3)
FrictWork, & ! work done by MKE dissipation mechanisms (W/m2)
FrictWork_diss, & ! negative definite work done by MKE dissipation mechanisms (W/m2)
FrictWorkMax, & ! maximum possible work done by MKE dissipation mechanisms (W/m2)
FrictWork_GME, & ! work done by GME (W/m2)
target_FrictWork_GME, & ! target amount of work for GME to do (W/m2)
FrictWork, & ! energy flux by parameterized shear production (W/m2)
FrictWork_diss, & ! MKE dissipated by parameterized shear production (m3 s-3)
FrictWorkMax, & ! maximum possible energy dissipated by lateral friction (m3 s-3)
FrictWork_GME, & ! MKE added by parameterized shear production in GME (m3 s-3)
target_FrictWork_GME, & ! target amount of energy to add via GME (m3 s-3)
div_xx_h ! horizontal divergence (s-1)
!real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
Expand Down Expand Up @@ -334,7 +326,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
logical :: find_FrictWork
logical :: apply_OBC = .false.
logical :: use_MEKE_Ku
logical :: use_MEKE_Au
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
integer :: i, j, k, n
real :: inv_PI3, inv_PI6
Expand Down Expand Up @@ -375,13 +366,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

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

!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h, &
!$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, &
!$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
!$OMP find_FrictWork,FrictWork,use_MEKE_Ku,
!$OMP use_MEKE_Au, MEKE, hq, &
!$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE, hq, &
!$OMP mod_Leith, legacy_bound, div_xx_h, vort_xy_q) &
!$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
!$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
Expand Down Expand Up @@ -443,8 +432,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
endif

! Get thickness diffusivity for use in GME
! call thickness_diffuse_get_KH(thickness_diffuse, KH_u_GME, KH_v_GME, G)

! call thickness_diffuse_get_KH(thickness_diffuse, KH_u_GME, KH_v_GME, G)
do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
grad_vel_mag_bt_h(i,j) = 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 @@ -464,11 +452,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
enddo ; enddo


! halo updates (presently not used since GME is now hooked to MEKE)
! halo updates (presently not used since GME is now hooked to MEKE)
! call pass_vector(KH_u_GME, KH_v_GME, G%Domain)
! call pass_vector(VarMix%slope_x, VarMix%slope_y, G%Domain)
! call pass_vector(VarMix%N2_u, VarMix%N2_v, G%Domain)

endif ! use_GME

do k=1,nz
Expand All @@ -493,14 +480,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
enddo ; enddo


if ((find_FrictWork) .or. (CS%use_GME)) then
do j=js,je ; do i=is,ie
grad_vel_mag_h(i,j) = (dudx(i,j)**2 + dvdy(i,j)**2 + &
(0.25*(dvdx(I,J)+dvdx(I-1,J)+dvdx(I,J-1)+dvdx(I-1,J-1)) )**2 + &
(0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2)
enddo ; enddo
endif

! Interpolate the thicknesses to velocity points.
! The extra wide halos are to accommodate the cross-corner-point projections
! in OBCs, which are not ordinarily be necessary, and might not be necessary
Expand Down Expand Up @@ -855,7 +834,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

str_xx(i,j) = -Kh * sh_xx(i,j)
else ! not Laplacian
Kh_h(i,j,k) = 0.0
str_xx(i,j) = 0.0
endif ! Laplacian

Expand Down Expand Up @@ -887,8 +865,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
Ah = CS%Ah_bg_xx(i,j)
endif ! Smagorinsky_Ah or Leith_Ah

if (use_MEKE_Au) Ah = Ah + MEKE%Au(i,j) ! *Add* the MEKE contribution

if (CS%better_bound_Ah) then
Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xx(i,j))
endif
Expand All @@ -905,8 +881,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
CS%DX_dyT(i,j) *(G%IdxCv(i,J)*v0(i,J) - G%IdxCv(i,J-1)*v0(i,J-1)))
bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j))

else
Ah_h(i,j,k) = 0.0
endif ! biharmonic

enddo ; enddo
Expand Down Expand Up @@ -1053,12 +1027,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
else
Ah = CS%Ah_bg_xy(I,J)
endif ! Smagorinsky_Ah or Leith_Ah

if (use_MEKE_Au) then ! *Add* the MEKE contribution
Ah = Ah + 0.25*( (MEKE%Au(I,J)+MEKE%Au(I+1,J+1)) &
+(MEKE%Au(I+1,J)+MEKE%Au(I,J+1)) )
endif

if (CS%better_bound_Ah) then
Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xy(I,J))
endif
Expand All @@ -1085,26 +1053,26 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

do j=js,je ; do i=is,ie
! Diagnose -Kh * |del u|^2 - Ah * |del^2 u|^2
diss_rate(i,j,k) = -Kh_h(i,j,k) * grad_vel_mag_h(i,j) - &
FrictWork_diss(i,j,k) = -Kh_h(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 * &
(dudx(i,j)**2 + dvdy(i,j)**2 + &
(0.25*(dvdx(I,J)+dvdx(I-1,J)+dvdx(I,J-1)+dvdx(I-1,J-1)) )**2 + &
(0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2) - &
Ah_h(i,j,k) * ((0.5*(u0(I,j) + u0(I-1,j)))**2 + &
(0.5*(v0(i,J) + v0(i,J-1)))**2)
FrictWork_diss(i,j,k) = diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2

if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then
! This is the maximum possible amount of energy that can be converted
! per unit time, according to theory (multiplied by h)
max_diss_rate(i,j,k) = 2.0*MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j)) * (MIN(G%bathyT(i,j)/H0,1.0)**2)

FrictWorkMax(i,j,k) = max_diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2

! Determine how much work GME needs to do to reach the "target" ratio between
! the amount of work actually done and the maximum allowed by theory. Note that
! we need to add the FrictWork done by the dissipation operators, since this work
! is done only for numerical stability and is therefore spurious
! Determine how much work GME needs to do to reach the "target" ratio between
! the amount of work actually done and the maximum allowed by theory. Note that
! we need to add the FrictWork done by the dissipation operators, since this work
! is done only for numerical stability and is therefore spurious
if (CS%use_GME) then
target_diss_rate_GME(i,j,k) = FWfrac * max_diss_rate(i,j,k) - diss_rate(i,j,k)
target_FrictWork_GME(i,j,k) = target_diss_rate_GME(i,j,k) * h(i,j,k) * GV%H_to_kg_m2
endif
target_FrictWork_GME(i,j,k) = FWfrac * FrictWorkMax(i,j,k) - FrictWork_diss(i,j,k)
endif

endif ; endif

Expand All @@ -1131,6 +1099,9 @@ 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 = 2e5 ! 1e6

! simple way to limit this coeff
GME_coeff = MIN(GME_coeff,GME_coeff_limiter)

if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff
Expand Down Expand Up @@ -1302,9 +1273,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
else
do j=js,je ; do i=is,ie
! MEKE%mom_src now is sign definite because it only uses the dissipation
MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + MAX(-FrictWorkMax(i,j,k),FrictWork_diss(i,j,k))
MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork_diss(i,j,k)
enddo ; enddo

if (CS%use_GME) then
if (associated(MEKE%GME_snk)) then
do j=js,je ; do i=is,ie
Expand Down

0 comments on commit 4b6ef52

Please sign in to comment.