From 3a097dc57dbb3c8a863fe16536ce3a52de70fc18 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 28 May 2025 10:55:43 -0400 Subject: [PATCH 01/24] MOM is updated with gradient model for lateral mixing --- src/core/MOM.F90 | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 156a397ff6..89dcdd9fc6 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 [UH ~> m2 s-1 or kg m-1 s-1] + vh => NULL() ! vh : layer thickness times v [VH ~> m2 s-1 or kg m-1 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) From dfafb3fe0bfda241c73a37f6bcb661b3ef23704a Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 28 May 2025 10:56:27 -0400 Subject: [PATCH 02/24] Use gradient model for lateral mixing --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 122 +++++++++++++++++- 1 file changed, 115 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index f2f476b0c8..5dd8980ace 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 (Khani) 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]. @@ -149,6 +154,7 @@ module MOM_lateral_mixing_coeffs logical :: use_Visbeck !< Use Visbeck formulation for thickness diffusivity 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 + real :: grad_Khani_scale !< Fixed length scale in Gradient formula [non-dimension] !! factor [nondim] relating this length scale squared to the cell area real :: Eady_GR_D_scale !< Depth over which to average SN [Z ~> m] real :: Res_coef_khth !< A coefficient [nondim] that determines the function @@ -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,13 @@ 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 [UH ~> m2 s-1 or kg m-1 s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),intent(inout) :: vh !< Layer thickness times v [VH ~> m2 s-1 or kg m-1 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 @@ -644,7 +653,7 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) 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) else - call calc_slope_functions_using_just_e(h, G, GV, US, CS, e) + call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) endif endif @@ -655,8 +664,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 @@ -967,10 +980,13 @@ end subroutine calc_Eady_growth_rate_2D !> The original calc_slope_function() that calculated slopes using !! interface positions only, not accounting for density variations. -subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) +!! Computes UH_grad and VH_grad for gradient (Khani) model (Khani & Dawson, JAMES 2023) +subroutine calc_slope_functions_using_just_e(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 [ZU ~> m2 s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vh !< Interface height times v [ZU ~> m2 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] @@ -980,12 +996,18 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics) real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m] ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m] + real :: U_xH_x(SZIB_(G), SZJ_(G)) ! X-slope of U and H [T-1 ~> s-1] + real :: U_yH_y(SZI_(G), SZJB_(G)) ! Y-slope of U and H [T-1 ~> s-1] + real :: V_xH_x(SZIB_(G), SZJ_(G)) ! X-slope of V and H [T-1 ~> s-1] + real :: V_yH_y(SZI_(G), SZJB_(G)) ! Y-slope of V and H [T-1 ~> s-1] real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2] real :: dZ_cutoff ! A minimum water column depth 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 :: S2 ! Interface slope squared [Z2 L-2 ~> nondim] real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2] + real :: gradUH ! Gradient model frequency, zonal transport [T-1 ~> s-1] + real :: gradVH ! Gradient model frequency, merid transport [T-1 ~> s-1] 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]. real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times @@ -994,6 +1016,9 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) ! the buoyancy frequency squared at v-points [Z T-2 ~> m s-2] logical :: use_dztot ! If true, use the total water column thickness rather than the ! bathymetric depth for certain calculations. + real :: UH_grad_local(SZIB_(G), SZJ_(G),SZK_(GV)) ! The depth integral of grad slopes for UH at u-points + real :: VH_grad_local(SZI_(G), SZJB_(G),SZK_(GV)) ! The depth integral of grad slopes for VH at v-points + real :: Lgrid ! Grid lengthscale for the grad model [H ~> m] integer :: is, ie, js, je, nz integer :: i, j, k integer :: l_seg @@ -1007,6 +1032,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:"// & "%SN_v is not associated with use_variable_mixing.") + 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 @@ -1031,23 +1062,60 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) ! enddo endif + ! To set length scale for gradient (Khani) model (Khani & Dawson, JAMES 2023) ! To set the length scale based on the deformation radius, use wave_speed to ! calculate the first-mode gravity wave speed and then blend the equatorial ! 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) + + ! Set the length scale at u-points. + do j=js,je ; do I=is-1,ie + Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) + CS%L2grad_u(I,j) = 1.0 * Lgrid**2 + enddo ; enddo + ! Set length scale at v-points + do J=js-1,je ; do i=is,ie + Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) + CS%L2grad_v(i,J) = 1.0 * Lgrid**2 + enddo ; enddo + do k=nz,CS%VarMix_Ktop,-1 ! Calculate the interface slopes E_x and E_y and u- and v- points respectively do j=js-1,je+1 ; do I=is-1,ie E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j) ! Mask slopes where interface intersects topography - if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) E_x(I,j) = 0. + if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) E_x(I,j) = 0. enddo ; enddo do J=js-1,je ; do i=is-1,ie+1 E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J) ! Mask slopes where interface intersects topography - if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0. + if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0. + enddo ; enddo + + ! Calculate the gradient slopes U_xH_x, V_xH_x, U_yH_y, V_yH_y on u- and v-points respectively + do j=js-1,je+1 ; do I=is-1,ie + U_xH_x(I,j) =1.0*(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) * (1.0*(h(I+1,j,K) - h(I,j,K))/( & + h(I+1,j,K) + h(I,j,K) + h_neglect)) + V_xH_x(I,j) =1.0*(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) * (1.0*(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) U_xH_x(I,j) = 0. + if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) V_xH_x(I,j) = 0. + enddo ; enddo + do J=js-1,je ; do i=is-1,ie+1 + U_yH_y(i,J) =1.0*(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) * (1.0*(h(i,J+1,K) - h(i,J,K))/( & + h(i,J+1,K) + h(i,J,K) + h_neglect)) + V_yH_y(i,J) =1.0*(G%IdyCv(i,J+1)*G%IdxCv(i,J+1)*vh(i,J,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) * (1.0*(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) U_yH_y(I,j) = 0. + if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) V_yH_y(I,j) = 0. enddo ; enddo ! Calculate N*S*h from this layer and add to the sum @@ -1061,6 +1129,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) H_geom = sqrt(Hdn*Hup) ! 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) ) + gradUH = U_xH_x(I,j) + 0.25*(U_yH_y(I,j)+U_yH_y(I,j-1)+U_yH_y(I+1,j)+U_yH_y(I+1,j-1)) + UH_grad_local(I,j,k) = gradUH enddo ; enddo do J=js-1,je ; do i=is,ie S2 = ( E_y(i,J)**2 + 0.25*( & @@ -1072,6 +1142,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) H_geom = sqrt(Hdn*Hup) ! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) S2N2_v_local(i,J,k) = (H_geom * S2) * (GV%g_prime(k) / (max(Hdn, Hup, CS%h_min_N2))) + gradVH = 0.25*(V_xH_x(i,J)+V_xH_x(i-1,J)+V_xH_x(i,J+1)+V_xH_x(i-1,J+1))+V_yH_y(i,J) + VH_grad_local(i,J,k) = gradVH enddo ; enddo enddo ! k @@ -1081,21 +1153,26 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) do I=is-1,ie ; CS%SN_u(I,j) = 0.0 ; enddo do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie CS%SN_u(I,j) = CS%SN_u(I,j) + S2N2_u_local(I,j,k) + CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) +!! print*, "UH_grad= ", CS%UH_grad(I,j,k) 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) ) +!! CS%UH_grad(I,j) = G%OBCmaskCu(I,j) * ( CS%UH_grad(I,j) / (max(G%bathyT(I,j), G%bathyT(I+1,j)) + G%Z_ref) ) enddo else do I=is-1,ie if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / & (max(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref) ) +!! CS%UH_grad(I,j) = G%OBCmaskCu(I,j) * ( CS%UH_grad(I,j) / (max(G%bathyT(I,j), G%bathyT(I+1,j)) + G%Z_ref) ) else CS%SN_u(I,j) = 0.0 +!! CS%UH_grad(I,j) = 0.0 endif enddo endif @@ -1105,11 +1182,14 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) do i=is,ie ; CS%SN_v(i,J) = 0.0 ; enddo 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) + CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) +!! print*, "VH_grad= ", CS%VH_grad(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) ) +!! CS%VH_grad(i,J) = G%OBCmaskCv(i,J) * (CS%VH_grad(i,J) / (max(G%bathyT(i,J), G%bathyT(i,J+1)) + G%Z_ref) ) enddo else do i=is,ie @@ -1119,8 +1199,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) if ( min(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref > dZ_cutoff ) then CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / & (max(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref) ) +!! CS%VH_grad(i,J) = G%OBCmaskCv(i,J) * (CS%VH_grad(i,J) / (max(G%bathyT(i,J), G%bathyT(i,J+1)) + G%Z_ref) ) else CS%SN_v(i,J) = 0.0 +!! CS%VH_grad(i,J) = 0.0 endif enddo endif @@ -1310,6 +1392,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 +1441,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 +1450,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 +1708,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%axesCu1, Time, & + 'Inverse gradient eddy time-scale, U_xH_x+U_yH_y, at u-points', 's^-1') + CS%id_VH_grad = register_diag_field('ocean_model', 'VH_grad', diag%axesCv1, Time, & + 'Inverse gradient eddy time-scale, V_xH_x+V_yH_y, 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', 'm^2') + CS%id_L2grad_v = register_diag_field('ocean_model', 'L2grad_v', diag%axesCv1, Time, & + 'Length scale squared for gradient coefficient, at v-points', 'm^2') + 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 & From 2903a0a297f6f7fabbaebef770de8c24d7101b8e Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 28 May 2025 10:57:27 -0400 Subject: [PATCH 03/24] Use gradient model for lateral mixing --- .../lateral/MOM_thickness_diffuse.F90 | 32 +++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index cf06be4ed2..a747c99bd9 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) = 1.0*CS%Grad_Khani_Scale*VarMix%L2grad_u(I,j)*VarMix%UH_grad(I,j,k) +!! print*, "KH_u= ", KH_u(I,j,k) + enddo ; enddo ; enddo + 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) = 1.0*CS%Grad_Khani_Scale*VarMix%L2grad_v(i,J)*VarMix%VH_grad(i,J,k) +!! print*, "KH_v=", KH_v(i,J,k) + enddo ; enddo ; enddo + 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) From 176340016aa31a041085bbb3e98d6e61045c0689 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 13:42:03 -0400 Subject: [PATCH 04/24] Use gradient model for lateral mixing --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 115 ++++++++++-------- 1 file changed, 66 insertions(+), 49 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 5dd8980ace..ad0f0cf35f 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -154,7 +154,7 @@ module MOM_lateral_mixing_coeffs logical :: use_Visbeck !< Use Visbeck formulation for thickness diffusivity 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 - real :: grad_Khani_scale !< Fixed length scale in Gradient formula [non-dimension] + real :: GRAD_KHANI_SCALE !< Fixed length scale in Gradient formula [non-dimension] !! factor [nondim] relating this length scale squared to the cell area real :: Eady_GR_D_scale !< Depth over which to average SN [Z ~> m] real :: Res_coef_khth !< A coefficient [nondim] that determines the function @@ -996,10 +996,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics) real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m] ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m] - real :: U_xH_x(SZIB_(G), SZJ_(G)) ! X-slope of U and H [T-1 ~> s-1] - real :: U_yH_y(SZI_(G), SZJB_(G)) ! Y-slope of U and H [T-1 ~> s-1] - real :: V_xH_x(SZIB_(G), SZJ_(G)) ! X-slope of V and H [T-1 ~> s-1] - real :: V_yH_y(SZI_(G), SZJB_(G)) ! Y-slope of V and H [T-1 ~> s-1] + real :: Ux_Hx(SZIB_(G), SZJ_(G)) ! X-slope of U and H [T-1 ~> s-1] + real :: Uy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of U and H [T-1 ~> s-1] + real :: Vx_Hx(SZIB_(G), SZJ_(G)) ! X-slope of V and H [T-1 ~> s-1] + real :: Vy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of V and H [T-1 ~> s-1] real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2] real :: dZ_cutoff ! A minimum water column depth for masking [H ~> m or kg m-2] real :: h_neglect ! A thickness that is so small it is usually lost @@ -1069,16 +1069,18 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) - ! Set the length scale at u-points. - do j=js,je ; do I=is-1,ie - Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) - CS%L2grad_u(I,j) = 1.0 * Lgrid**2 - enddo ; enddo - ! Set length scale at v-points - do J=js-1,je ; do i=is,ie - Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) - CS%L2grad_v(i,J) = 1.0 * Lgrid**2 - enddo ; enddo + ! Set the length scale at u-points for the gradient model + if (CS%use_gradient_model) then + do j=js,je ; do I=is-1,ie + Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) + CS%L2grad_u(I,j) = 1.0 * Lgrid**2 + enddo ; enddo + ! Set length scale at v-points for the gradient model + do J=js-1,je ; do i=is,ie + Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) + CS%L2grad_v(i,J) = 1.0 * Lgrid**2 + enddo ; enddo + endif do k=nz,CS%VarMix_Ktop,-1 @@ -1093,30 +1095,32 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) ! Mask slopes where interface intersects topography if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0. enddo ; enddo - - ! Calculate the gradient slopes U_xH_x, V_xH_x, U_yH_y, V_yH_y on u- and v-points respectively - do j=js-1,je+1 ; do I=is-1,ie - U_xH_x(I,j) =1.0*(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))*( & + + if (CS%use_gradient_model) then + ! 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) =1.0*(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) * (1.0*(h(I+1,j,K) - h(I,j,K))/( & h(I+1,j,K) + h(I,j,K) + h_neglect)) - V_xH_x(I,j) =1.0*(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))*( & + Vx_Hx(I,j) =1.0*(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) * (1.0*(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) U_xH_x(I,j) = 0. - if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) V_xH_x(I,j) = 0. - enddo ; enddo - do J=js-1,je ; do i=is-1,ie+1 - U_yH_y(i,J) =1.0*(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))*( & + ! 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) =1.0*(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) * (1.0*(h(i,J+1,K) - h(i,J,K))/( & h(i,J+1,K) + h(i,J,K) + h_neglect)) - V_yH_y(i,J) =1.0*(G%IdyCv(i,J+1)*G%IdxCv(i,J+1)*vh(i,J,K) - G%IdyCv(i,J)*G%IdxCv(i,J)*vh(i,J,k))*( & + Vy_Hy(i,J) =1.0*(G%IdyCv(i,J+1)*G%IdxCv(i,J+1)*vh(i,J,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) * (1.0*(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) U_yH_y(I,j) = 0. - if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) V_yH_y(I,j) = 0. - enddo ; enddo + ! 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 + endif ! Calculate N*S*h from this layer and add to the sum do j=js,je ; do I=is-1,ie @@ -1129,9 +1133,17 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) H_geom = sqrt(Hdn*Hup) ! 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) ) - gradUH = U_xH_x(I,j) + 0.25*(U_yH_y(I,j)+U_yH_y(I,j-1)+U_yH_y(I+1,j)+U_yH_y(I+1,j-1)) - UH_grad_local(I,j,k) = gradUH + ! gradUH = Ux_Hx(I,j) + 0.25*(Uy_Hy(I,j)+Uy_Hy(I,j-1)+Uy_Hy(I+1,j)+Uy_Hy(I+1,j-1)) + ! UH_grad_local(I,j,k) = gradUH enddo ; enddo + + if (CS%use_gradient_model) then + do j=js,je ; do I=is-1,ie + gradUH = Ux_Hx(I,j) + 0.25*(Uy_Hy(I,j)+Uy_Hy(I,j-1)+Uy_Hy(I+1,j)+Uy_Hy(I+1,j-1)) + UH_grad_local(I,j,k) = gradUH + enddo ; enddo + endif + 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)) ) ) @@ -1142,10 +1154,17 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) H_geom = sqrt(Hdn*Hup) ! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) S2N2_v_local(i,J,k) = (H_geom * S2) * (GV%g_prime(k) / (max(Hdn, Hup, CS%h_min_N2))) - gradVH = 0.25*(V_xH_x(i,J)+V_xH_x(i-1,J)+V_xH_x(i,J+1)+V_xH_x(i-1,J+1))+V_yH_y(i,J) - VH_grad_local(i,J,k) = gradVH + ! gradVH = 0.25*(Vx_Hx(i,J)+Vx_Hx(i-1,J)+Vx_Hx(i,J+1)+Vx_Hx(i-1,J+1))+Vy_Hy(i,J) + ! VH_grad_local(i,J,k) = gradVH enddo ; enddo + if (CS%use_gradient_model) then + do J=js-1,je ; do i=is,ie + gradVH = 0.25*(Vx_Hx(i,J)+Vx_Hx(i-1,J)+Vx_Hx(i,J+1)+Vx_Hx(i-1,J+1))+Vy_Hy(i,J) + VH_grad_local(i,J,k) = gradVH + enddo ; enddo + endif + enddo ! k !$OMP parallel do default(shared) @@ -1153,8 +1172,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) do I=is-1,ie ; CS%SN_u(I,j) = 0.0 ; enddo do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie CS%SN_u(I,j) = CS%SN_u(I,j) + S2N2_u_local(I,j,k) - CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) -!! print*, "UH_grad= ", CS%UH_grad(I,j,k) + if (CS%use_gradient_model) then + CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) + !! print*, "UH_grad= ", CS%UH_grad(I,j,k) + endif enddo ; enddo ! SN above contains S^2*N^2*H, convert to vertical average of S*N @@ -1162,17 +1183,14 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) 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) ) -!! CS%UH_grad(I,j) = G%OBCmaskCu(I,j) * ( CS%UH_grad(I,j) / (max(G%bathyT(I,j), G%bathyT(I+1,j)) + G%Z_ref) ) enddo else do I=is-1,ie if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / & (max(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref) ) -!! CS%UH_grad(I,j) = G%OBCmaskCu(I,j) * ( CS%UH_grad(I,j) / (max(G%bathyT(I,j), G%bathyT(I+1,j)) + G%Z_ref) ) else CS%SN_u(I,j) = 0.0 -!! CS%UH_grad(I,j) = 0.0 endif enddo endif @@ -1182,14 +1200,15 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) do i=is,ie ; CS%SN_v(i,J) = 0.0 ; enddo 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) - CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) -!! print*, "VH_grad= ", CS%VH_grad(I,j,k) + if (CS%use_gradient_model) then + CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) + !! print*, "VH_grad= ", CS%VH_grad(I,j,k) + endif 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) ) -!! CS%VH_grad(i,J) = G%OBCmaskCv(i,J) * (CS%VH_grad(i,J) / (max(G%bathyT(i,J), G%bathyT(i,J+1)) + G%Z_ref) ) enddo else do i=is,ie @@ -1199,10 +1218,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) if ( min(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref > dZ_cutoff ) then CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / & (max(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref) ) -!! CS%VH_grad(i,J) = G%OBCmaskCv(i,J) * (CS%VH_grad(i,J) / (max(G%bathyT(i,J), G%bathyT(i,J+1)) + G%Z_ref) ) else CS%SN_v(i,J) = 0.0 -!! CS%VH_grad(i,J) = 0.0 endif enddo endif @@ -1392,7 +1409,7 @@ 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 + 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] @@ -1710,7 +1727,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) if (CS%use_gradient_model) then in_use = .true. - call get_param(param_file, mdl, "GRAD_Khani_SCALE", CS%grad_Khani_scale, & + 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) @@ -1718,9 +1735,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 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%axesCu1, Time, & - 'Inverse gradient eddy time-scale, U_xH_x+U_yH_y, at u-points', 's^-1') + '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%axesCv1, Time, & - 'Inverse gradient eddy time-scale, V_xH_x+V_yH_y, at v-points', 's^-1') + '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', 'm^2') CS%id_L2grad_v = register_diag_field('ocean_model', 'L2grad_v', diag%axesCv1, Time, & From 9eeee4412d7cde029adba28b74970bf904634414 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 13:42:25 -0400 Subject: [PATCH 05/24] Use gradient model for lateral mixing --- .../lateral/MOM_thickness_diffuse.F90 | 28 ++++++++++--------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index a747c99bd9..474790ecc1 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -39,7 +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 :: 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 @@ -209,7 +209,7 @@ 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.) .or. (CS%Grad_Khani_Scale > 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 @@ -337,12 +337,13 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (use_VarMix) then if (use_gradient_model) then !< Gradient model (Khani & Dawson, JAMES 2023) -!! if (CS%Grad_Khani_Scale > 0.0) then + 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) = 1.0*CS%Grad_Khani_Scale*VarMix%L2grad_u(I,j)*VarMix%UH_grad(I,j,k) -!! print*, "KH_u= ", KH_u(I,j,k) - enddo ; enddo ; enddo + do k=1,nz ; do j=js,je ; do I=is-1,ie + KH_u(I,j,k) = 1.0*CS%GRAD_KHANI_SCALE*VarMix%L2grad_u(I,j)*VarMix%UH_grad(I,j,k) + !! print*, "KH_u= ", KH_u(I,j,k) + enddo ; enddo ; enddo + endif endif endif @@ -453,12 +454,13 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (use_VarMix) then if (use_gradient_model) then !< Gradient model (Khani & Dawson, JAMES 2023) -!! if (CS%Grad_Khani_Scale > 0.0) then + 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) = 1.0*CS%Grad_Khani_Scale*VarMix%L2grad_v(i,J)*VarMix%VH_grad(i,J,k) -!! print*, "KH_v=", KH_v(i,J,k) - enddo ; enddo ; enddo + do k=1,nz ; do J=js-1,je ; do i=is,ie + KH_v(i,J,k) = 1.0*CS%GRAD_KHANI_SCALE*VarMix%L2grad_v(i,J)*VarMix%VH_grad(i,J,k) + !! print*, "KH_v=", KH_v(i,J,k) + enddo ; enddo ; enddo + endif endif endif @@ -2262,7 +2264,7 @@ 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, & + 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, & From 0e766d89bcbed84f35411ce9f55cacb3f1286763 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 13:47:59 -0400 Subject: [PATCH 06/24] Use gradient model for lateral mixing --- src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index ad0f0cf35f..5fdccd66db 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1141,7 +1141,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) do j=js,je ; do I=is-1,ie gradUH = Ux_Hx(I,j) + 0.25*(Uy_Hy(I,j)+Uy_Hy(I,j-1)+Uy_Hy(I+1,j)+Uy_Hy(I+1,j-1)) UH_grad_local(I,j,k) = gradUH - enddo ; enddo + enddo ; enddo endif do J=js-1,je ; do i=is,ie From 50c57181108c331d1da9bfc2f24805e08937fd40 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 13:54:17 -0400 Subject: [PATCH 07/24] Use gradient model for lateral mixing --- src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 5fdccd66db..56b86e3cbd 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1095,7 +1095,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) ! Mask slopes where interface intersects topography if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0. enddo ; enddo - + if (CS%use_gradient_model) then ! 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 From a19eb6b65ca6e4cf5d51c42d0726bdaa20316f03 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 14:09:34 -0400 Subject: [PATCH 08/24] Use gradient model for lateral mixing --- src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 56b86e3cbd..49a895cb9e 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1070,7 +1070,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) ! Set the length scale at u-points for the gradient model - if (CS%use_gradient_model) then + use_gradient_model = CS%use_gradient_model + if (use_gradient_model) then do j=js,je ; do I=is-1,ie Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) CS%L2grad_u(I,j) = 1.0 * Lgrid**2 From 50ca1044f2e34c294550d1e7c80d6a649dd43d7b Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 14:20:36 -0400 Subject: [PATCH 09/24] Use gradient model for lateral mixing --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 49a895cb9e..6777d609ab 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1018,7 +1018,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) ! bathymetric depth for certain calculations. real :: UH_grad_local(SZIB_(G), SZJ_(G),SZK_(GV)) ! The depth integral of grad slopes for UH at u-points real :: VH_grad_local(SZI_(G), SZJB_(G),SZK_(GV)) ! The depth integral of grad slopes for VH at v-points - real :: Lgrid ! Grid lengthscale for the grad model [H ~> m] + real :: Lgrid ! Grid lengthscale for the gradient model [H ~> m] integer :: is, ie, js, je, nz integer :: i, j, k integer :: l_seg @@ -1070,8 +1070,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) ! Set the length scale at u-points for the gradient model - use_gradient_model = CS%use_gradient_model - if (use_gradient_model) then + !if (CS%use_gradient_model) then do j=js,je ; do I=is-1,ie Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) CS%L2grad_u(I,j) = 1.0 * Lgrid**2 @@ -1081,7 +1080,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) CS%L2grad_v(i,J) = 1.0 * Lgrid**2 enddo ; enddo - endif + !endif do k=nz,CS%VarMix_Ktop,-1 From e68af61aa9237023bf6999db45bf59bf673a97b1 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 14:38:05 -0400 Subject: [PATCH 10/24] Use gradient model for lateral mixing --- src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 6777d609ab..f5567493b3 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1070,17 +1070,19 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) ! Set the length scale at u-points for the gradient model - !if (CS%use_gradient_model) then + if (CS%use_gradient_model) then do j=js,je ; do I=is-1,ie Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) CS%L2grad_u(I,j) = 1.0 * Lgrid**2 enddo ; enddo + endif ! Set length scale at v-points for the gradient model + if (CS%use_gradient_model) then do J=js-1,je ; do i=is,ie Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) CS%L2grad_v(i,J) = 1.0 * Lgrid**2 enddo ; enddo - !endif + endif do k=nz,CS%VarMix_Ktop,-1 From 3c7b4b63ea340c1129a3070c896c1b25a02fb0a3 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 6 Jun 2025 14:44:24 -0400 Subject: [PATCH 11/24] Use gradient model for lateral mixing --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index f5567493b3..70d81d3287 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1070,19 +1070,19 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) ! Set the length scale at u-points for the gradient model - if (CS%use_gradient_model) then + !! if (CS%use_gradient_model) then do j=js,je ; do I=is-1,ie Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) CS%L2grad_u(I,j) = 1.0 * Lgrid**2 enddo ; enddo - endif + !! endif ! Set length scale at v-points for the gradient model - if (CS%use_gradient_model) then + !! if (CS%use_gradient_model) then do J=js-1,je ; do i=is,ie Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) CS%L2grad_v(i,J) = 1.0 * Lgrid**2 enddo ; enddo - endif + !! endif do k=nz,CS%VarMix_Ktop,-1 From 9f696ababc6972b575ff00dfa73f33598c3937c0 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 13 Jun 2025 14:44:20 -0400 Subject: [PATCH 12/24] Use gradient model for lateral mixing --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 90 +++++++++---------- 1 file changed, 44 insertions(+), 46 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 70d81d3287..bbd32801b4 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -29,7 +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 (Khani) model (Khani & Dawson, JAMES 2023). + 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. @@ -980,7 +980,7 @@ end subroutine calc_Eady_growth_rate_2D !> The original calc_slope_function() that calculated slopes using !! interface positions only, not accounting for density variations. -!! Computes UH_grad and VH_grad for gradient (Khani) model (Khani & Dawson, JAMES 2023) +!! Computes UH_grad and VH_grad for gradient model (Khani & Dawson, JAMES 2023) subroutine calc_slope_functions_using_just_e(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 @@ -1006,8 +1006,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim] real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2] - real :: gradUH ! Gradient model frequency, zonal transport [T-1 ~> s-1] - real :: gradVH ! Gradient model frequency, merid transport [T-1 ~> s-1] + real :: graduh ! Gradient model frequency, zonal transport [T-1 ~> s-1] + real :: gradvh ! Gradient model frequency, merid transport [T-1 ~> s-1] 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]. real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times @@ -1034,9 +1034,9 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) 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.") + "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.") + "VH_grad is not associated with use_gradient_model.") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -1062,7 +1062,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) ! enddo endif - ! To set length scale for gradient (Khani) model (Khani & Dawson, JAMES 2023) ! To set the length scale based on the deformation radius, use wave_speed to ! calculate the first-mode gravity wave speed and then blend the equatorial ! and midlatitude deformation radii, using calc_resoln_function as a template. @@ -1070,19 +1069,15 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) ! Set the length scale at u-points for the gradient model - !! if (CS%use_gradient_model) then - do j=js,je ; do I=is-1,ie - Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) - CS%L2grad_u(I,j) = 1.0 * Lgrid**2 - enddo ; enddo - !! endif + do j=js,je ; do I=is-1,ie + Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) + CS%L2grad_u(I,j) = 1.0 * Lgrid**2 + enddo ; enddo ! Set length scale at v-points for the gradient model - !! if (CS%use_gradient_model) then - do J=js-1,je ; do i=is,ie - Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) - CS%L2grad_v(i,J) = 1.0 * Lgrid**2 - enddo ; enddo - !! endif + do J=js-1,je ; do i=is,ie + Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) + CS%L2grad_v(i,J) = 1.0 * Lgrid**2 + enddo ; enddo do k=nz,CS%VarMix_Ktop,-1 @@ -1098,13 +1093,13 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0. enddo ; enddo - if (CS%use_gradient_model) then + if (CS%use_gradient_model) then ! 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) =1.0*(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) * (1.0*(h(I+1,j,K) - h(I,j,K))/( & + Ux_Hx(I,j) = 1.0 * (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) * (1.0 * (h(I+1,j,K) - h(I,j,K))/( & h(I+1,j,K) + h(I,j,K) + h_neglect)) - Vx_Hx(I,j) =1.0*(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))*( & + Vx_Hx(I,j) = 1.0 * (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) * (1.0*(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 @@ -1112,10 +1107,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) 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) =1.0*(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))*( & + Uy_Hy(i,J) = 1.0 * (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) * (1.0*(h(i,J+1,K) - h(i,J,K))/( & h(i,J+1,K) + h(i,J,K) + h_neglect)) - Vy_Hy(i,J) =1.0*(G%IdyCv(i,J+1)*G%IdxCv(i,J+1)*vh(i,J,K) - G%IdyCv(i,J)*G%IdxCv(i,J)*vh(i,J,k))*( & + Vy_Hy(i,J) = 1.0 * (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) * (1.0*(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 @@ -1135,14 +1130,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) H_geom = sqrt(Hdn*Hup) ! 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) ) - ! gradUH = Ux_Hx(I,j) + 0.25*(Uy_Hy(I,j)+Uy_Hy(I,j-1)+Uy_Hy(I+1,j)+Uy_Hy(I+1,j-1)) - ! UH_grad_local(I,j,k) = gradUH enddo ; enddo if (CS%use_gradient_model) then do j=js,je ; do I=is-1,ie - gradUH = Ux_Hx(I,j) + 0.25*(Uy_Hy(I,j)+Uy_Hy(I,j-1)+Uy_Hy(I+1,j)+Uy_Hy(I+1,j-1)) - UH_grad_local(I,j,k) = gradUH + graduh = 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)) ) + UH_grad_local(I,j,k) = graduh enddo ; enddo endif @@ -1156,14 +1149,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) H_geom = sqrt(Hdn*Hup) ! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) S2N2_v_local(i,J,k) = (H_geom * S2) * (GV%g_prime(k) / (max(Hdn, Hup, CS%h_min_N2))) - ! gradVH = 0.25*(Vx_Hx(i,J)+Vx_Hx(i-1,J)+Vx_Hx(i,J+1)+Vx_Hx(i-1,J+1))+Vy_Hy(i,J) - ! VH_grad_local(i,J,k) = gradVH enddo ; enddo if (CS%use_gradient_model) then do J=js-1,je ; do i=is,ie - gradVH = 0.25*(Vx_Hx(i,J)+Vx_Hx(i-1,J)+Vx_Hx(i,J+1)+Vx_Hx(i-1,J+1))+Vy_Hy(i,J) - VH_grad_local(i,J,k) = gradVH + gradvh = 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) + VH_grad_local(i,J,k) = gradvh enddo ; enddo endif @@ -1174,13 +1165,16 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) do I=is-1,ie ; CS%SN_u(I,j) = 0.0 ; enddo do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie CS%SN_u(I,j) = CS%SN_u(I,j) + S2N2_u_local(I,j,k) - if (CS%use_gradient_model) then - CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) - !! print*, "UH_grad= ", CS%UH_grad(I,j,k) - endif enddo ; enddo ! SN above contains S^2*N^2*H, convert to vertical average of S*N + if (CS%use_gradient_model) then + do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie + CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) + !! print*, "UH_grad = ", CS%UH_grad(I,j,k) + enddo ; enddo + endif + if (use_dztot) then do I=is-1,ie CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / & @@ -1202,11 +1196,15 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) do i=is,ie ; CS%SN_v(i,J) = 0.0 ; enddo 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) - if (CS%use_gradient_model) then - CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) - !! print*, "VH_grad= ", CS%VH_grad(I,j,k) - endif enddo ; enddo + + if (CS%use_gradient_model) then + do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie + CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) + !! print*, "VH_grad = ", CS%VH_grad(I,j,k) + enddo ; enddo + endif + if (use_dztot) then do i=is,ie CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / & @@ -1729,7 +1727,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) if (CS%use_gradient_model) then in_use = .true. - call get_param(param_file, mdl, "GRAD_KHANI_SCALE", CS%GRAD_KHANI_SCALE, & + 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) @@ -1737,13 +1735,13 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 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%axesCu1, Time, & - 'Inverse gradient eddy time-scale, Ux_Hx+Uy_Hy, at u-points', 's^-1') + '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%axesCv1, Time, & - 'Inverse gradient eddy time-scale, Vx_Hx+Vy_Hy, at v-points', 's^-1') + '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', 'm^2') + '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', 'm^2') + '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, & From 4185bdf119ebd3b7a2765de210474921ac2110c5 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 13 Jun 2025 14:44:31 -0400 Subject: [PATCH 13/24] Use gradient model for lateral mixing --- .../lateral/MOM_thickness_diffuse.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 474790ecc1..08a769dc60 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -209,7 +209,7 @@ 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. .or. CS%GRAD_KHANI_SCALE > 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 @@ -337,11 +337,11 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (use_VarMix) then if (use_gradient_model) then !< Gradient model (Khani & Dawson, JAMES 2023) - if (CS%GRAD_KHANI_SCALE > 0.0) then + 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) = 1.0*CS%GRAD_KHANI_SCALE*VarMix%L2grad_u(I,j)*VarMix%UH_grad(I,j,k) - !! print*, "KH_u= ", KH_u(I,j,k) + KH_u(I,j,k) = 1.0 * CS%grad_Khani_scale * VarMix%L2grad_u(I,j) * VarMix%UH_grad(I,j,k) + !! print*, "KH_u = ", KH_u(I,j,k) enddo ; enddo ; enddo endif endif @@ -454,11 +454,11 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (use_VarMix) then if (use_gradient_model) then !< Gradient model (Khani & Dawson, JAMES 2023) - if (CS%GRAD_KHANI_SCALE > 0.0) then + 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) = 1.0*CS%GRAD_KHANI_SCALE*VarMix%L2grad_v(i,J)*VarMix%VH_grad(i,J,k) - !! print*, "KH_v=", KH_v(i,J,k) + KH_v(i,J,k) = 1.0 * CS%grad_Khani_scale * VarMix%L2grad_v(i,J) * VarMix%VH_grad(i,J,k) + !! print*, "KH_v =", KH_v(i,J,k) enddo ; enddo ; enddo endif endif @@ -2264,7 +2264,7 @@ 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, & + 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, & From c3929dbc56a9da493e0347d96d1b3bcb49b98f7b Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 13 Jun 2025 14:55:04 -0400 Subject: [PATCH 14/24] Use gradient model for lateral mixing --- src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index bbd32801b4..962c9f5344 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1093,7 +1093,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0. enddo ; enddo - if (CS%use_gradient_model) then + if (CS%use_gradient_model) then ! 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) = 1.0 * (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)) * ( & From fa2f67391c0f0f24eaf3057c46bd0585a0d5d88f Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 13 Jun 2025 15:04:22 -0400 Subject: [PATCH 15/24] Use gradient model for lateral mixing --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 962c9f5344..ffefa770ca 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1096,22 +1096,22 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) if (CS%use_gradient_model) then ! 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) = 1.0 * (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) * (1.0 * (h(I+1,j,K) - h(I,j,K))/( & + 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) = 1.0 * (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) * (1.0*(h(I+1,j,K) - h(I,j,K))/( & + 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) = 1.0 * (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) * (1.0*(h(i,J+1,K) - h(i,J,K))/( & + 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) = 1.0 * (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) * (1.0*(h(i,J+1,K) - h(i,J,K))/( & + 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. From b5e70bd98e5e51b687e0bc5639466a3c28ec6f78 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 9 Jul 2025 15:44:41 -0400 Subject: [PATCH 16/24] Use gradient model for lateral mixing --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index ffefa770ca..f1e9fd40d1 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -154,7 +154,7 @@ module MOM_lateral_mixing_coeffs logical :: use_Visbeck !< Use Visbeck formulation for thickness diffusivity 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 - real :: GRAD_KHANI_SCALE !< Fixed length scale in Gradient formula [non-dimension] + real :: Grad_Khani_Scale !< Fixed length scale in Gradient formula [non-dimension] !! factor [nondim] relating this length scale squared to the cell area real :: Eady_GR_D_scale !< Depth over which to average SN [Z ~> m] real :: Res_coef_khth !< A coefficient [nondim] that determines the function @@ -1097,25 +1097,25 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) ! 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)) + 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)) + 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. + 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)) + 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)) + 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. + 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 endif @@ -1409,7 +1409,7 @@ 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 + 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] @@ -1727,7 +1727,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) if (CS%use_gradient_model) then in_use = .true. - call get_param(param_file, mdl, "GRAD_KHANI_SCALE", CS%grad_Khani_scale, & + 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) From ddefb6052beb8b2a3adf3e3cfe50944627bc6c5a Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 9 Jul 2025 15:44:53 -0400 Subject: [PATCH 17/24] Use gradient model for lateral mixing --- src/parameterizations/lateral/MOM_thickness_diffuse.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 08a769dc60..40e18b0671 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -39,7 +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 :: 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 @@ -2264,7 +2264,7 @@ 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, & + 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, & From 720a0c04c736306015f5f38b4b894215734c1347 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Mon, 14 Jul 2025 11:27:14 -0400 Subject: [PATCH 18/24] Center indices in h(i,j,k) --- src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index f1e9fd40d1..b5927511dd 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1085,12 +1085,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) do j=js-1,je+1 ; do I=is-1,ie E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j) ! Mask slopes where interface intersects topography - if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) E_x(I,j) = 0. + if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) E_x(I,j) = 0. enddo ; enddo do J=js-1,je ; do i=is-1,ie+1 E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J) ! Mask slopes where interface intersects topography - if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0. + if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(I,j) = 0. enddo ; enddo if (CS%use_gradient_model) then From 59e59b033c701cacf676219e754228d7d69a640c Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 30 Jul 2025 11:13:53 -0400 Subject: [PATCH 19/24] A subroutine for the gradient model caculations --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 211 +++++++++++------- 1 file changed, 127 insertions(+), 84 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index b5927511dd..bd951370ae 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -154,8 +154,8 @@ module MOM_lateral_mixing_coeffs logical :: use_Visbeck !< Use Visbeck formulation for thickness diffusivity 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 - real :: Grad_Khani_Scale !< Fixed length scale in Gradient formula [non-dimension] !! 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: @@ -652,11 +652,17 @@ subroutine calc_slope_functions(h, uh, vh, 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, uh, vh) + call calc_slope_functions_using_just_e(h, G, GV, US, CS, e) endif endif +!! if (CS%use_gradient_model) then +!! call calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) +!! endif + if (query_averaging_enabled(CS%diag)) then if (CS%id_dzu > 0) call post_data(CS%id_dzu, dzu, CS%diag) if (CS%id_dzv > 0) call post_data(CS%id_dzv, dzv, CS%diag) @@ -980,13 +986,10 @@ end subroutine calc_Eady_growth_rate_2D !> The original calc_slope_function() that calculated slopes using !! interface positions only, not accounting for density variations. -!! Computes UH_grad and VH_grad for gradient model (Khani & Dawson, JAMES 2023) -subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) +subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) 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 [ZU ~> m2 s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vh !< Interface height times v [ZU ~> m2 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] @@ -996,18 +999,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics) real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m] ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m] - real :: Ux_Hx(SZIB_(G), SZJ_(G)) ! X-slope of U and H [T-1 ~> s-1] - real :: Uy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of U and H [T-1 ~> s-1] - real :: Vx_Hx(SZIB_(G), SZJ_(G)) ! X-slope of V and H [T-1 ~> s-1] - real :: Vy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of V and H [T-1 ~> s-1] real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2] real :: dZ_cutoff ! A minimum water column depth 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 :: S2 ! Interface slope squared [Z2 L-2 ~> nondim] real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2] - real :: graduh ! Gradient model frequency, zonal transport [T-1 ~> s-1] - real :: gradvh ! Gradient model frequency, merid transport [T-1 ~> s-1] 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]. real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times @@ -1016,9 +1013,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) ! the buoyancy frequency squared at v-points [Z T-2 ~> m s-2] logical :: use_dztot ! If true, use the total water column thickness rather than the ! bathymetric depth for certain calculations. - real :: UH_grad_local(SZIB_(G), SZJ_(G),SZK_(GV)) ! The depth integral of grad slopes for UH at u-points - real :: VH_grad_local(SZI_(G), SZJB_(G),SZK_(GV)) ! The depth integral of grad slopes for VH at v-points - real :: Lgrid ! Grid lengthscale for the gradient model [H ~> m] integer :: is, ie, js, je, nz integer :: i, j, k integer :: l_seg @@ -1031,13 +1025,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) "%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.") - - 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 @@ -1068,16 +1055,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2) - ! Set the length scale at u-points for the gradient model - do j=js,je ; do I=is-1,ie - Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) - CS%L2grad_u(I,j) = 1.0 * Lgrid**2 - enddo ; enddo - ! Set length scale at v-points for the gradient model - do J=js-1,je ; do i=is,ie - Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) - CS%L2grad_v(i,J) = 1.0 * Lgrid**2 - enddo ; enddo do k=nz,CS%VarMix_Ktop,-1 @@ -1090,35 +1067,9 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) do J=js-1,je ; do i=is-1,ie+1 E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J) ! Mask slopes where interface intersects topography - if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(I,j) = 0. + if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0. enddo ; enddo - if (CS%use_gradient_model) then - ! 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 - endif - ! Calculate N*S*h from this layer and add to the sum do j=js,je ; do I=is-1,ie S2 = ( E_x(I,j)**2 + 0.25*( & @@ -1132,12 +1083,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) S2N2_u_local(I,j,k) = (H_geom * S2) * (GV%g_prime(k) / max(Hdn, Hup, CS%h_min_N2) ) enddo ; enddo - if (CS%use_gradient_model) then - do j=js,je ; do I=is-1,ie - graduh = 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)) ) - UH_grad_local(I,j,k) = graduh - enddo ; enddo - endif do J=js-1,je ; do i=is,ie S2 = ( E_y(i,J)**2 + 0.25*( & @@ -1151,13 +1096,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) S2N2_v_local(i,J,k) = (H_geom * S2) * (GV%g_prime(k) / (max(Hdn, Hup, CS%h_min_N2))) enddo ; enddo - if (CS%use_gradient_model) then - do J=js-1,je ; do i=is,ie - gradvh = 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) - VH_grad_local(i,J,k) = gradvh - enddo ; enddo - endif - enddo ! k !$OMP parallel do default(shared) @@ -1168,12 +1106,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) enddo ; enddo ! SN above contains S^2*N^2*H, convert to vertical average of S*N - if (CS%use_gradient_model) then - do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie - CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) - !! print*, "UH_grad = ", CS%UH_grad(I,j,k) - enddo ; enddo - endif if (use_dztot) then do I=is-1,ie @@ -1198,12 +1130,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) CS%SN_v(i,J) = CS%SN_v(i,J) + S2N2_v_local(i,J,k) enddo ; enddo - if (CS%use_gradient_model) then - do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie - CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) - !! print*, "VH_grad = ", CS%VH_grad(I,j,k) - enddo ; enddo - endif if (use_dztot) then do i=is,ie @@ -1227,6 +1153,123 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, uh, vh) 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 [ZU ~> m2 s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vh !< Interface height times v [ZU ~> m2 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 :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m] + ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m] + real :: Ux_Hx(SZIB_(G), SZJ_(G)) ! X-slope of U and H [T-1 ~> s-1] + real :: Uy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of U and H [T-1 ~> s-1] + real :: Vx_Hx(SZIB_(G), SZJ_(G)) ! X-slope of V and H [T-1 ~> s-1] + real :: Vy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of V and H [T-1 ~> s-1] + real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2] + !real :: dZ_cutoff ! A minimum water column depth 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 :: graduh ! Gradient model frequency, zonal transport [T-1 ~> s-1] + real :: gradvh ! Gradient model frequency, merid transport [T-1 ~> s-1] + 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]. + !logical :: use_dztot ! If true, use the total water column thickness rather than the + ! bathymetric depth for certain calculations. + real :: UH_grad_local(SZIB_(G), SZJ_(G),SZK_(GV)) ! The depth integral of grad slopes for UH at u-points + real :: VH_grad_local(SZI_(G), SZJB_(G),SZK_(GV)) ! The depth integral of grad slopes for VH at v-points + 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 + Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) + CS%L2grad_u(I,j) = 1.0 * Lgrid**2 + enddo ; enddo + ! Set length scale at v-points for the gradient model + do J=js-1,je ; do i=is,ie + Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) + CS%L2grad_v(i,J) = 1.0 * Lgrid**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 + graduh = 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)) ) + UH_grad_local(I,j,k) = graduh + enddo ; enddo + + do J=js-1,je ; do i=is,ie + gradvh = 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) + VH_grad_local(i,J,k) = gradvh + enddo ; enddo + + enddo ! k + + !$OMP parallel do default(shared) + do j=js,je + + do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie + CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) + enddo ; enddo + + enddo + !$OMP parallel do default(shared) + do J=js-1,je + + do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie + CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) + enddo ; enddo + + enddo + +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) From 4f5378e3233737f4a0e1f0112212f04520dfb629 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 27 Aug 2025 14:25:38 -0400 Subject: [PATCH 20/24] Cleaning up code blocks related to the gradient model --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 58 +++++-------------- 1 file changed, 16 insertions(+), 42 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index bd951370ae..d129be2c3c 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -623,8 +623,8 @@ subroutine calc_slope_functions(h, uh, vh, tv, dt, G, GV, US, CS, OBC) 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 [UH ~> m2 s-1 or kg m-1 s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),intent(inout) :: vh !< Layer thickness times v [VH ~> m2 s-1 or kg m-1 s-1] + 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 @@ -659,10 +659,6 @@ subroutine calc_slope_functions(h, uh, vh, tv, dt, G, GV, US, CS, OBC) endif endif -!! if (CS%use_gradient_model) then -!! call calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) -!! endif - if (query_averaging_enabled(CS%diag)) then if (CS%id_dzu > 0) call post_data(CS%id_dzu, dzu, CS%diag) if (CS%id_dzv > 0) call post_data(CS%id_dzv, dzv, CS%diag) @@ -1055,7 +1051,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e) !$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 @@ -1158,8 +1153,8 @@ 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 [ZU ~> m2 s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vh !< Interface height times v [ZU ~> m2 s-1] + 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 times v [L2 G 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] @@ -1167,22 +1162,22 @@ subroutine calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) ! Local variables !real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m] ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m] - real :: Ux_Hx(SZIB_(G), SZJ_(G)) ! X-slope of U and H [T-1 ~> s-1] - real :: Uy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of U and H [T-1 ~> s-1] - real :: Vx_Hx(SZIB_(G), SZJ_(G)) ! X-slope of V and H [T-1 ~> s-1] - real :: Vy_Hy(SZI_(G), SZJB_(G)) ! Y-slope of V and H [T-1 ~> s-1] + 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 :: dZ_cutoff ! A minimum water column depth 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 :: graduh ! Gradient model frequency, zonal transport [T-1 ~> s-1] - real :: gradvh ! Gradient model frequency, merid transport [T-1 ~> s-1] +! real :: graduh ! Gradient model frequency, zonal transport [T-1 ~> s-1] +! real :: gradvh ! Gradient model frequency, merid transport [T-1 ~> s-1] 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]. !logical :: use_dztot ! If true, use the total water column thickness rather than the ! bathymetric depth for certain calculations. - real :: UH_grad_local(SZIB_(G), SZJ_(G),SZK_(GV)) ! The depth integral of grad slopes for UH at u-points - real :: VH_grad_local(SZI_(G), SZJB_(G),SZK_(GV)) ! The depth integral of grad slopes for VH at v-points +! real :: UH_grad_local(SZIB_(G), SZJ_(G),SZK_(GV)) ! The depth integral of grad slopes for UH at u-points +! real :: VH_grad_local(SZI_(G), SZJB_(G),SZK_(GV)) ! The depth integral of grad slopes for VH at v-points real :: Lgrid ! Grid lengthscale for the gradient model [H ~> m] integer :: is, ie, js, je, nz integer :: i, j, k @@ -1205,13 +1200,11 @@ subroutine calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) ! Set the length scale at u-points for the gradient model do j=js,je ; do I=is-1,ie - Lgrid = sqrt(G%dxCu(I,j)**2 + G%dyCu(I,j)**2) - CS%L2grad_u(I,j) = 1.0 * Lgrid**2 + 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 - Lgrid = sqrt(G%dxCv(i,J)**2 + G%dyCv(i,J)**2) - CS%L2grad_v(i,J) = 1.0 * Lgrid**2 + CS%L2grad_v(i,J) = (G%dxCv(i,J)**2) + (G%dyCv(i,J)**2) enddo ; enddo do k=nz,CS%VarMix_Ktop,-1 @@ -1241,34 +1234,15 @@ subroutine calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) enddo ; enddo do j=js,je ; do I=is-1,ie - graduh = 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)) ) - UH_grad_local(I,j,k) = graduh + 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 - gradvh = 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) - VH_grad_local(i,J,k) = gradvh + 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 - !$OMP parallel do default(shared) - do j=js,je - - do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie - CS%UH_grad(I,j,k) = UH_grad_local(I,j,k) - enddo ; enddo - - enddo - !$OMP parallel do default(shared) - do J=js-1,je - - do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie - CS%VH_grad(i,J,k) = VH_grad_local(i,J,k) - enddo ; enddo - - enddo - end subroutine calc_slope_functions_with_gradient_model !> Calculates and returns isopycnal slopes with wider halos for use in finding QG viscosity. From 0ece3d789d01233c576f14e6b71cceee4dd5e785 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 27 Aug 2025 14:25:50 -0400 Subject: [PATCH 21/24] Cleaning up code blocks related to the gradient model --- src/parameterizations/lateral/MOM_thickness_diffuse.F90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 40e18b0671..539f1e2325 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -340,8 +340,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp 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) = 1.0 * CS%grad_Khani_scale * VarMix%L2grad_u(I,j) * VarMix%UH_grad(I,j,k) - !! print*, "KH_u = ", KH_u(I,j,k) + 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 @@ -457,8 +456,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp 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) = 1.0 * CS%grad_Khani_scale * VarMix%L2grad_v(i,J) * VarMix%VH_grad(i,J,k) - !! print*, "KH_v =", KH_v(i,J,k) + 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 From bd2815af108dab9a080cca3b8c75f93e339de182 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 27 Aug 2025 14:29:29 -0400 Subject: [PATCH 22/24] Correcting the units of h, uh and vh --- src/core/MOM.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 89dcdd9fc6..a434d50020 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1191,8 +1191,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, & 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] - uh => NULL(), & ! uh : layer thickness times u [UH ~> m2 s-1 or kg m-1 s-1] - vh => NULL() ! vh : layer thickness times v [VH ~> m2 s-1 or kg m-1 s-1] + 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. From b36afef46553d0854869b31278daa111edef9042 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Wed, 27 Aug 2025 15:25:37 -0400 Subject: [PATCH 23/24] Cleaning up code blocks related to the gradient model --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index d129be2c3c..07048f8d02 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -623,8 +623,10 @@ subroutine calc_slope_functions(h, uh, vh, tv, dt, G, GV, US, CS, OBC) 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] + 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 @@ -1153,8 +1155,10 @@ 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 times v [L2 G T-1 ~> m3 s-1 or kg s-1] + 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] From d1daf4705ed49dcd5770e1d1d8b872e52360a651 Mon Sep 17 00:00:00 2001 From: Sina Khani Date: Fri, 12 Dec 2025 14:16:49 -0500 Subject: [PATCH 24/24] Correction on 3-d diag registration plus cleaning up the gradient model subroutine --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 07048f8d02..e42cb4c5d1 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1164,24 +1164,16 @@ subroutine calc_slope_functions_with_gradient_model(h, G, GV, US, CS, e, uh, vh) 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 :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m] - ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m] 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 :: dZ_cutoff ! A minimum water column depth 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 :: graduh ! Gradient model frequency, zonal transport [T-1 ~> s-1] -! real :: gradvh ! Gradient model frequency, merid transport [T-1 ~> s-1] 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]. - !logical :: use_dztot ! If true, use the total water column thickness rather than the ! bathymetric depth for certain calculations. -! real :: UH_grad_local(SZIB_(G), SZJ_(G),SZK_(GV)) ! The depth integral of grad slopes for UH at u-points -! real :: VH_grad_local(SZI_(G), SZJB_(G),SZK_(GV)) ! The depth integral of grad slopes for VH at v-points real :: Lgrid ! Grid lengthscale for the gradient model [H ~> m] integer :: is, ie, js, je, nz integer :: i, j, k @@ -1755,9 +1747,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 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%axesCu1, Time, & + 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%axesCv1, Time, & + 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')