From d5f5d3203137a0cde9efc87d8ae84122fa9ab996 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 Jan 2019 11:26:43 -0700 Subject: [PATCH] Adds GME and smoothing function --- .../lateral/MOM_hor_visc.F90 | 224 +++++++++++++++++- 1 file changed, 211 insertions(+), 13 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 243d52abbc..f0f296a0a5 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -5,11 +5,13 @@ module MOM_hor_visc use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type -use MOM_domains, only : pass_var +use MOM_domains, only : pass_var, CORNER use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_Leith_viscosity +use MOM_barotropic, only : barotropic_CS, barotropic_get_tav +use MOM_thickness_diffuse, only : thickness_diffuse_CS, thickness_diffuse_get_KH use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE @@ -73,6 +75,7 @@ module MOM_hor_visc real :: Kh_aniso !< The anisotropic viscosity in m2 s-1. 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 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 @@ -159,6 +162,7 @@ module MOM_hor_visc integer :: id_diffu = -1, id_diffv = -1 integer :: id_Ah_h = -1, id_Ah_q = -1 integer :: id_Kh_h = -1, id_Kh_q = -1 + integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1 integer :: id_vort_xy_q = -1, id_div_xx_h = -1 integer :: id_FrictWork = -1, id_FrictWorkIntz = -1 !!@} @@ -179,7 +183,9 @@ module MOM_hor_visc !! u[is-2:ie+2,js-2:je+2] !! 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, CS, OBC) + +subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, & + thickness_diffuse, G, GV, CS, OBC) 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)), & @@ -199,24 +205,34 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, !! related to Mesoscale Eddy Kinetic Energy. type(VarMix_CS), pointer :: VarMix !< Pointer to a structure with fields that !! specify the spatially variable viscosities - type(hor_visc_CS), pointer :: CS !< Pontrol structure returned by a previous + type(barotropic_CS), pointer :: Barotropic !< Pointer to a structure containing + !! barotropic velocities + type(thickness_diffuse_CS), pointer :: thickness_diffuse !< Pointer to a structure containing + !! interface height diffusivities + type(hor_visc_CS), pointer :: CS !< Control structure returned by a previous !! call to hor_visc_init. type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type + ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & u0, & ! Laplacian of u (m-1 s-1) h_u, & ! Thickness interpolated to u points, in H. vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1) - div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1) + div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1) + ubtav ! zonal barotropic vel. ave. over baroclinic time-step (m s-1) real, dimension(SZI_(G),SZJB_(G)) :: & v0, & ! Laplacian of v (m-1 s-1) h_v, & ! Thickness interpolated to v points, in H. vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 s-1) - div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1) + div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1) + vbtav ! meridional barotropic vel. ave. over baroclinic time-step (m s-1) real, dimension(SZI_(G),SZJ_(G)) :: & + dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension (s-1) div_xx, & ! Estimate of horizontal divergence at h-points (s-1) sh_xx, & ! horizontal tension (du/dx - dv/dy) (1/sec) including metric terms + sh_xx_bt, & ! barotropic horizontal tension (du/dx - dv/dy) (1/sec) including metric terms str_xx,& ! str_xx is the diagonal term in the stress tensor (H m2 s-2) + str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME (H m2 s-2) bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2) FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction (W/m2) Leith_Kh_h, & ! Leith Laplacian viscosity at h-points (m2 s-1) @@ -227,8 +243,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain (s-1) + dvdx_bt, dudy_bt, & ! components in the barotropic shearing strain (s-1) sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) (1/sec) including metric terms + sh_xy_bt, & ! barotropic horizontal shearing strain (du/dy + dv/dx) (1/sec) inc. metric terms str_xy, & ! str_xy is the cross term in the stress tensor (H m2 s-2) + str_xy_GME, & ! smoothed cross term in the stress tensor from GME (H m2 s-2) bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2) vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (s-1) Leith_Kh_q, & ! Leith Laplacian viscosity at q-points (m2 s-1) @@ -238,15 +257,22 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, grad_div_mag_q ! Magnitude of divergence gradient at q-points (m-1 s-1) real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: & - Ah_q, & ! biharmonic viscosity at corner points (m4/s) - Kh_q, & ! Laplacian viscosity at corner points (m2/s) - vort_xy_q ! vertical vorticity at corner points (s-1) - + Ah_q, & ! biharmonic viscosity at corner points (m4/s) + Kh_q, & ! Laplacian viscosity at corner points (m2/s) + vort_xy_q, & ! vertical vorticity at corner points (s-1) + GME_coeff_q !< GME coeff. at q-points (m2 s-1) + + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & + KH_u_GME !< interface height diffusivities in u-columns (m2 s-1) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & + KH_v_GME !< interface height diffusivities in v-columns (m2 s-1) 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) FrictWork, & ! energy dissipated by lateral friction (W/m2) - div_xx_h ! horizontal divergence (s-1) + div_xx_h, & ! horizontal divergence (s-1) + KH_t_GME, & !< interface height diffusivities in t-columns (m2 s-1) + GME_coeff_h !< GME coeff. at h-points (m2 s-1) real :: Ah ! biharmonic viscosity (m4/s) real :: Kh ! Laplacian viscosity (m2/s) @@ -278,6 +304,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, real :: FatH ! abs(f) at h-point for MEKE source term (s-1) 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 :: DY_dxBu, DX_dyBu logical :: rescale_Kh, legacy_bound logical :: find_FrictWork @@ -330,12 +357,46 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, !$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, & !$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE, & !$OMP mod_Leith, legacy_bound, div_xx_h, vort_xy_q) & -!$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, & +!$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, & !$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, & +!$OMP sh_xx_bt, sh_xy_bt, dvdx_bt, dudy_bt, & !$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, & !$OMP vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, & !$OMP div_xx, div_xx_dx, div_xx_dy, & !$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min) + + + if (CS%use_GME) then + call barotropic_get_tav(Barotropic, ubtav, vbtav, G) + + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 + dudx_bt(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * ubtav(I,j) - & + G%IdyCu(I-1,j) * ubtav(I-1,j)) + dvdy_bt(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * vbtav(i,J) - & + G%IdxCv(i,J-1) * vbtav(i,J-1)) + sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j) + enddo ; enddo + + ! Components for the barotropic shearing strain + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + dvdx_bt(I,J) = CS%DY_dxBu(I,J)*(vbtav(i+1,J)*G%IdyCv(i+1,J) & + - vbtav(i,J)*G%IdyCv(i,J)) + dudy_bt(I,J) = CS%DX_dyBu(I,J)*(ubtav(I,j+1)*G%IdxCu(I,j+1) & + - ubtav(I,j)*G%IdxCu(I,j)) + enddo ; enddo + + if (CS%no_slip) then + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) ) + enddo ; enddo + else + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + sh_xy_bt(I,J) = G%mask2dBu(I,J) * ( dvdx_bt(I,J) + dudy_bt(I,J) ) + enddo ; enddo + endif + + endif ! use_GME + do k=1,nz ! The following are the forms of the horizontal tension and horizontal @@ -636,8 +697,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif ! CS%Leith_Kh - - do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + & @@ -691,6 +750,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, str_xx(i,j) = 0.0 endif ! Laplacian + if (CS%use_GME) then + call thickness_diffuse_get_KH(thickness_diffuse, KH_t_GME, KH_u_GME, KH_v_GME, G) + GME_coeff = KH_t_GME(i,j,k) * & + 0.5*(VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)) * & + ( (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) + + str_xx_GME(i,j) = GME_coeff * sh_xx_bt(i,j) + + endif + + if (CS%id_GME_coeff_h>0) GME_coeff_h(i,j,k) = GME_coeff + if (CS%anisotropic) then ! Shearing-strain averaged to h-points local_strain = 0.25 * ( (sh_xy(I,J) + sh_xy(I-1,J-1)) + (sh_xy(I-1,J) + sh_xy(I,J-1)) ) @@ -740,6 +816,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo + ! applying GME diagonal term + if (CS%use_GME) then + call smooth_GME(CS,G,GME_flux_h=str_xx_GME) + do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx(i,j) = str_xx(i,j) + str_xx_GME(i,j) + enddo ; enddo + endif + if (CS%biharmonic) then ! Gradient of Laplacian, for use in bi-harmonic term do J=js-1,Jeq ; do I=is-1,Ieq @@ -854,6 +938,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, str_xy(I,J) = 0.0 endif ! Laplacian + if (CS%use_GME) then + GME_coeff = ( 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*(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.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(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)) + + str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J) + endif + + if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff + if (CS%anisotropic) then ! Horizontal-tension averaged to q-points local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) ) @@ -902,6 +1003,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif enddo ; enddo + ! applying GME diagonal term + if (CS%use_GME) then + call smooth_GME(CS,G,GME_flux_q=str_xy_GME) + do J=js-1,Jeq ; do I=is-1,Ieq + if (CS%no_slip) then + str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq * CS%reduction_xy(I,J)) + else + str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) + endif + enddo ; enddo + endif + ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent. do j=js,je ; do I=Isq,Ieq diffu(I,j,k) = ((G%IdyCu(I,j)*(CS%DY2h(i,j) *str_xx(i,j) - & @@ -1021,6 +1134,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, if (CS%id_Ah_q>0) call post_data(CS%id_Ah_q, Ah_q, CS%diag) if (CS%id_Kh_h>0) call post_data(CS%id_Kh_h, Kh_h, CS%diag) if (CS%id_Kh_q>0) call post_data(CS%id_Kh_q, Kh_q, CS%diag) + if (CS%id_GME_coeff_h > 0) call post_data(CS%id_GME_coeff_h, GME_coeff_h, CS%diag) + if (CS%id_GME_coeff_q > 0) call post_data(CS%id_GME_coeff_q, GME_coeff_q, CS%diag) if (CS%id_FrictWorkIntz > 0) then do j=js,je @@ -1309,6 +1424,9 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) "viscosities. The final viscosity is the maximum of the other "//& "terms and this background value.", default=.false.) + call get_param(param_file, mdl, "USE_GME", CS%use_GME, & + "If true, use the GM+E backscatter scheme in association \n"//& + "with the Gent and McWilliams parameterization.", default=.false.) if (CS%bound_Kh .or. CS%bound_Ah .or. CS%better_bound_Kh .or. CS%better_bound_Ah) & call get_param(param_file, mdl, "DT", dt, & @@ -1696,6 +1814,14 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) endif + if (CS%use_GME) then + CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, & + 'GME coefficient at h Points', 'm^2 s-1') + + CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, & + 'GME coefficient at q Points', 'm^2 s-1') + endif + CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,& 'Integral work done by lateral friction terms', 'W m-2') @@ -1730,6 +1856,78 @@ subroutine align_aniso_tensor_to_grid(CS, n1, n2) end subroutine align_aniso_tensor_to_grid +!> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any +!! horizontal two-grid-point noise +subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) + ! Arguments + type(hor_visc_CS), pointer :: CS !< Control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid + real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: GME_flux_h!< + real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: GME_flux_q!< + + ! local variables + real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original + real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original + real :: wc, ww, we, wn, ws ! averaging weights for smoothing + integer :: i, j, k, s + + !do s=1,CS%n_smooth + do s=1,1 + + ! Update halos + if (present(GME_flux_h)) then + call pass_var(GME_flux_h, G%Domain) + GME_flux_h_original = GME_flux_h + ! apply smoothing on GME + do j = G%jsc, G%jec + do i = G%isc, G%iec + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + + ! compute weights + ww = 0.125 * G%mask2dT(i-1,j) + we = 0.125 * G%mask2dT(i+1,j) + ws = 0.125 * G%mask2dT(i,j-1) + wn = 0.125 * G%mask2dT(i,j+1) + wc = 1.0 - (ww+we+wn+ws) + + GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) & + + ww * GME_flux_h_original(i-1,j) & + + we * GME_flux_h_original(i+1,j) & + + ws * GME_flux_h_original(i,j-1) & + + wn * GME_flux_h_original(i,j+1) + enddo; enddo + endif + + ! Update halos + if (present(GME_flux_q)) then + call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.) + GME_flux_q_original = GME_flux_q + ! apply smoothing on GME + do J = G%JscB, G%JecB + do I = G%IscB, G%IecB + ! skip land points + if (G%mask2dBu(I,J)==0.) cycle + + ! compute weights + ww = 0.125 * G%mask2dBu(I-1,J) + we = 0.125 * G%mask2dBu(I+1,J) + ws = 0.125 * G%mask2dBu(I,J-1) + wn = 0.125 * G%mask2dBu(I,J+1) + wc = 1.0 - (ww+we+wn+ws) + + GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) & + + ww * GME_flux_q_original(I-1,J) & + + we * GME_flux_q_original(I+1,J) & + + ws * GME_flux_q_original(I,J-1) & + + wn * GME_flux_q_original(I,J+1) + enddo; enddo + endif + + enddo ! s-loop + +end subroutine smooth_GME + !> Deallocates any variables allocated in hor_visc_init. subroutine hor_visc_end(CS) type(hor_visc_CS), pointer :: CS !< The control structure returned by a