diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 156a397ff6..a434d50020 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1190,7 +1190,9 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, & real, dimension(:,:,:), pointer :: & u => NULL(), & ! u : zonal velocity component [L T-1 ~> m s-1] v => NULL(), & ! v : meridional velocity component [L T-1 ~> m s-1] - h => NULL() ! h : layer thickness [H ~> m or kg m-2] + h => NULL(), & ! h : layer thickness [H ~> m or kg m-2] + uh => NULL(), & ! uh : layer thickness times u [L2 H T-1 ~> m3 s-1 or kg s-1] + vh => NULL() ! vh : layer thickness times v [L2 H T-1 ~> m3 s-1 or kg s-1] logical :: calc_dtbt ! Indicates whether the dynamically adjusted ! barotropic time step needs to be updated. @@ -1204,7 +1206,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, & Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - u => CS%u ; v => CS%v ; h => CS%h + u => CS%u ; v => CS%v ; h => CS%h ; uh => CS%uh ; vh => CS%vh showCallTree = callTree_showQuery() call cpu_clock_begin(id_clock_dynamics) @@ -1225,7 +1227,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, & if (CS%thickness_diffuse) then call cpu_clock_begin(id_clock_thick_diff) if (CS%VarMix%use_variable_mixing) & - call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) + call calc_slope_functions(h, uh, vh, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt_tr_adv, G, GV, US, & CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, & CS%stoch_CS) @@ -1378,7 +1380,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, & if (CS%thickness_diffuse) then call cpu_clock_begin(id_clock_thick_diff) if (CS%VarMix%use_variable_mixing) & - call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) + call calc_slope_functions(h, uh, vh, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, US, & CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, CS%stoch_CS) @@ -2067,7 +2069,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call pass_var(CS%h, G%Domain) call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline) call calc_depth_function(G, CS%VarMix) - call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) + call calc_slope_functions(CS%h, CS%uh, CS%vh, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, CS%visc, G, GV, US, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) @@ -2094,7 +2096,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call pass_var(CS%h, G%Domain) call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline) call calc_depth_function(G, CS%VarMix) - call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) + call calc_slope_functions(CS%h, CS%uh, CS%vh, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, CS%visc, G, GV, US, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index f2f476b0c8..e42cb4c5d1 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -29,6 +29,7 @@ module MOM_lateral_mixing_coeffs type, public :: VarMix_CS logical :: initialized = .false. !< True if this control structure has been initialized. logical :: use_variable_mixing !< If true, use the variable mixing. + logical :: use_gradient_model !< If true, use the gradient model (Khani & Dawson, JAMES 2023). logical :: Resoln_scaling_used !< If true, a resolution function is used somewhere to scale !! away one of the viscosities or diffusivities when the !! deformation radius is well resolved. @@ -87,8 +88,12 @@ module MOM_lateral_mixing_coeffs real, allocatable :: SN_u(:,:) !< S*N at u-points [T-1 ~> s-1] real, allocatable :: SN_v(:,:) !< S*N at v-points [T-1 ~> s-1] + real, allocatable :: UH_grad(:,:,:) !< Grad model at u-points [T-1 ~> s-1] + real, allocatable :: VH_grad(:,:,:) !< Grad model at v-points [T-1 ~> s-1] real, allocatable :: L2u(:,:) !< Length scale^2 at u-points [L2 ~> m2] real, allocatable :: L2v(:,:) !< Length scale^2 at v-points [L2 ~> m2] + real, allocatable :: L2grad_u(:,:) !< Grad length scale^2 at u-points [L2 ~> m2] + real, allocatable :: L2grad_v(:,:) !< Grad length scale^2 at v-points [L2 ~> m2] real, allocatable :: cg1(:,:) !< The first baroclinic gravity wave speed [L T-1 ~> m s-1]. real, allocatable :: Res_fn_h(:,:) !< Non-dimensional function of the ratio the first baroclinic !! deformation radius to the grid spacing at h points [nondim]. @@ -150,6 +155,7 @@ module MOM_lateral_mixing_coeffs integer :: VarMix_Ktop !< Top layer to start downward integrals real :: Visbeck_L_scale !< Fixed length scale in Visbeck formula [L ~> m], or if negative a scaling !! factor [nondim] relating this length scale squared to the cell area + real :: Grad_Khani_Scale !< A scaling factor in Gradient formula [nondim] real :: Eady_GR_D_scale !< Depth over which to average SN [Z ~> m] real :: Res_coef_khth !< A coefficient [nondim] that determines the function !! of resolution, used for thickness and tracer mixing, as: @@ -176,6 +182,7 @@ module MOM_lateral_mixing_coeffs !>@{ !! Diagnostic identifier integer :: id_SN_u=-1, id_SN_v=-1, id_L2u=-1, id_L2v=-1, id_Res_fn = -1 + integer :: id_UH_grad=-1, id_VH_grad=-1, id_L2grad_u=-1, id_L2grad_v=-1 integer :: id_N2_u=-1, id_N2_v=-1, id_S2_u=-1, id_S2_v=-1 integer :: id_dzu=-1, id_dzv=-1, id_dzSxN=-1, id_dzSyN=-1 integer :: id_Rd_dx=-1, id_KH_u_QG = -1, id_KH_v_QG = -1 @@ -611,11 +618,15 @@ end subroutine calc_sqg_struct !> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al. !! style scaling of diffusivity -subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) +subroutine calc_slope_functions(h, uh, vh, tv, dt, G, GV, US, CS, OBC) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)),intent(inout) :: uh !< Layer thickness + !times u [L2 H T-1 ~> m3 s-1 or kg s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),intent(inout) :: vh !< Layer thickness + !times v [L2 H T-1 ~> m3 s-1 or kg s-1] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables real, intent(in) :: dt !< Time increment [T ~> s] type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure @@ -643,6 +654,8 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, halo=1, OBC=OBC) call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS) + elseif (CS%use_gradient_model) then + call calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) else call calc_slope_functions_using_just_e(h, G, GV, US, CS, e) endif @@ -655,8 +668,12 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) if (CS%id_dzSyN > 0) call post_data(CS%id_dzSyN, dzSyN, CS%diag) if (CS%id_SN_u > 0) call post_data(CS%id_SN_u, CS%SN_u, CS%diag) if (CS%id_SN_v > 0) call post_data(CS%id_SN_v, CS%SN_v, CS%diag) + if (CS%id_UH_grad > 0) call post_data(CS%id_UH_grad, CS%UH_grad, CS%diag) + if (CS%id_VH_grad > 0) call post_data(CS%id_VH_grad, CS%VH_grad, CS%diag) if (CS%id_L2u > 0) call post_data(CS%id_L2u, CS%L2u, CS%diag) if (CS%id_L2v > 0) call post_data(CS%id_L2v, CS%L2v, CS%diag) + if (CS%id_L2grad_u > 0) call post_data(CS%id_L2grad_u, CS%L2grad_u, CS%diag) + if (CS%id_L2grad_v > 0) call post_data(CS%id_L2grad_v, CS%L2grad_v, CS%diag) if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag) if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag) endif @@ -1006,7 +1023,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) "%SN_u is not associated with use_variable_mixing.") if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:"// & "%SN_v is not associated with use_variable_mixing.") - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke h_neglect = GV%H_subroundoff @@ -1036,6 +1052,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) ! and midlatitude deformation radii, using calc_resoln_function as a template. !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) + do k=nz,CS%VarMix_Ktop,-1 ! Calculate the interface slopes E_x and E_y and u- and v- points respectively @@ -1062,6 +1079,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) ! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) S2N2_u_local(I,j,k) = (H_geom * S2) * (GV%g_prime(k) / max(Hdn, Hup, CS%h_min_N2) ) enddo ; enddo + + do J=js-1,je ; do i=is,ie S2 = ( E_y(i,J)**2 + 0.25*( & ((E_x(I,j)**2) + (E_x(I-1,j+1)**2)) + ((E_x(I,j+1)**2) + (E_x(I-1,j)**2)) ) ) @@ -1084,10 +1103,11 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) enddo ; enddo ! SN above contains S^2*N^2*H, convert to vertical average of S*N + if (use_dztot) then do I=is-1,ie CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / & - max(dz_tot(i,j), dz_tot(i+1,j), GV%dz_subroundoff) ) + max(dz_tot(i,j), dz_tot(i+1,j), GV%dz_subroundoff) ) enddo else do I=is-1,ie @@ -1106,10 +1126,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie CS%SN_v(i,J) = CS%SN_v(i,J) + S2N2_v_local(i,J,k) enddo ; enddo + + if (use_dztot) then do i=is,ie CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / & - max(dz_tot(i,j), dz_tot(i,j+1), GV%dz_subroundoff) ) + max(dz_tot(i,j), dz_tot(i,j+1), GV%dz_subroundoff) ) enddo else do i=is,ie @@ -1128,6 +1150,96 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) end subroutine calc_slope_functions_using_just_e +!! Computes UH_grad and VH_grad for gradient model (Khani & Dawson, JAMES 2023) +subroutine calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: uh !< Interface height times u + ![L2 H T-1 ~> m3 s-1 or kg s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vh !< Interface height time v + ![L2 H T-1 ~> m3 s-1 or kg s-1] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m] + ! type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + ! Local variables + real :: Ux_Hx(SZIB_(G), SZJ_(G)) ! X-slope of U and H [H L-1 T-1 ~> s-1 or kg m-3 s-1] + real :: Uy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of U and H [H L-1 T-1 ~> s-1 or kg m-3 s-1] + real :: Vx_Hx(SZIB_(G), SZJ_(G)) ! X-slope of V and H [H L-1 T-1 ~> s-1 or kg m-3 s-1] + real :: Vy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of V and H [H L-1 T-1 ~> s-1 or kg m-3 s-1] + real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2] + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected [H ~> m or kg m-2]. + real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2] + real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2]. + ! bathymetric depth for certain calculations. + real :: Lgrid ! Grid lengthscale for the gradient model [H ~> m] + integer :: is, ie, js, je, nz + integer :: i, j, k + integer :: l_seg + + if (.not. CS%initialized) call MOM_error(FATAL, "calc_slope_functions_with_gradient_model: "// & + "Module must be initialized before it is used.") + + if (.not. CS%use_gradient_model) return + if (.not. allocated(CS%UH_grad)) call MOM_error(FATAL, "calc_slope_function:"// & + "UH_grad is not associated with use_gradient_model.") + if (.not. allocated(CS%VH_grad)) call MOM_error(FATAL, "calc_slope_function:"// & + "VH_grad is not associated with use_gradient_model.") + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + h_neglect = GV%H_subroundoff + H_cutoff = real(2*nz) * (GV%Angstrom_H + h_neglect) + + + ! Set the length scale at u-points for the gradient model + do j=js,je ; do I=is-1,ie + CS%L2grad_u(I,j) = (G%dxCu(I,j)**2) + (G%dyCu(I,j)**2) + enddo ; enddo + ! Set length scale at v-points for the gradient model + do J=js-1,je ; do i=is,ie + CS%L2grad_v(i,J) = (G%dxCv(i,J)**2) + (G%dyCv(i,J)**2) + enddo ; enddo + + do k=nz,CS%VarMix_Ktop,-1 + + ! Calculate the gradient slopes Ux_Hx, Vx_Hx, Uy_Hy, Vy_Hy on u- and v-points respectively + do j=js-1,je+1 ; do I=is-1,ie + Ux_Hx(I,j) = ((G%IdxCu(I+1,j)*G%IdyCu(I+1,j)*uh(I+1,j,K)) - (G%IdxCu(I,j)*G%IdyCu(I,j)*uh(I,j,k))) * ( & + G%IareaT(I+1,j) + G%IareaT(I,j)) * G%dyT(I,j) * ((h(i+1,j,k) - h(i,j,k))/( & + h(i+1,j,k) + h(i,j,k) + h_neglect)) + Vx_Hx(I,j) = ((G%IdxCv(I+1,j)*G%IdxCv(I+1,j)*vh(I+1,j,K)) - (G%IdxCv(I,j)*G%IdxCv(I,j)*vh(I,j,k))) * ( & + G%IareaT(I+1,j) + G%IareaT(I,j)) * G%dyT(I,j) * ((h(i+1,j,k) - h(i,j,k))/( & + h(i+1,j,k) + h(i,j,k) + h_neglect)) + ! Mask slopes where interface intersects topography + if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) Ux_Hx(I,j) = 0. + if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) Vx_Hx(I,j) = 0. + enddo ; enddo + do J=js-1,je ; do i=is-1,ie+1 + Uy_Hy(i,J) = ((G%IdyCu(i,J+1)*G%IdyCu(i,J+1)*uh(i,J+1,K)) - (G%IdyCu(i,J)*G%IdyCu(i,J)*uh(i,J,k))) * ( & + G%IareaT(i,J+1) + G%IareaT(i,J)) * G%dxT(i,J) * ((h(i,j+1,k) - h(i,j,k))/( & + h(i,j+1,k) + h(i,j,k) + h_neglect)) + Vy_Hy(i,J) = ((G%IdyCv(i,J+1)*G%IdxCv(i,J+1)*vh(i,J+1,K)) - (G%IdyCv(i,J)*G%IdxCv(i,J)*vh(i,J,k))) * ( & + G%IareaT(i,J+1) + G%IareaT(i,J)) * G%dxT(I,j) * ((h(i,j+1,k) - h(i,j,k))/( & + h(i,j+1,k) + h(i,j,k) + h_neglect)) + ! Mask slopes where interface intersects topography + if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) Uy_Hy(I,j) = 0. + if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) Vy_Hy(I,j) = 0. + enddo ; enddo + + do j=js,je ; do I=is-1,ie + CS%UH_grad(I,j,k) = Ux_Hx(I,j) + 0.25 * ( (Uy_Hy(I,j) + Uy_Hy(I+1,j-1)) + (Uy_Hy(I,j-1) + Uy_Hy(I+1,j)) ) + enddo ; enddo + + do J=js-1,je ; do i=is,ie + CS%VH_grad(i,J,k) = 0.25 * ( (Vx_Hx(i,J) + Vx_Hx(i-1,J+1)) + (Vx_Hx(i-1,J) + Vx_Hx(i,J+1)) ) + Vy_Hy(i,J) + enddo ; enddo + + enddo ! k + +end subroutine calc_slope_functions_with_gradient_model !> Calculates and returns isopycnal slopes with wider halos for use in finding QG viscosity. subroutine calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, CS, OBC) @@ -1310,6 +1422,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) ! for the epipycnal tracer diffusivity [nondim] real :: KhTh_Slope_Cff ! The nondimensional coefficient in the Visbeck formula ! for the interface depth diffusivity [nondim] + real :: Grad_Khani_Scale ! The nondimensional coefficient in the gradient formula + ! for the depth diffusivity [nondim] real :: oneOrTwo ! A variable that may be 1 or 2, depending on which form ! of the equatorial deformation radius us used [nondim] real :: N2_filter_depth ! A depth below which stratification is treated as monotonic when @@ -1357,6 +1471,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) CS%use_simpler_Eady_growth_rate = .false. CS%full_depth_Eady_growth_rate = .false. CS%calculate_depth_fns = .false. + CS%use_gradient_model = .false. ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "USE_VARIABLE_MIXING", CS%use_variable_mixing,& @@ -1365,6 +1480,10 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//& "this is set to true regardless of what is in the "//& "parameter file.", default=.false.) + call get_param(param_file, mdl, "USE_GRADIENT_MODEL", CS%use_gradient_model,& + "If true, use the gradient model formula for eddy diffusivity. This "//& + "allows diagnostics to be created even if the scheme is "//& + "not used.", default=.false.) call get_param(param_file, mdl, "USE_VISBECK", CS%use_Visbeck,& "If true, use the Visbeck et al. (1997) formulation for \n"//& "thickness diffusivity.", default=.false.) @@ -1619,6 +1738,25 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 'm2', conversion=US%L_to_m**2) endif + if (CS%use_gradient_model) then + in_use = .true. + call get_param(param_file, mdl, "Grad_Khani_Scale", CS%grad_Khani_scale, & + "The fixed length scale in the gradient formula.", units="nondim", & + default=1.0) + allocate(CS%UH_grad(IsdB:IedB,jsd:jed,GV%ke), source=0.0) + allocate(CS%VH_grad(isd:ied,JsdB:JedB,GV%ke), source=0.0) + allocate(CS%L2grad_u(IsdB:IedB,jsd:jed), source=0.0) + allocate(CS%L2grad_v(isd:ied,JsdB:JedB), source=0.0) + CS%id_UH_grad = register_diag_field('ocean_model', 'UH_grad', diag%axesCuL, Time, & + 'Inverse gradient eddy time-scale, Ux_Hx+Uy_Hy, at u-points', 's-1') + CS%id_VH_grad = register_diag_field('ocean_model', 'VH_grad', diag%axesCvL, Time, & + 'Inverse gradient eddy time-scale, Vx_Hx+Vy_Hy, at v-points', 's-1') + CS%id_L2grad_u = register_diag_field('ocean_model', 'L2grad_u', diag%axesCu1, Time, & + 'Length scale squared for gradient coefficient, at u-points', 'm2') + CS%id_L2grad_v = register_diag_field('ocean_model', 'L2grad_v', diag%axesCv1, Time, & + 'Length scale squared for gradient coefficient, at v-points', 'm2') + endif + CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, & 'Vertical structure of SQG mode', 'nondim') if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index cf06be4ed2..539f1e2325 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -39,6 +39,7 @@ module MOM_thickness_diffuse logical :: initialized = .false. !< True if this control structure has been initialized. real :: Khth !< Background isopycnal depth diffusivity [L2 T-1 ~> m2 s-1] real :: Khth_Slope_Cff !< Slope dependence coefficient of Khth [nondim] + real :: Grad_Khani_Scale !< Gradient model coefficient [nondim] real :: max_Khth_CFL !< Maximum value of the diffusive CFL for isopycnal height diffusion [nondim] real :: Khth_Min !< Minimum value of Khth [L2 T-1 ~> m2 s-1] real :: Khth_Max !< Maximum value of Khth [L2 T-1 ~> m2 s-1], or 0 for no max @@ -186,7 +187,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp real :: KH_v_lay(SZI_(G),SZJ_(G)) ! Diagnostic of isopycnal height diffusivities at v-points averaged ! to layer centers [L2 T-1 ~> m2 s-1] logical :: use_VarMix, Resoln_scaled, Depth_scaled, use_stored_slopes, khth_use_vert_struct, use_Visbeck - logical :: use_QG_Leith + logical :: use_QG_Leith, use_gradient_model integer :: i, j, k, is, ie, js, je, nz if (.not. CS%initialized) call MOM_error(FATAL, "MOM_thickness_diffuse: "//& @@ -208,13 +209,14 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp Depth_scaled = .false. if (VarMix%use_variable_mixing) then - use_VarMix = VarMix%use_variable_mixing .and. (CS%KHTH_Slope_Cff > 0.) + use_VarMix = VarMix%use_variable_mixing .and. (CS%KHTH_Slope_Cff > 0. .or. CS%grad_Khani_scale > 0.) Resoln_scaled = VarMix%Resoln_scaled_KhTh Depth_scaled = VarMix%Depth_scaled_KhTh use_stored_slopes = VarMix%use_stored_slopes khth_use_vert_struct = allocated(VarMix%khth_struct) use_Visbeck = VarMix%use_Visbeck use_QG_Leith = VarMix%use_QG_Leith_GM + use_gradient_model = VarMix%use_gradient_model if (allocated(VarMix%cg1)) cg1 => VarMix%cg1 else cg1 => null() @@ -333,6 +335,18 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif endif + if (use_VarMix) then + if (use_gradient_model) then !< Gradient model (Khani & Dawson, JAMES 2023) + if (CS%grad_Khani_scale > 0.0) then + !$OMP do + do k=1,nz ; do j=js,je ; do I=is-1,ie + KH_u(I,j,k) = CS%grad_Khani_scale * VarMix%L2grad_u(I,j) * VarMix%UH_grad(I,j,k) + enddo ; enddo ; enddo + endif + endif + endif + + if (CS%use_GME_thickness_diffuse) then !$OMP do do k=1,nz+1 ; do j=js,je ; do I=is-1,ie @@ -437,6 +451,17 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif endif + if (use_VarMix) then + if (use_gradient_model) then !< Gradient model (Khani & Dawson, JAMES 2023) + if (CS%grad_Khani_scale > 0.0) then + !$OMP do + do k=1,nz ; do J=js-1,je ; do i=is,ie + KH_v(i,J,k) = CS%grad_Khani_scale * VarMix%L2grad_v(i,J) * VarMix%VH_grad(i,J,k) + enddo ; enddo ; enddo + endif + endif + endif + if (CS%use_GME_thickness_diffuse) then !$OMP do do k=1,nz+1 ; do J=js-1,je ; do i=is,ie @@ -2237,6 +2262,9 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) call get_param(param_file, mdl, "KHTH_SLOPE_CFF", CS%KHTH_Slope_Cff, & "The nondimensional coefficient in the Visbeck formula for "//& "the interface depth diffusivity", units="nondim", default=0.0) + call get_param(param_file, mdl, "Grad_Khani_Scale", CS%grad_Khani_scale, & + "The nondimensional coefficient in the Gradient model for "//& + "the thickness depth diffusivity", units="nondim", default=1.0) call get_param(param_file, mdl, "KHTH_MIN", CS%KHTH_Min, & "The minimum horizontal thickness diffusivity.", & default=0.0, units="m2 s-1", scale=US%m_to_L**2*US%T_to_s)