diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index f5a2b6c1f9..8106ebd8bf 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -11,6 +11,7 @@ module MOM_MEKE use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_domains, only : create_group_pass, do_group_pass use MOM_domains, only : group_pass_type +use MOM_domains, only : pass_var, pass_vector use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, MOM_mesg use MOM_file_parser, only : read_param, get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -43,6 +44,8 @@ module MOM_MEKE real :: MEKE_Ct !< Coefficient in the \f$\gamma_{bt}\f$ expression [nondim] 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 :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC + !! framework (Marshall et al., 2012) 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 @@ -57,8 +60,11 @@ module MOM_MEKE real :: MEKE_K4 !< Background bi-harmonic diffusivity (of MEKE) [m4 s-1] real :: KhMEKE_Fac !< A factor relating MEKE%Kh to the diffusivity used for !! MEKE itself [nondim]. - real :: viscosity_coeff !< The scaling coefficient in the expression for - !! viscosity used to parameterize lateral momentum mixing + real :: viscosity_coeff_Ku !< The scaling coefficient in the expression for + !! viscosity used to parameterize lateral harmonic momentum mixing + !! by unresolved eddies represented by MEKE. + real :: viscosity_coeff_Au !< The scaling coefficient in the expression for + !! viscosity used to parameterize lateral biharmonic momentum mixing !! by unresolved eddies represented by MEKE. real :: Lfixed !< Fixed mixing length scale [m]. real :: aDeform !< Weighting towards deformation scale of mixing length [nondim] @@ -89,6 +95,7 @@ module MOM_MEKE integer :: id_clock_pass !< Clock for group pass calls 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_Kh_diff !< 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 @@ -104,8 +111,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. - real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [s-1]. - real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [s-1]. + real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: SN_u !< Eady growth rate at u-points [s-1]. + real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: SN_v !< Eady growth rate at v-points [s-1]. type(vertvisc_type), intent(in) :: visc !< The vertical viscosity type. real, intent(in) :: dt !< Model(baroclinic) time-step [s]. type(MEKE_CS), pointer :: CS !< MEKE control structure. @@ -126,6 +133,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h 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 [kg m2 s-3]. Kh_u, & ! The zonal diffusivity that is actually used [m2 s-1]. @@ -425,6 +433,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Limit Kh to avoid CFL violations. if (associated(MEKE%Kh)) & Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh(i,j)+MEKE%Kh(i+1,j)) + if (associated(MEKE%Kh_diff)) & + Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i+1,j)) Inv_Kh_max = 2.0*sdt * ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * & max(G%IareaT(i,j),G%IareaT(i+1,j))) if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max @@ -438,6 +448,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h do J=js-1,je ; do i=is,ie if (associated(MEKE%Kh)) & Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh(i,j)+MEKE%Kh(i,j+1)) + if (associated(MEKE%Kh_diff)) & + Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i,j+1)) Inv_Kh_max = 2.0*sdt * ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * & max(G%IareaT(i,j),G%IareaT(i,j+1))) if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max @@ -491,7 +503,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h 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%MEKE(i,j) = MEKE%MEKE(i,j) -sdt_damp*drag_rate(i,j) MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) enddo ; enddo else @@ -508,6 +520,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP do endif endif ! MEKE_KH>=0 + + do j=js,je ; do i=is,ie + MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0) + enddo ; enddo + !$OMP end parallel call cpu_clock_begin(CS%id_clock_pass) call do_group_pass(CS%pass_MEKE, G%Domain) @@ -515,43 +532,53 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Calculate diffusivity for main model to use if (CS%MEKE_KhCoeff>0.) then - if (CS%use_old_lscale) then - if (CS%Rd_as_max_scale) then - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - MEKE%Kh(i,j) = (CS%MEKE_KhCoeff & + if (.not.CS%MEKE_GEOMETRIC) then + if (CS%use_old_lscale) then + if (CS%Rd_as_max_scale) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE%Kh(i,j) = (CS%MEKE_KhCoeff & * sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j))) & * min(MEKE%Rd_dx_h(i,j), 1.0) - enddo ; enddo + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE%Kh(i,j) = CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j)) + enddo ; enddo + endif else !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - MEKE%Kh(i,j) = CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j)) + MEKE%Kh(i,j) = (CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j)))*LmixScale(i,j)) enddo ; enddo endif - else - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - MEKE%Kh(i,j) = (CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j)))*LmixScale(i,j)) - enddo ; enddo - endif - call cpu_clock_begin(CS%id_clock_pass) - call do_group_pass(CS%pass_Kh, G%Domain) - call cpu_clock_end(CS%id_clock_pass) + call cpu_clock_begin(CS%id_clock_pass) + call do_group_pass(CS%pass_Kh, G%Domain) + call cpu_clock_end(CS%id_clock_pass) + endif endif ! Calculate viscosity for the main model to use - if (CS%viscosity_coeff/=0.) then + if (CS%viscosity_coeff_Ku /=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 + MEKE%Ku(i,j) = CS%viscosity_coeff_Ku*sqrt(2.*max(0.,MEKE%MEKE(i,j)))*LmixScale(i,j) enddo ; enddo call cpu_clock_begin(CS%id_clock_pass) call do_group_pass(CS%pass_Ku, G%Domain) + call cpu_clock_end(CS%id_clock_pass) + endif + + if (CS%viscosity_coeff_Au /=0.) then + do j=js,je ; do i=is,ie + MEKE%Au(i,j) = CS%viscosity_coeff_Au*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_Au, G%Domain) call cpu_clock_end(CS%id_clock_pass) endif + ! Offer fields for averaging. if (CS%id_MEKE>0) call post_data(CS%id_MEKE, MEKE%MEKE, CS%diag) if (CS%id_Ue>0) call post_data(CS%id_Ue, sqrt(max(0.,2.0*MEKE%MEKE)), CS%diag) @@ -897,6 +924,9 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) "into MEKE by the thickness mixing parameterization. \n"//& "If MEKE_GMCOEFF is negative, this conversion is not \n"//& "used or calculated.", units="nondim", default=-1.0) + call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & + "If MEKE_GEOMETRIC is true, uses the GM coefficient formulation \n"//& + "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.) call get_param(param_file, mdl, "MEKE_FRCOEFF", CS%MEKE_FrCoeff, & "The efficiency of the conversion of mean energy into \n"//& "MEKE. If MEKE_FRCOEFF is negative, this conversion \n"//& @@ -955,9 +985,15 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) "If true, the length scale used by MEKE is the minimum of\n"//& "the deformation radius or grid-spacing. Only used if\n"//& "MEKE_OLD_LSCALE=True", units="nondim", default=.false.) - call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF", CS%viscosity_coeff, & + call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_KU", CS%viscosity_coeff_Ku, & "If non-zero, is the scaling coefficient in the expression for\n"//& - "viscosity used to parameterize lateral momentum mixing by\n"//& + "viscosity used to parameterize harmonic lateral momentum mixing by\n"//& + "unresolved eddies represented by MEKE. Can be negative to\n"//& + "represent backscatter from the unresolved eddies.", & + units="nondim", default=0.0) + call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_AU", CS%viscosity_coeff_Au, & + "If non-zero, is the scaling coefficient in the expression for\n"//& + "viscosity used to parameterize biharmonic lateral momentum mixing by\n"//& "unresolved eddies represented by MEKE. Can be negative to\n"//& "represent backscatter from the unresolved eddies.", & units="nondim", default=0.0) @@ -989,7 +1025,7 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) "If true, initialize EKE to zero. Otherwise a local equilibrium solution\n"//& "is used as an initial condition for EKE.", default=.false.) call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", MEKE%backscatter_Ro_c, & - "The coefficient in the Rossby number function for scaling the buharmonic\n"//& + "The coefficient in the Rossby number function for scaling the biharmonic\n"//& "frictional energy source. Setting to non-zero enables the Rossby number function.", & units="nondim", default=0.0) call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_POW", MEKE%backscatter_Ro_pow, & @@ -1014,8 +1050,11 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) 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_Ku/=0. .and. .not. laplacian) call MOM_error(FATAL, & + "LAPLACIAN must be true if MEKE_VISCOSITY_COEFF_KU is true.") + + if (CS%viscosity_coeff_Au/=0. .and. .not. biharmonic) call MOM_error(FATAL, & + "BIHARMONIC must be true if MEKE_VISCOSITY_COEFF_AU is true.") call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.) @@ -1033,6 +1072,10 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) call create_group_pass(CS%pass_Kh, MEKE%Kh, G%Domain) call do_group_pass(CS%pass_Kh, G%Domain) endif + if (associated(MEKE%Kh_diff)) then + call create_group_pass(CS%pass_Kh_diff, MEKE%Kh_diff, G%Domain) + call do_group_pass(CS%pass_Kh_diff, G%Domain) + endif if (associated(MEKE%Ku)) then call create_group_pass(CS%pass_Ku, MEKE%Ku, G%Domain) call do_group_pass(CS%pass_Ku, G%Domain) @@ -1118,7 +1161,8 @@ 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_GMECoeff, MEKE_KHCoeff, MEKE_viscCoeff + real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff, MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au + logical :: Use_KH_in_MEKE logical :: useMEKE integer :: isd, ied, jsd, jed @@ -1130,8 +1174,9 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) 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) - + MEKE_viscCoeff_Ku =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku) + MEKE_viscCoeff_Au =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au) + Use_KH_in_MEKE = .false.; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE) ! Allocate control structure if (associated(MEKE)) then call MOM_error(WARNING, "MEKE_alloc_register_restart called with an associated "// & @@ -1150,8 +1195,8 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) call register_restart_field(MEKE%MEKE, vd, .false., restart_CS) if (MEKE_GMcoeff>=0.) then allocate(MEKE%GM_src(isd:ied,jsd:jed)) ; MEKE%GM_src(:,:) = 0.0 - endif - if (MEKE_FrCoeff>=0.) then + endif + if (MEKE_FrCoeff>=0. .or. MEKE_GMECoeff>=0.) then allocate(MEKE%mom_src(isd:ied,jsd:jed)) ; MEKE%mom_src(:,:) = 0.0 endif if (MEKE_GMECoeff>=0.) then @@ -1164,12 +1209,20 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) call register_restart_field(MEKE%Kh, vd, .false., restart_CS) endif allocate(MEKE%Rd_dx_h(isd:ied,jsd:jed)) ; MEKE%Rd_dx_h(:,:) = 0.0 - if (MEKE_viscCoeff/=0.) then + if (MEKE_viscCoeff_Ku/=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', & longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy") call register_restart_field(MEKE%Ku, vd, .false., restart_CS) + endif + if (Use_Kh_in_MEKE) then + allocate(MEKE%Kh_diff(isd:ied,jsd:jed)) ; MEKE%Kh_diff(:,:) = 0.0 + vd = var_desc("MEKE_Kh_diff", "m2 s-1",hor_grid='h',z_grid='1', & + longname="Copy of thickness diffusivity for diffusing MEKE") + call register_restart_field(MEKE%Kh_diff, vd, .false., restart_CS) + endif + if (MEKE_viscCoeff_Au/=0.) then 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") @@ -1193,6 +1246,7 @@ subroutine MEKE_end(MEKE, CS) 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%Kh_diff)) deallocate(MEKE%Kh_diff) if (associated(MEKE%Ku)) deallocate(MEKE%Ku) if (associated(MEKE%Au)) deallocate(MEKE%Au) if (allocated(CS%del2MEKE)) deallocate(CS%del2MEKE) diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90 index 2a5bd3a92f..0f4c58b68c 100644 --- a/src/parameterizations/lateral/MOM_MEKE_types.F90 +++ b/src/parameterizations/lateral/MOM_MEKE_types.F90 @@ -13,6 +13,7 @@ module MOM_MEKE_types 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. + Kh_diff => NULL(), & !< Uses the non-MEKE-derived thickness diffusion coefficient to diffuse MEKE, 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. real, dimension(:,:), pointer :: Ku => NULL() !< The MEKE-derived lateral viscosity coefficient [m2 s-1]. diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index f84df14b69..2a2625c599 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -167,7 +167,7 @@ module MOM_hor_visc 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 - integer :: id_FrictWorkMax = -1, id_target_FrictWork_GME = -1 + integer :: id_FrictWorkMax = -1 integer :: id_FrictWork_diss = -1, id_FrictWork_GME !!@} @@ -287,12 +287,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, 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) + target_diss_rate_GME, & ! the maximum theoretical dissipation plus the amount spuriously dissipated + ! by friction (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) div_xx_h ! horizontal divergence (s-1) !real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: & real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & @@ -350,6 +350,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, inv_PI6 = inv_PI3**2 epsilon = 1.e-7 + Ah_h(:,:,:) = 0.0 + Kh_h(:,:,:) = 0.0 + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally apply_OBC = .true. @@ -394,10 +397,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, !$OMP div_xx, div_xx_dx, div_xx_dy,local_strain, & !$OMP Shear_mag, h2uq, h2vq, Kh_scale, hrat_min) + do j=js,je ; do i=is,ie + boundary_mask(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1)) + enddo ; enddo + if (CS%use_GME) then ! GME tapers off above this depth H0 = 1000.0 - FWfrac = 0.1 + FWfrac = 1.0 GME_coeff_limiter = 1e7 ! initialize diag. array with zeros @@ -406,11 +413,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, str_xx_GME(:,:) = 0.0 str_xy_GME(:,:) = 0.0 - do j=js,je ; do i=is,ie - boundary_mask(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1)) - enddo ; enddo - - call pass_var(boundary_mask, G%Domain, complete=.true.) +! call pass_var(boundary_mask, G%Domain, complete=.true.) ! Get barotropic velocities and their gradients call barotropic_get_tav(Barotropic, ubtav, vbtav, G) @@ -472,12 +475,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2) enddo ; enddo - - ! 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 @@ -812,7 +809,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then - !vert_vort_mag = MIN(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j), beta_h(i,j)*3) vert_vort_mag = MIN(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j),3*grad_vort_mag_h_2d(i,j)) else vert_vort_mag = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j) @@ -855,7 +851,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 @@ -905,8 +900,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 @@ -951,15 +944,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then vert_vort_mag = MIN(grad_vort_mag_q(I,J) + grad_div_mag_q(I,J), 3*grad_vort_mag_q_2d(I,J)) - !vert_vort_mag = MIN(grad_vort_mag_q(I,J) + grad_div_mag_q(I,J), beta_q(I,J)*3) else vert_vort_mag = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J) endif endif h2uq = 4.0 * h_u(I,j) * h_u(I,j+1) h2vq = 4.0 * h_v(i,J) * h_v(i+1,J) - !hq = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * & - ! ((h(i,j,k) + h(i+1,j+1,k)) + (h(i,j+1,k) + h(i+1,j,k)))) hq(I,J) = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * & ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J)))) @@ -1110,14 +1100,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! 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) - & Ah_h(i,j,k) * grad_d2vel_mag_h(i,j) - 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 + max_diss_rate(i,j,k) = 2.0*MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j)) + FrictWork_diss(i,j,k) = diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 + 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 @@ -1125,9 +1114,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! 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 + else + FrictWork_diss(i,j,k) = diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 endif ; endif enddo ; enddo @@ -1139,13 +1129,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, if (.not. (associated(MEKE))) call MOM_error(FATAL, & "MEKE must be enabled for GME to be used.") - if (.not. (associated(MEKE%mom_src))) call MOM_error(FATAL, & - "MEKE%mom_src must be enabled for GME to be used.") - 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 = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_h(i,j) +! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*target_diss_rate_GME(i,j,k) / grad_vel_mag_bt_h(i,j) else GME_coeff = 0.0 endif @@ -1153,7 +1141,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! apply mask GME_coeff = GME_coeff * boundary_mask(i,j) - GME_coeff = MIN(GME_coeff,GME_coeff_limiter) + 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 @@ -1166,6 +1154,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, 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) +! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*target_diss_rate_GME(i,j,k) / grad_vel_mag_bt_q(I,J) else GME_coeff = 0.0 endif @@ -1173,7 +1162,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! apply mask GME_coeff = GME_coeff * boundary_mask(i,j) - GME_coeff = MIN(GME_coeff,GME_coeff_limiter) + GME_coeff = MIN(GME_coeff, GME_coeff_limiter) if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J) @@ -1325,31 +1314,20 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, enddo ; enddo 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 now is sign definite because it only uses the dissipation + MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + MAX(FrictWork_diss(i,j,k), FrictWorkMax(i,j,k)) enddo ; enddo + endif ! MEKE%backscatter - if (CS%use_GME) then - if (associated(MEKE%GME_snk)) then - do j=js,je ; do i=is,ie - MEKE%GME_snk(i,j) = MEKE%GME_snk(i,j) + FrictWork_GME(i,j,k) - enddo ; enddo - endif + if (CS%use_GME .and. associated(MEKE)) then + if (associated(MEKE%GME_snk)) then + do j=js,je ; do i=is,ie + MEKE%GME_snk(i,j) = MEKE%GME_snk(i,j) + FrictWork_GME(i,j,k) + enddo ; enddo endif - -! 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)) -! 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 -! endif - endif - endif ; endif + + endif ; endif ! find_FrictWork and associated(mom_src) enddo ! end of k loop @@ -1360,7 +1338,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, if (CS%id_FrictWorkMax>0) call post_data(CS%id_FrictWorkMax, FrictWorkMax, CS%diag) if (CS%id_FrictWork_diss>0) call post_data(CS%id_FrictWork_diss, FrictWork_diss, CS%diag) if (CS%id_FrictWork_GME>0) call post_data(CS%id_FrictWork_GME, FrictWork_GME, CS%diag) - if (CS%id_target_FrictWork_GME>0) call post_data(CS%id_target_FrictWork_GME, target_FrictWork_GME, CS%diag) if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag) if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag) if (CS%id_vort_xy_q>0) call post_data(CS%id_vort_xy_q, vort_xy_q, CS%diag) @@ -1583,7 +1560,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS) "the grid spacing to calculate the biharmonic viscosity. \n"//& "The final viscosity is the largest of this scaled \n"//& "viscosity, the Smagorinsky and Leith viscosities, and AH.", & - units="m s-1", default=0.0) + units="m s-1", default=0.1) call get_param(param_file, mdl, "SMAGORINSKY_AH", CS%Smagorinsky_Ah, & "If true, use a biharmonic Smagorinsky nonlinear eddy \n"//& "viscosity.", default=.false.) @@ -2055,9 +2032,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS) 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') - CS%id_target_FrictWork_GME = register_diag_field('ocean_model','target_FrictWork_GME',diag%axesTL,Time,& - 'Target for the amount of integral work done by lateral friction terms in GME', 'W m-2') - CS%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,Time,& 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', 'W m-2') endif diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index f7d008f366..52ca7a7b19 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -61,6 +61,11 @@ module MOM_thickness_diffuse logical :: debug !< write verbose checksums for debugging purposes logical :: use_GME_thickness_diffuse !< If true, passes GM coefficients to MOM_hor_visc for use !! with GME closure. + logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC + !! framework (Marshall et al., 2012) + real :: MEKE_GEOMETRIC_alpha !< The nondimensional coefficient governing the efficiency of + !! the GEOMETRIC thickness difussion [nondim] + logical :: Use_KH_in_MEKE !! If true, uses the thickness diffusivity calculated here to diffuse MEKE. logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather !! than the streamfunction for the GM source term. type(diag_ctrl), pointer :: diag => NULL() !< structure used to regulate timing of diagnostics @@ -139,6 +144,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp real :: hv(SZI_(G), SZJ_(G)) ! v-thickness [H ~> m or kg m-2] real :: KH_u_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [m2 s-1] real :: KH_v_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [m2 s-1] + real :: epsilon if (.not. associated(CS)) call MOM_error(FATAL, "MOM_thickness_diffuse:"// & "Module must be initialized before it is used.") @@ -170,6 +176,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp cg1 => null() endif + epsilon = 1.e-6 + !$OMP parallel do default(none) shared(is,ie,js,je,KH_u_CFL,dt,G,CS) do j=js,je ; do I=is-1,ie KH_u_CFL(I,j) = (0.25*CS%max_Khth_CFL) / & @@ -205,9 +213,16 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then !$OMP do - do j=js,je ; do I=is-1,ie - Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) - enddo ; enddo + if (CS%MEKE_GEOMETRIC) then + do j=js,je ; do I=is-1,ie + Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + G%mask2dCu(I,j) * CS%MEKE_GEOMETRIC_alpha * 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)) / & + (VarMix%SN_u(I,j) + epsilon) + enddo ; enddo + else + do j=js,je ; do I=is-1,ie + Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) + enddo ; enddo + endif endif ; endif if (Resoln_scaled) then @@ -276,9 +291,16 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then !$OMP do - do J=js-1,je ; do i=is,ie - Khth_Loc(i,j) = Khth_Loc(i,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) - enddo ; enddo + if (CS%MEKE_GEOMETRIC) then + do j=js-1,je ; do I=is,ie + Khth_Loc(I,j) = Khth_Loc(I,j) + G%mask2dCv(i,J) * CS%MEKE_GEOMETRIC_alpha * 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i,j+1)) / & + (VarMix%SN_v(i,J) + epsilon) + enddo ; enddo + else + do J=js-1,je ; do i=is,ie + Khth_Loc(i,j) = Khth_Loc(i,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) + enddo ; enddo + endif endif ; endif if (Resoln_scaled) then @@ -335,6 +357,17 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo ; enddo ; enddo endif + if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then +!$OMP do + if (CS%MEKE_GEOMETRIC) then + do j=js,je ; do I=is,ie + MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / & + (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) + enddo ; enddo + endif + endif ; endif + + !$OMP do do K=1,nz+1 ; do j=js,je ; do I=is-1,ie ; int_slope_u(I,j,K) = 0.0 ; enddo ; enddo ; enddo !$OMP do @@ -394,7 +427,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp ! in the case where KH_u and KH_v are depth independent. Otherwise, ! if use thickness weighted average, the variations of thickness with ! depth will place a spurious depth dependence to the diagnosed KH_t. - if (CS%id_KH_t > 0 .or. CS%id_KH_t1 > 0) then + if (CS%id_KH_t > 0 .or. CS%id_KH_t1 > 0 .or. CS%Use_KH_in_MEKE) then + MEKE%Kh_diff(:,:) = 0.0 do k=1,nz ! thicknesses across u and v faces, converted to 0/1 mask ! layer average of the interface diffusivities KH_u and KH_v @@ -413,9 +447,14 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp KH_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j)+hu(I,j)*KH_u_lay(I,j)) & +(hv(i,J-1)*KH_v_lay(i,J-1)+hv(i,J)*KH_v_lay(i,J))) & / (hu(I-1,j)+hu(I,j)+hv(i,J-1)+hv(i,J)+h_neglect) + MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) + KH_t(i,j,k) * h(i,j,k) enddo ; enddo enddo + do j=js,je ; do i=is,ie + MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j)) + enddo ; enddo + if (CS%id_KH_t > 0) call post_data(CS%id_KH_t, KH_t, CS%diag) if (CS%id_KH_t1 > 0) call post_data(CS%id_KH_t1, KH_t(:,:,1), CS%diag) endif @@ -1869,12 +1908,23 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) - call get_param(param_file, mdl, "USE_GME", CS%use_GME_thickness_diffuse, & - "If true, use the GM+E backscatter scheme in association \n"//& - "with the Gent and McWilliams parameterization.", 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 GM source term.", default=.false.) + call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & + "If true, uses the GM coefficient formulation \n"//& + "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.) + call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, & + "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//& + "thickness diffusion.", units="nondim", default=0.05) + call get_param(param_file, mdl, "USE_KH_IN_MEKE", CS%Use_KH_in_MEKE, & + "If true, uses the thickness diffusivity calculated here to diffuse \n"//& + "MEKE.", default=.false.) + + call get_param(param_file, mdl, "USE_GME", CS%use_GME_thickness_diffuse, & + "If true, use the GM+E backscatter scheme in association \n"//& + "with the Gent and McWilliams parameterization.", default=.false.) + if (CS%use_GME_thickness_diffuse) then call safe_alloc_ptr(CS%KH_u_GME,G%IsdB,G%IedB,G%jsd,G%jed,G%ke+1) call safe_alloc_ptr(CS%KH_v_GME,G%isd,G%ied,G%JsdB,G%JedB,G%ke+1)