From 4b6ef52a94553024fabe56b97c39007f20e3ace7 Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Fri, 3 May 2019 15:06:04 -0600 Subject: [PATCH] Revert "Added biharmonic option for MEKE viscosity." This reverts commit 6b321744f269f725568634e68280150507f7633e. Conflicts: src/parameterizations/lateral/MOM_MEKE.F90 src/parameterizations/lateral/MOM_hor_visc.F90 --- src/parameterizations/lateral/MOM_MEKE.F90 | 30 ++----- .../lateral/MOM_MEKE_types.F90 | 1 - .../lateral/MOM_hor_visc.F90 | 78 ++++++------------- 3 files changed, 29 insertions(+), 80 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index a02eed93af..3882e2c974 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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. @@ -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.) @@ -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) @@ -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 @@ -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 @@ -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) diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90 index e5d0ce9072..d87d4b810f 100644 --- a/src/parameterizations/lateral/MOM_MEKE_types.F90 +++ b/src/parameterizations/lateral/MOM_MEKE_types.F90 @@ -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. diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index c6f9df0df0..b73c5227fa 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -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) @@ -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) @@ -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)) :: & @@ -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 @@ -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, & @@ -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 + & @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -1085,10 +1053,12 @@ 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 @@ -1096,15 +1066,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, 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 @@ -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 @@ -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