Skip to content

Commit

Permalink
Added Jansen et al. (2015) version of MEKE dissipation. Added an alte…
Browse files Browse the repository at this point in the history 
…rnative

way of calculating PE-to_MEKE energy conversion.
  • Loading branch information
@sdbachman
sdbachman committed Apr 30, 2019
1 parent 172ae02 commit c24bfb3
Show file tree
Hide file tree
Showing 2 changed files with 100 additions and 41 deletions.
92 changes: 60 additions & 32 deletions src/parameterizations/lateral/MOM_MEKE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,9 @@ module MOM_MEKE
real :: MEKE_min_gamma!< Minimum value of gamma_b^2 allowed (non-dim)
real :: MEKE_Ct !< Coefficient in the \f$\gamma_{bt}\f$ expression (non-dim)
logical :: visc_drag !< If true use the vertvisc_type to calculate bottom drag.
logical :: Jansen15_drag !< If true use the bottom drag formulation from Jansen et al. (2015)
logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather
!! than the streamfunction for the MEKE GM source term.
logical :: Rd_as_max_scale !< If true the length scale can not exceed the
!! first baroclinic deformation radius.
logical :: use_old_lscale !< Use the old formula for mixing length scale.
Expand Down Expand Up @@ -114,6 +117,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
LmixScale, & ! Square of eddy mixing length, in m2.
barotrFac2, & ! Ratio of EKE_barotropic / EKE (nondim)/
bottomFac2 ! Ratio of EKE_bottom / EKE (nondim)/

real, dimension(SZIB_(G),SZJ_(G)) :: &
MEKE_uflux, & ! The zonal diffusive flux of MEKE, in kg m2 s-3.
Kh_u, & ! The zonal diffusivity that is actually used, in m2 s-1.
Expand Down Expand Up @@ -289,9 +293,14 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)

if (associated(MEKE%GM_src)) then
!$OMP do
do j=js,je ; do i=is,ie
src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j)
enddo ; enddo
if (CS%GM_src_alt) then
do j=js,je ; do i=is,ie
src(i,j) = src(i,j) - CS%MEKE_GMcoeff*MEKE%GM_src(i,j) / MAX(1.0,G%bathyT(i,j))
enddo ; enddo
else
do j=js,je ; do i=is,ie
src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j)
enddo ; enddo
endif

! Increase EKE by a full time-steps worth of source
Expand All @@ -303,22 +312,37 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
if (use_drag_rate) then
! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
!$OMP do
do j=js,je ; do i=is,ie
drag_rate(i,j) = (Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
if (CS%Jansen15_drag) then
do j=js,je ; do i=is,ie
drag_rate(i,j) = (cdrag2/MAX(1.0,G%bathyT(i,j))) * sqrt(CS%MEKE_Uscale**2 + drag_rate_visc(i,j)**2 + &
2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) * 2.0 * bottomFac2(i,j)*MEKE%MEKE(i,j)
enddo ; enddo
else
do j=js,je ; do i=is,ie
drag_rate(i,j) = (Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
+ cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
enddo ; enddo
enddo ; enddo
endif
endif

! First stage of Strang splitting
!$OMP do
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j)<0.) ldamping = 0.
! notice that the above line ensures a damping only if MEKE is positive,
! while leaving MEKE unchanged if it is negative
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
enddo ; enddo
if (CS%Jansen15_drag) then
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j)
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) - MIN(MEKE%MEKE(i,j),sdt_damp*drag_rate(i,j))
MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
enddo ; enddo
else
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j)<0.) ldamping = 0.
! notice that the above line ensures a damping only if MEKE is positive,
! while leaving MEKE unchanged if it is negative
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
enddo ; enddo
endif
!$OMP end parallel

if (CS%MEKE_KH >= 0.0 .or. CS%KhMEKE_FAC > 0.0 .or. CS%MEKE_K4 >= 0.0) then
Expand Down Expand Up @@ -473,30 +497,28 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
! Recalculate the drag rate, since MEKE has changed.
if (use_drag_rate) then
!$OMP do
do j=js,je ; do i=is,ie
drag_rate(i,j) = (Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
+ cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
enddo ; enddo
if (CS%Jansen15_drag) then
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j)
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) - MIN(MEKE%MEKE(i,j),sdt_damp*drag_rate(i,j))
MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
enddo ; enddo
else
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j)<0.) ldamping = 0.
! notice that the above line ensures a damping only if MEKE is positive,
! while leaving MEKE unchanged if it is negative
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
enddo ; enddo
endif
endif
!$OMP do
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j)<0.) ldamping = 0.
! notice that the above line ensures a damping only if MEKE is positive,
! while leaving MEKE unchanged if it is negative
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
MEKE_decay(i,j) = 0.5 * G%mask2dT(i,j) * (MEKE_decay(i,j) + ldamping)
enddo ; enddo
endif
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


call cpu_clock_begin(CS%id_clock_pass)
call do_group_pass(CS%pass_MEKE, G%Domain)
call cpu_clock_end(CS%id_clock_pass)
Expand Down Expand Up @@ -888,6 +910,12 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
call get_param(param_file, mdl, "MEKE_USCALE", CS%MEKE_Uscale, &
"The background velocity that is combined with MEKE to \n"//&
"calculate the bottom drag.", units="m s-1", default=0.0)
call get_param(param_file, mdl, "MEKE_JANSEN15_DRAG", CS%Jansen15_drag, &
"If true, use the bottom drag formulation from Jansen et al. (2015) \n"//&
"to calculate the drag acting on MEKE.", default=.false.)
call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", CS%GM_src_alt, &
"If true, use the GM energy conversion form S^2*N^2*kappa rather \n"//&
"than the streamfunction for the MEKE GM source term.", default=.false.)
call get_param(param_file, mdl, "MEKE_VISC_DRAG", CS%visc_drag, &
"If true, use the vertvisc_type to calculate the bottom \n"//&
"drag acting on MEKE.", default=.true.)
Expand Down
49 changes: 40 additions & 9 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ module MOM_thickness_diffuse
use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
use MOM_diag_mediator, only : diag_update_remap_grids
use MOM_domains, only : pass_var, CORNER, pass_vector
use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
use MOM_EOS, only : calculate_density, calculate_density_derivs
use MOM_file_parser, only : get_param, log_version, param_file_type
Expand Down Expand Up @@ -244,7 +245,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS

!$OMP do
if (CS%use_GME_thickness_diffuse) then
do k=1,nz ; do j=js,je ; do I=is-1,ie
do k=1,nz+1 ; do j=js,je ; do I=is-1,ie
CS%KH_u_GME(I,j,k) = KH_u(I,j,k)
enddo ; enddo ; enddo
endif
Expand Down Expand Up @@ -318,7 +319,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS

!$OMP do
if (CS%use_GME_thickness_diffuse) then
do k=1,nz ; do j=js-1,je ; do I=is,ie
do k=1,nz+1 ; do j=js-1,je ; do I=is,ie
CS%KH_v_GME(I,j,k) = KH_v(I,j,k)
enddo ; enddo ; enddo
endif
Expand Down Expand Up @@ -484,6 +485,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
! by dt, in H m2 s-1.
h_frac ! The fraction of the mass in the column above the bottom
! interface of a layer that is within a layer, ND. 0<h_frac<=1
real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below, nondim)
hN2_y_PE ! thickness in m times Brunt-Vaisala freqeuncy at v-points (m s-2),
! used for calculating PE release
real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below, nondim)
hN2_x_PE ! thickness in m times Brunt-Vaisala freqeuncy at u-points (m s-2)
! used for calculating PE release
real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
pres, & ! The pressure at an interface, in Pa.
h_avail_rsum ! The running sum of h_avail above an interface, in H m2 s-1.
Expand All @@ -504,6 +513,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
real :: Work_u(SZIB_(G), SZJ_(G)) ! The work being done by the thickness
real :: Work_v(SZI_(G), SZJB_(G)) ! diffusion integrated over a cell, in W.
real :: Work_h ! The work averaged over an h-cell in W m-2.
real :: PE_release_h ! The amount of potential energy released by GM, averaged over an h-cell in m3 s-3.
! The calculation is equal to h * S^2 * N^2 * kappa_GM.
real :: I4dt ! 1 / 4 dt in s-1.
real :: drdiA, drdiB ! Along layer zonal- and meridional- potential density
real :: drdjA, drdjB ! gradients in the layers above (A) and below(B) the
Expand Down Expand Up @@ -576,6 +587,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

nk_linear = max(GV%nkml, 1)

Slope_x_PE(:,:,:) = 0.0
Slope_y_PE(:,:,:) = 0.0
hN2_x_PE(:,:,:) = 0.0
hN2_y_PE(:,:,:) = 0.0

find_work = .false.
if (associated(MEKE)) find_work = associated(MEKE%GM_src)
find_work = (associated(CS%GMwork) .or. find_work)
Expand Down Expand Up @@ -682,7 +698,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

if (k > nk_linear) then
if (use_EOS) then
if (CS%use_FGNV_streamfn .or. .not.present_slope_x) then
if (CS%use_FGNV_streamfn .or. find_work .or. .not.present_slope_x) then
hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
hg2R = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
Expand Down Expand Up @@ -740,6 +756,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
int_slope_u(I,j,K) * GV%Z_to_m*((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j))
slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K)
endif

Slope_x_PE(I,j,k) = MIN(Slope,CS%slope_max)
hN2_x_PE(I,j,k) = hN2_u(I,K) * GV%m_to_Z
if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope

! Estimate the streamfunction at each interface (m3 s-1).
Expand Down Expand Up @@ -928,7 +947,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

if (k > nk_linear) then
if (use_EOS) then
if (CS%use_FGNV_streamfn .or. .not.present_slope_y) then
if (CS%use_FGNV_streamfn .or. find_work .or. .not. present_slope_y) then
hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
hg2R = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
Expand Down Expand Up @@ -986,6 +1005,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
int_slope_v(i,J,K) * GV%Z_to_m*((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J))
slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K)
endif

Slope_y_PE(i,J,k) = MIN(Slope,CS%slope_max)
hN2_y_PE(i,J,k) = hN2_v(i,K) * GV%m_to_Z
if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope

! Estimate the streamfunction at each interface (m3 s-1).
Expand Down Expand Up @@ -1174,15 +1196,24 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
enddo
endif

if (find_work) then ; do j=js,je ; do i=is,ie

if (find_work) then ; do j=js,je ; do i=is,ie ; do k=nz,1,-1
! Note that the units of Work_v and Work_u are W, while Work_h is W m-2.
Work_h = 0.5 * G%IareaT(i,j) * &
((Work_u(I-1,j) + Work_u(I,j)) + (Work_v(i,J-1) + Work_v(i,J)))
PE_release_h = -0.25*(Kh_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + &
Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k) + &
Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + &
Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k))
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) + MIN(0.0,Work_h)
if (MEKE%GM_src_alt) then
MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h
else
MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h
endif
endif ; endif
enddo ; enddo ; endif
enddo ; enddo ; enddo ; endif

if (CS%id_slope_x > 0) call post_data(CS%id_slope_x, CS%diagSlopeX, CS%diag)
if (CS%id_slope_y > 0) call post_data(CS%id_slope_y, CS%diagSlopeY, CS%diag)
Expand Down Expand Up @@ -1899,11 +1930,11 @@ subroutine thickness_diffuse_get_KH(CS, KH_u_GME, KH_v_GME, G)
! Local variables
integer :: i,j,k

do k=1,G%ke ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
do k=1,G%ke+1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
KH_u_GME(I,j,k) = CS%KH_u_GME(I,j,k)
enddo ; enddo ; enddo

do k=1,G%ke ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
do k=1,G%ke+1 ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
KH_v_GME(i,J,k) = CS%KH_v_GME(i,J,k)
enddo ; enddo ; enddo

Expand Down

0 comments on commit c24bfb3

Please sign in to comment.