From 80c65c72421ecbb9945021ed8600a7ac305256c8 Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Fri, 3 May 2019 10:35:37 -0600 Subject: [PATCH] Added a GME sink term to the MEKE budget. Cleaned up MOM_hor_visc.F90. --- src/parameterizations/lateral/MOM_MEKE.F90 | 27 ++- .../lateral/MOM_MEKE_types.F90 | 1 + .../lateral/MOM_hor_visc.F90 | 167 ++++-------------- 3 files changed, 61 insertions(+), 134 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index d44bf88c32..a02eed93af 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -30,6 +30,7 @@ module MOM_MEKE ! Parameters real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE (non-dim) real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE (non-dim) + real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME (non-dim) real :: MEKE_damping !< Local depth-independent MEKE dissipation rate in s-1. real :: MEKE_Cd_scale !< The ratio of the bottom eddy velocity to the column mean !! eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1 @@ -73,7 +74,7 @@ module MOM_MEKE !>@{ Diagnostic handles integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1 integer :: id_Ub = -1, id_Ut = -1 - integer :: id_GM_src = -1, id_mom_src = -1, id_decay = -1 + integer :: id_GM_src = -1, id_mom_src = -1, id_GME_snk = -1, id_decay = -1 integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1 integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1 integer :: id_Lrhines = -1, id_Leady = -1 @@ -113,6 +114,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv) MEKE_decay, & ! The MEKE decay timescale, in s-1. MEKE_GM_src, & ! The MEKE source from thickness mixing, in m2 s-3. MEKE_mom_src, & ! The MEKE source from momentum, in m2 s-3. + MEKE_GME_snk, & ! The MEKE sink from GME backscatter, in m2 s-3. drag_rate_visc, & drag_rate, & ! The MEKE spindown timescale due to bottom drag, in s-1. LmixScale, & ! Square of eddy mixing length, in m2. @@ -165,6 +167,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv) if (CS%debug) then if (associated(MEKE%mom_src)) call hchksum(MEKE%mom_src, 'MEKE mom_src',G%HI) + if (associated(MEKE%GME_snk)) call hchksum(MEKE%GME_snk, 'MEKE GME_snk',G%HI) if (associated(MEKE%GM_src)) call hchksum(MEKE%GM_src, 'MEKE GM_src',G%HI) if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE',G%HI) call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI) @@ -292,6 +295,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv) enddo ; enddo endif + if (associated(MEKE%GME_snk)) then +!$OMP do + do j=js,je ; do i=is,ie + src(i,j) = src(i,j) - CS%MEKE_GMECoeff*I_mass(i,j)*MEKE%GME_snk(i,j) + enddo ; enddo + endif + if (associated(MEKE%GM_src)) then !$OMP do if (CS%GM_src_alt) then @@ -578,6 +588,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv) if (CS%id_decay>0) call post_data(CS%id_decay, MEKE_decay, CS%diag) if (CS%id_GM_src>0) call post_data(CS%id_GM_src, MEKE%GM_src, CS%diag) if (CS%id_mom_src>0) call post_data(CS%id_mom_src, MEKE%mom_src, CS%diag) + if (CS%id_GME_snk>0) call post_data(CS%id_GME_snk, MEKE%GME_snk, CS%diag) if (CS%id_Le>0) call post_data(CS%id_Le, LmixScale, CS%diag) if (CS%id_gamma_b>0) then do j=js,je ; do i=is,ie @@ -892,6 +903,10 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) "The efficiency of the conversion of mean energy into \n"//& "MEKE. If MEKE_FRCOEFF is negative, this conversion \n"//& "is not used or calculated.", units="nondim", default=-1.0) + call get_param(param_file, mdl, "MEKE_GMECOEFF", CS%MEKE_GMECoeff, & + "The efficiency of the conversion of MEKE into mean energy \n"//& + "by GME. If MEKE_GMECOEFF is negative, this conversion \n"//& + "is not used or calculated.", units="nondim", default=-1.0) call get_param(param_file, mdl, "MEKE_BGSRC", CS%MEKE_BGsrc, & "A background energy source for MEKE.", units="W kg-1", & default=0.0) @@ -1065,6 +1080,9 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) CS%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, Time, & 'MEKE energy available from momentum', 'W m-2') if (.not. associated(MEKE%mom_src)) CS%id_mom_src = -1 + CS%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, Time, & + 'MEKE energy lost to GME backscatter', 'W m-2') + if (.not. associated(MEKE%GME_snk)) CS%id_GME_snk = -1 CS%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, Time, & 'Eddy mixing length used in the MEKE derived eddy diffusivity', 'm') CS%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, Time, & @@ -1097,7 +1115,7 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) type(MOM_restart_CS), pointer :: restart_CS !< Restart control structure for MOM_MEKE. ! Local variables type(vardesc) :: vd - real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_KHCoeff, MEKE_viscCoeff + real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff, MEKE_KHCoeff, MEKE_viscCoeff logical :: useMEKE integer :: isd, ied, jsd, jed @@ -1107,6 +1125,7 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) ! Read these parameters to determine what should be in the restarts MEKE_GMcoeff =-1.; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff) MEKE_FrCoeff =-1.; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff) + MEKE_GMEcoeff =-1.; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff) MEKE_KhCoeff =1.; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff) MEKE_viscCoeff =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF",MEKE_viscCoeff) @@ -1132,6 +1151,9 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) if (MEKE_FrCoeff>=0.) then allocate(MEKE%mom_src(isd:ied,jsd:jed)) ; MEKE%mom_src(:,:) = 0.0 endif + if (MEKE_GMECoeff>=0.) then + allocate(MEKE%GME_snk(isd:ied,jsd:jed)) ; MEKE%GME_snk(:,:) = 0.0 + endif if (MEKE_KhCoeff>=0.) then allocate(MEKE%Kh(isd:ied,jsd:jed)) ; MEKE%Kh(:,:) = 0.0 vd = var_desc("MEKE_Kh", "m2 s-1",hor_grid='h',z_grid='1', & @@ -1166,6 +1188,7 @@ subroutine MEKE_end(MEKE, CS) if (associated(MEKE%MEKE)) deallocate(MEKE%MEKE) if (associated(MEKE%GM_src)) deallocate(MEKE%GM_src) if (associated(MEKE%mom_src)) deallocate(MEKE%mom_src) + 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) diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90 index 26257360b5..e5d0ce9072 100644 --- a/src/parameterizations/lateral/MOM_MEKE_types.F90 +++ b/src/parameterizations/lateral/MOM_MEKE_types.F90 @@ -11,6 +11,7 @@ module MOM_MEKE_types MEKE => NULL(), & !< Vertically averaged eddy kinetic energy, in m2 s-2. GM_src => NULL(), & !< MEKE source due to thickness mixing (GM), in W m-2. mom_src => NULL(),& !< MEKE source from lateral friction in the momentum equations, in W m-2. + GME_snk => NULL(),& !< MEKE sink from GME backscatter in the momentum equations, in W m-2. Kh => NULL(), & !< The MEKE-derived lateral mixing coefficient in m2 s-1. Rd_dx_h => NULL() !< The deformation radius compared with the grid spacing, nondim. !! Rd_dx_h is copied from VarMix_CS. diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 94094759ed..a8e5005205 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -76,7 +76,7 @@ module MOM_hor_visc logical :: dynamic_aniso !< If true, the anisotropic viscosity is recomputed as a function !! of state. This is set depending on ANISOTROPIC_MODE. logical :: use_GME !< If true, use GME backscatter scheme. - real :: GME_dt + real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Kh_bg_xx !< The background Laplacian viscosity at h points, in units !! of m2 s-1. The actual viscosity may be the larger of this @@ -326,11 +326,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, real :: local_strain ! Local variable for interpolating computed strain rates (s-1). real :: epsilon real :: GME_coeff ! The GME (negative) viscosity coefficient (m2 s-1) - real :: GME_coeff_limiter - real :: FWfrac + real :: GME_coeff_limiter ! Maximum permitted value of the GME coefficient (m2 s-1) + real :: FWfrac ! Fraction of maximum theoretical energy transfer to use when scaling GME coefficient real :: DY_dxBu, DX_dyBu - real :: H0 - real :: dr, dv + real :: H0 ! Depth used to scale down GME coefficient in shallow areas (m) logical :: rescale_Kh, legacy_bound logical :: find_FrictWork logical :: apply_OBC = .false. @@ -396,6 +395,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! GME tapers off above this depth H0 = 1000.0 FWfrac = 0.1 + GME_coeff_limiter = 1e7 + ! initialize diag. array with zeros GME_coeff_h(:,:,:) = 0.0 GME_coeff_q(:,:,:) = 0.0 @@ -1091,7 +1092,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, 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)) + 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 @@ -1111,64 +1112,24 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, if (CS%use_GME) then + + if (.not. (associated(MEKE))) call MOM_error(FATAL, & + "MEKE must be enabled for GME to be used.") - do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (.not. (associated(MEKE%mom_src))) call MOM_error(FATAL, & + "MEKE%mom_src must be enabled for GME to be used.") -! 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)) * & -! (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) - - if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then -! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_FrictWork_GME(i,j,k) / ( 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 = (MIN(G%bathyT(i,j)/H0,1.0)**2) * max_diss_rate(i,j,k) * (G%areaT(i,j)/0.01) - -! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * MEKE%MEKE(i,j) / MAX(0.25*(VarMix%SN_u(I,j)+VarMIX%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)),epsilon) - GME_coeff = 1e10; dr = 1e6 - -! if (G%mask2dT(i,j)) -! do while (dr >= max_diss_rate(i,j,k)) -! GME_coeff = GME_coeff / 2.0 - dr = GME_coeff * grad_vel_mag_h(i,j) + do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_h(i,j)>0) ) then - GME_coeff = FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_h(i,j) -! GME_coeff = max_diss_rate_bt(i,j) / grad_vel_mag_bt_h(i,j) - else - GME_coeff = 0.0 - endif - -! GME_coeff = (MIN(G%bathyT(i,j)*G%Zd_to_m/H0,1.0)**2)*GME_coeff - -! GME_coeff = 10.0 -! enddo -! else -! GME_coeff = 0.0 -! endif - - - endif ; endif + if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_h(i,j)>0) ) then + GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_h(i,j) + else + GME_coeff = 0.0 + endif ! 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 = 1e7 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 @@ -1180,74 +1141,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, do J=js-1,Jeq ; do I=is-1,Ieq -! 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)) * & -! ( 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 ) / & -! ( 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 * (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) - - if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then -! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_FrictWork_GME(i,j,k) / (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 = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_diss_rate_GME(i,j,k) * G%areaT(i,j) / & -! (0.1**2) -! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_diss_rate_GME(i,j,k) * sqrt(G%areaT(i,j))/0.0001 -! -! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * MEKE%MEKE(i,j) / MAX(0.25*(VarMix%SN_u(I,j)+VarMIX%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)),epsilon) -! - -! if (G%mask2dT(i,j)) - GME_coeff = 1e10; dr = 1e6 - - if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_q(i,j)>0) ) then - GME_coeff = FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_q(I,J) -! GME_coeff = max_diss_rate_bt(i,j) / grad_vel_mag_bt_q(I,J) - else - GME_coeff = 0.0 - endif - -! dr = GME_coeff * bt_grad_mag_q(I,J) -! if ((max_diss_rate(i,j,k) > 0) .and. (dr>0) ) then -! dv = dr / (FWfrac*max_diss_rate(i,j,k)) -! GME_coeff = GME_coeff / dv -! else -! GME_coeff = 0.0 -! endif + if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_q(i,j)>0) ) then + GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_q(I,J) + else + GME_coeff = 0.0 + endif -! GME_coeff = (MIN(G%bathyT(i,j)*G%Zd_to_m/H0,1.0)**2)*GME_coeff -! do while (dr >= max_diss_rate(i,j,k)) -! GME_coeff = GME_coeff / 2.0 -! dr = GME_coeff * (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) -! enddo -! else -! GME_coeff = 0.0 -! endif -! GME_coeff = 10.0 - - endif ; endif - ! 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 = 1e7 GME_coeff = MIN(GME_coeff,GME_coeff_limiter) if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff @@ -1272,11 +1174,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, endif enddo ; enddo - if (find_FrictWork) 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_kg_m2 * grad_vel_mag_bt_h(i,j) - enddo ; enddo - endif + 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_kg_m2 * grad_vel_mag_bt_h(i,j) + enddo ; enddo + endif else ! use_GME do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1401,11 +1303,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! 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)) enddo ; enddo + if (CS%use_GME) then - 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) + FrictWork_GME(i,j,k) - enddo ; enddo + if (associated(MEKE%GME_snk)) then + do j=js,je ; do i=is,ie + ! MEKE%mom_src now is sign definite because it only uses the dissipation + MEKE%GME_snk(i,j) = MEKE%GME_snk(i,j) + FrictWork_GME(i,j,k) + enddo ; enddo + endif endif endif endif ; endif @@ -1864,8 +1769,6 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) CS%reduction_xy(I,J) = G%dx_Cv(i+1,J) / G%dxCv(i+1,J) enddo ; enddo - CS%GME_dt = dt - if (CS%Laplacian) then ! The 0.3 below was 0.4 in MOM1.10. The change in hq requires ! this to be less than 1/3, rather than 1/2 as before.