From 90c438f2911e379ba4aa90589a6a4ab0c9b94d1c Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Mon, 3 Mar 2025 15:16:18 -0500 Subject: [PATCH 1/8] Updates MOM_kappa_shear averaging and diagnostics. - Update to horizontal averaging procedure in MOM_kappa_shear for VERTEX_SHEAR = True. - Adds geometric and thickness weighted averages for moving diffusivity from vertices to tracer points. - Adds new diagnostics of N2, S2, Kd, and TKE both with and without VERTEX_SHEAR, useful for debugging. --- .../vertical/MOM_kappa_shear.F90 | 319 ++++++++++++++++-- 1 file changed, 287 insertions(+), 32 deletions(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index b3f8fb4da1..8a92b3a525 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -97,8 +97,12 @@ module MOM_kappa_shear !! time average TKE when there is mass in all layers. Otherwise always !! report the time-averaged TKE, as is currently done when there !! are some massless layers. + integer :: Kd_tracer_mean_answer_date !< Answer date for Kd tracer point horizontal averaging logical :: VS_viscosity_bug !< If true, use a bug in the calculation of the viscosity that sets !! it to zero for all vertices that are on a coastline. + logical :: VS_GeometricMean !< If true use geometric averaging for Kd from vertices to tracer points + logical :: VS_ThicknessMean !< If true use thickness weighting when averaging Kd from vertices to + !!tracer points logical :: restrictive_tolerance_check !< If false, uses the less restrictive tolerance check to !! determine if a timestep is acceptable for the KS_it outer iteration !! loop, as the code was originally written. True uses the more @@ -109,7 +113,8 @@ module MOM_kappa_shear type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. !>@{ Diagnostic IDs - integer :: id_Kd_shear = -1, id_TKE = -1 + integer :: id_Kd_shear = -1, id_TKE = -1, id_Kd_vertex = -1, & + id_S2_init = -1, id_N2_init = -1, id_S2_mean = -1, id_N2_mean = -1 !>@} end type Kappa_shear_CS @@ -151,6 +156,11 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & !! call to kappa_shear_init. ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & + diag_N2_init, & ! Diagnostic of N2 as provided to this routine [T-2 ~> s-2] + diag_S2_init, & ! Diagnostic of S2 as provided to this routine [T-2 ~> s-2] + diag_N2_mean, & ! Diagnostic of N2 averaged over the timestep applied [T-2 ~> s-2] + diag_S2_mean ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2] real, dimension(SZI_(G),SZK_(GV)) :: & h_2d, & ! A 2-D version of h [H ~> m or kg m-2]. dz_2d, & ! Vertical distance between interface heights [Z ~> m]. @@ -172,7 +182,12 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s] tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2]. kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s] - tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. + tke_avg, & ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. + N2_init, & ! N2 as provided to this routine [T-2 ~> s-2]. + S2_init, & ! S2 as provided to this routine [T-2 ~> s-2]. + N2_mean, & ! The time-weighted average of N2 [T-2 ~> s-2]. + S2_mean ! The time-weighted average of S2 [T-2 ~> s-2]. + real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2]. real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa]. @@ -196,6 +211,11 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & k0dt = dt*CS%kappa_0 dz_massless = 0.1*sqrt((US%Z_to_m*GV%m_to_H)*k0dt) + if (CS%id_N2_init>0) diag_N2_init(:,:,:) = 0.0 + if (CS%id_S2_init>0) diag_S2_init(:,:,:) = 0.0 + if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0 + if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 + !$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,tv,G,GV,US, & !$OMP CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io) do j=js,je @@ -284,7 +304,6 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & nzc = nz do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo endif - f2 = 0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + & (G%Coriolis2Bu(I,J-1) + G%Coriolis2Bu(I-1,J))) surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = p_surf(i,j) @@ -297,7 +316,8 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & h_lay, dz_lay, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, & - tke_avg, tv, CS, GV, US) + tke_avg, N2_init, S2_init, N2_mean, S2_mean, & + tv, CS, GV, US) ! call cpu_clock_begin(id_clock_setup) ! Extrapolate from the vertically reduced grid back to the original layers. @@ -310,6 +330,18 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & tke_2d(i,K) = tke_avg(K) endif enddo + if (CS%id_N2_mean>0) then ; do K=1,nz+1 + diag_N2_mean(i,j,K) = N2_mean(K) + enddo ; endif + if (CS%id_S2_mean>0) then ; do K=1,nz+1 + diag_S2_mean(i,j,K) = S2_mean(K) + enddo ; endif + if (CS%id_N2_init>0) then ; do K=1,nz+1 + diag_N2_init(i,j,K) = N2_init(K) + enddo ; endif + if (CS%id_S2_init>0) then ; do K=1,nz+1 + diag_S2_init(i,j,K) = S2_init(K) + enddo ; endif else do K=1,nz+1 if (kf(K) == 0.0) then @@ -322,6 +354,31 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & kf(K) * tke_avg(kc(K)+1) endif enddo + do K=1,nz+1 + if (kf(K) == 0.0) then + if (CS%id_N2_mean>0) diag_N2_mean(i,j,K) = N2_mean(kc(K)) + if (CS%id_S2_mean>0) diag_S2_mean(i,j,K) = S2_mean(kc(K)) + if (CS%id_N2_init>0) diag_N2_init(i,j,K) = N2_init(kc(K)) + if (CS%id_S2_init>0) diag_S2_init(i,j,K) = S2_init(kc(K)) + else + if (CS%id_N2_mean>0) then + diag_N2_mean(i,j,K) = (1.0-kf(K)) * N2_mean(kc(K)) & + + kf(K) *N2_mean(kc(K)+1) + endif + if (CS%id_S2_mean>0) then + diag_S2_mean(i,j,K) = (1.0-kf(K)) * S2_mean(kc(K)) & + + kf(K) *S2_mean(kc(K)+1) + endif + if (CS%id_N2_init>0) then + diag_N2_init(i,j,K) = (1.0-kf(K)) * N2_init(kc(K)) & + + kf(K) *N2_init(kc(K)+1) + endif + if (CS%id_S2_init>0) then + diag_S2_init(i,j,K) = (1.0-kf(K)) * S2_init(kc(K)) & + + kf(K) *S2_init(kc(K)+1) + endif + endif + enddo endif ! call cpu_clock_end(id_clock_setup) else ! Land points, still inside the i-loop. @@ -334,6 +391,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & kappa_io(i,j,K) = G%mask2dT(i,j) * kappa_2d(i,K) tke_io(i,j,K) = G%mask2dT(i,j) * tke_2d(i,K) kv_io(i,j,K) = ( G%mask2dT(i,j) * kappa_2d(i,K) ) * CS%Prandtl_turb + enddo ; enddo enddo ! end of j-loop @@ -345,6 +403,10 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag) if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag) + if (CS%id_N2_init > 0) call post_data(CS%id_N2_init, diag_N2_init, CS%diag) + if (CS%id_S2_init > 0) call post_data(CS%id_S2_init, diag_S2_init, CS%diag) + if (CS%id_N2_mean > 0) call post_data(CS%id_N2_mean, diag_N2_mean, CS%diag) + if (CS%id_S2_mean > 0) call post_data(CS%id_S2_mean, diag_S2_mean, CS%diag) end subroutine Calculate_kappa_shear @@ -387,15 +449,23 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ !! call to kappa_shear_init. ! Local variables + real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: & + diag_N2_init, & ! Diagnostic of N2 as provided to this routine [T-2 ~> s-2] + diag_S2_init, & ! Diagnostic of S2 as provided to this routine [T-2 ~> s-2] + diag_N2_mean, & ! Diagnostic of N2 averaged over the timestep applied [T-2 ~> s-2] + diag_S2_mean, & ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2] + diag_Kd_vertex ! Diagnostic of Kd at vertex [H Z T-1 ~> m2 s-1 or Pa s] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & dz_3d ! Vertical distance between interface heights [Z ~> m]. real, dimension(SZIB_(G),SZK_(GV)) :: & - h_2d, & ! A 2-D version of h [H ~> m or kg m-2]. dz_2d, & ! Vertical distance between interface heights [Z ~> m]. u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1]. T_2d, S_2d, rho_2d ! 2-D versions of T [C ~> degC], S [S ~> ppt], and rho [R ~> kg m-3]. real, dimension(SZIB_(G),SZK_(GV)+1,2) :: & kappa_2d ! Quasi 2-D versions of kappa_io [H Z T-1 ~> m2 s-1 or Pa s] + real, dimension(SZIB_(G),SZK_(GV),2) :: & + h_2d ! Quasi 2-D version of h [H ~> m or kg m-2]. real, dimension(SZIB_(G),SZK_(GV)+1) :: & tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2]. real, dimension(SZK_(GV)) :: & @@ -410,7 +480,12 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s] tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2]. kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s] - tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. + tke_avg, & ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. + N2_init, & ! N2 as provided to this routine [T-2 ~> s-2]. + S2_init, & ! S2 as provided to this routine [T-2 ~> s-2]. + N2_mean, & ! The time-weighted average of N2 [T-2 ~> s-2]. + S2_mean ! The time-weighted average of S2 [T-2 ~> s-2]. + real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2]. real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa]. @@ -427,11 +502,18 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! merged into nearby massive layers. real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for ! interpolating back to the original index space [nondim]. + real :: a1, a2, a3, a4, b1, b2, b3, b4, bT ! Terms in horizontal averaging of diffusivity to tracer location integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc, J2, J2m1 ! Diagnostics that should be deleted? isB = G%isc-1 ; ieB = G%iecB ; jsB = G%jsc-1 ; jeB = G%jecB ; nz = GV%ke + if (CS%id_N2_init>0) diag_N2_init(:,:,:) = 0.0 + if (CS%id_S2_init>0) diag_S2_init(:,:,:) = 0.0 + if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0 + if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 + if (CS%id_Kd_vertex>0) diag_Kd_vertex(:,:,:) = 0.0 + use_temperature = associated(tv%T) k0dt = dt*CS%kappa_0 @@ -469,7 +551,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ (G%mask2dT(i+1,j) * (h(i+1,j,k) * S_in(i+1,j,k)) + & G%mask2dT(i,j+1) * (h(i,j+1,k) * S_in(i,j+1,k))) ) * I_hwt endif - h_2d(I,k) = ((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + & + h_2d(I,k,J2) = ((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + & (G%mask2dT(i+1,j) * h(i+1,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) ) / & ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + & (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 ) @@ -503,23 +585,23 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Add a new layer if this one has mass. ! if ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless)) nzc = nzc+1 if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. & - (h_2d(I,k) > dz_massless)) nzc = nzc+1 + (h_2d(I,k,J2) > dz_massless)) nzc = nzc+1 ! Only merge clusters of massless layers. ! if ((h_lay(nzc) > dz_massless) .or. & ! ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless))) nzc = nzc+1 kc(k) = nzc - h_lay(nzc) = h_lay(nzc) + h_2d(I,k) + h_lay(nzc) = h_lay(nzc) + h_2d(I,k,J2) dz_lay(nzc) = dz_lay(nzc) + dz_2d(I,k) - u0xdz(nzc) = u0xdz(nzc) + u_2d(I,k)*h_2d(I,k) - v0xdz(nzc) = v0xdz(nzc) + v_2d(I,k)*h_2d(I,k) + u0xdz(nzc) = u0xdz(nzc) + u_2d(I,k)*h_2d(I,k,J2) + v0xdz(nzc) = v0xdz(nzc) + v_2d(I,k)*h_2d(I,k,J2) if (use_temperature) then - T0xdz(nzc) = T0xdz(nzc) + T_2d(I,k)*h_2d(I,k) - S0xdz(nzc) = S0xdz(nzc) + S_2d(I,k)*h_2d(I,k) + T0xdz(nzc) = T0xdz(nzc) + T_2d(I,k)*h_2d(I,k,J2) + S0xdz(nzc) = S0xdz(nzc) + S_2d(I,k)*h_2d(I,k,J2) else - T0xdz(nzc) = T0xdz(nzc) + rho_2d(I,k)*h_2d(I,k) - S0xdz(nzc) = S0xdz(nzc) + rho_2d(I,k)*h_2d(I,k) + T0xdz(nzc) = T0xdz(nzc) + rho_2d(I,k)*h_2d(I,k,J2) + S0xdz(nzc) = S0xdz(nzc) + rho_2d(I,k)*h_2d(I,k,J2) endif enddo kc(nz+1) = nzc+1 @@ -529,18 +611,18 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Now determine kf, the fractional weight of interface kc when ! interpolating between interfaces kc and kc+1. - kf(1) = 0.0 ; dz_in_lay = h_2d(I,1) + kf(1) = 0.0 ; dz_in_lay = h_2d(I,1,J2) do k=2,nz if (kc(k) > kc(k-1)) then - kf(k) = 0.0 ; dz_in_lay = h_2d(I,k) + kf(k) = 0.0 ; dz_in_lay = h_2d(I,k,J2) else - kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(I,k) + kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(I,k,J2) endif enddo kf(nz+1) = 0.0 else do k=1,nz - h_lay(k) = h_2d(I,k) + h_lay(k) = h_2d(I,k,J2) dz_lay(k) = dz_2d(I,k) u0xdz(k) = u_2d(I,k)*h_lay(k) ; v0xdz(k) = v_2d(I,k)*h_lay(k) enddo @@ -579,7 +661,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & h_lay, dz_lay, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, & - tke_avg, tv, CS, GV, US) + tke_avg, N2_init, S2_init, N2_mean, S2_mean, tv, CS, GV, US) ! call cpu_clock_begin(Id_clock_setup) ! Extrapolate from the vertically reduced grid back to the original layers. if (nz == nzc) then @@ -591,6 +673,21 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ tke_2d(i,K) = tke_avg(K) endif enddo + if (CS%id_N2_mean>0) then ; do K=1,nz+1 + diag_N2_mean(i,j,K) = N2_mean(K) + enddo ; endif + if (CS%id_S2_mean>0) then ; do K=1,nz+1 + diag_S2_mean(i,j,K) = S2_mean(K) + enddo ; endif + if (CS%id_N2_init>0) then ; do K=1,nz+1 + diag_N2_init(i,j,K) = N2_init(K) + enddo ; endif + if (CS%id_S2_init>0) then ; do K=1,nz+1 + diag_S2_init(i,j,K) = S2_init(K) + enddo ; endif + if (CS%id_Kd_vertex>0) then ; do K=1,nz+1 + diag_Kd_vertex(i,j,K) = kappa_avg(K) + enddo ; endif else do K=1,nz+1 if (kf(K) == 0.0) then @@ -601,6 +698,36 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ tke_2d(I,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1) endif enddo + do K=1,nz+1 + if (kf(K) == 0.0) then + if (CS%id_N2_mean>0) diag_N2_mean(I,J,K) = N2_mean(kc(K)) + if (CS%id_S2_mean>0) diag_S2_mean(I,J,K) = S2_mean(kc(K)) + if (CS%id_N2_init>0) diag_N2_init(I,J,K) = N2_init(kc(K)) + if (CS%id_S2_init>0) diag_S2_init(I,J,K) = S2_init(kc(K)) + if (CS%id_Kd_vertex>0) diag_Kd_vertex(I,J,K) = kappa_avg(kc(K)) + else + if (CS%id_N2_mean>0) then + diag_N2_mean(I,J,K) = (1.0-kf(K)) * N2_mean(kc(K)) & + + kf(K) *N2_mean(kc(K)+1) + endif + if (CS%id_S2_mean>0) then + diag_S2_mean(I,J,K) = (1.0-kf(K)) * S2_mean(kc(K)) & + + kf(K) *S2_mean(kc(K)+1) + endif + if (CS%id_N2_init>0) then + diag_N2_init(I,J,K) = (1.0-kf(K)) * N2_init(kc(K)) & + + kf(K) *N2_init(kc(K)+1) + endif + if (CS%id_S2_init>0) then + diag_S2_init(I,J,K) = (1.0-kf(K)) * S2_init(kc(K)) & + + kf(K) *S2_init(kc(K)+1) + endif + if (CS%id_Kd_vertex>0) then + diag_Kd_vertex(I,J,K) = (1.0-kf(K)) * kappa_avg(kc(K)) & + + kf(K) *kappa_avg(kc(K)+1) + endif + endif + enddo endif ! call cpu_clock_end(Id_clock_setup) else ! Land points, still inside the i-loop. @@ -620,13 +747,75 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ kv_io(I,J,K) = kappa_2d(I,K,J2) * CS%Prandtl_turb enddo; enddo endif - if (J>=G%jsc) then ; do K=1,nz+1 ; do i=G%isc,G%iec - ! Set the diffusivities in tracer columns from the values at vertices. - kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * & - ((kappa_2d(I-1,K,J2m1) + kappa_2d(I,K,J2)) + & - (kappa_2d(I-1,K,J2) + kappa_2d(I,K,J2m1))) - enddo ; enddo ; endif - + if (J>=G%jsc) then + if (CS%Kd_tracer_mean_answer_date < 20250304) then + ! This is a non-thickness weighted arithmetic mean, but it is not bitwise identical to the more general + ! version of the mean below the "else" that includes a mathematically equivalent form. + do K=1,nz+1 ; do i=G%isc,G%iec + ! Set the diffusivities in tracer columns from the values at vertices. + kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * & + ((kappa_2d(I-1,K,J2m1)+kappa_2d(I,K,J2)) +& + (kappa_2d(I-1,K,J2)+kappa_2d(I,K,J2m1))) + enddo; enddo + else + do i=G%isc,G%iec + ! It doesn't matter what these values are so make them zero. + kappa_io(i,1,j) = 0.0 + kappa_io(i,nz+1,j) = 0.0 + enddo + if (CS%VS_GeometricMean) then + do K=2,nz ; do i=G%isc,G%iec + ! Set the diffusivities in tracer columns from the values at vertices. + ! The geometric mean is zero if any component is zero, hence the need + ! and sensitivity to a value of kappa_trunc in regions on boundaries of shear zones. + a1 = max(kappa_2d(I-1,K,J2m1),CS%kappa_trunc) ! Diffusivity, lower left + a2 = max(kappa_2d(I,K,J2),CS%kappa_trunc) ! Diffusivity, upper right + a3 = max(kappa_2d(I-1,K,J2),CS%kappa_trunc) ! Diffusivity, upper left + a4 = max(kappa_2d(I,K,J2m1),CS%kappa_trunc) ! Difusivity, lower right + if (CS%VS_ThicknessMean) then + b1 = h_2d(I-1,k,J2m1)+h_2d(I-1,k-1,J2m1) ! 2x thickness, lower left + b2 = h_2d(I,k,J2)+h_2d(I,k-1,J2) ! 2x thickness, upper right + b3 = h_2d(I-1,k,J2)+h_2d(I-1,k-1,J2) ! 2x thickness, upper left + b4 = h_2d(I,k,J2m1)+h_2d(I,k-1,J2m1) ! 2x thickness, upper right + else + b1 = 1. + b2 = 1. + b3 = 1. + b4 = 1. + endif + bT = (b1 + b2) + (b3 + b4) + ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: + if (bT>GV%H_subroundoff) then + kappa_io(i,j,K) = G%mask2dT(i,j) * & + ( (a1**(b1/bT)*a2**(b2/bT)) * (a3**(b3/bT)*a4**(b4/bT)) ) + endif + enddo; enddo + else !arithmetic mean + do K=2,nz ; do i=G%isc,G%iec + a1 = kappa_2d(I-1,K,J2m1) ! Diffusivity, lower left + a2 = kappa_2d(I,K,J2) ! Diffusivity, upper right + a3 = kappa_2d(I-1,K,J2) ! Diffusivity, upper left + a4 = kappa_2d(I,K,J2m1) ! Difusivity, lower right + if (CS%VS_ThicknessMean) then + b1 = h_2d(I-1,k,J2m1)+h_2d(I-1,k-1,J2m1) ! 2x thickness, lower left + b2 = h_2d(I,k,J2)+h_2d(I,k-1,J2) ! 2x thickness, upper right + b3 = h_2d(I-1,k,J2)+h_2d(I-1,k-1,J2) ! 2x thickness, upper left + b4 = h_2d(I,k,J2m1)+h_2d(I,k-1,J2m1) ! 2x thickness, upper right + else + b1 = 1. + b2 = 1. + b3 = 1. + b4 = 1. + endif + bT = (b1+b2) + (b3+b4) + ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: + if (bT>GV%H_subroundoff) then + kappa_io(i,j,K) = G%mask2dT(i,j)*((a1*b1+a2*b2)+(a3*b3+a4*b4))/bT + endif + enddo; enddo + endif + endif + endif enddo ! end of J-loop if (CS%debug) then @@ -636,13 +825,19 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag) if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag) + if (CS%id_Kd_vertex > 0) call post_data(CS%id_Kd_vertex, diag_Kd_vertex, CS%diag) + if (CS%id_N2_init > 0) call post_data(CS%id_N2_init, diag_N2_init, CS%diag) + if (CS%id_S2_init > 0) call post_data(CS%id_S2_init, diag_S2_init, CS%diag) + if (CS%id_N2_mean > 0) call post_data(CS%id_N2_mean, diag_N2_mean, CS%diag) + if (CS%id_S2_mean > 0) call post_data(CS%id_S2_mean, diag_S2_mean, CS%diag) end subroutine Calc_kappa_shear_vertex !> This subroutine calculates shear-driven diffusivity and TKE in a single column subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, hlay, dz_lay, & - u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, tke_avg, tv, CS, GV, US) + u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, tke_avg, N2_init, S2_init, & + N2_mean, S2_mean, tv, CS, GV, US ) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZK_(GV)+1), & intent(inout) :: kappa !< The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s] @@ -669,6 +864,14 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, hlay, dz_la intent(out) :: kappa_avg !< The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s] real, dimension(SZK_(GV)+1), & intent(out) :: tke_avg !< The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. + real, dimension(SZK_(GV)+1), & + intent(out) :: N2_mean !< The time-weighted average of N2 [Z2 T-2 ~> m2 s-2]. + real, dimension(SZK_(GV)+1), & + intent(out) :: S2_mean !< The time-weighted average of S2 [Z2 T-2 ~> m2 s-2]. + real, dimension(SZK_(GV)+1), & + intent(out) :: N2_init !< The initial value of N2 [Z2 T-2 ~> m2 s-2]. + real, dimension(SZK_(GV)+1), & + intent(out) :: S2_init !< The initial value of S2 [Z2 T-2 ~> m2 s-2]. real, intent(in) :: dt !< Time increment [T ~> s]. type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any !! available thermodynamic fields. Absent fields @@ -935,13 +1138,19 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, hlay, dz_la ! This call just calculates N2 and S2. call calculate_projected_state(kappa, u, v, T, Sal, 0.0, nzc, hlay, I_dz_int, dbuoy_dT, dbuoy_dS, & CS%vel_underflow, u, v, T, Sal, N2, S2, GV, US) + do K=1,nzc+1 + N2_init(K) = N2(K) + S2_init(K) = S2(K) + enddo + + ! ---------------------------------------------------- ! Iterate ! ---------------------------------------------------- dt_rem = dt do K=1,nzc+1 K_Q(K) = 0.0 - kappa_avg(K) = 0.0 ; tke_avg(K) = 0.0 + kappa_avg(K) = 0.0 ; tke_avg(K) = 0.0 ; N2_mean(K) = 0.0 ; S2_mean(K) = 0.0 local_src_avg(K) = 0.0 ! Use the grid spacings to scale errors in the source. if ( h_Int(K) > 0.0 ) & @@ -1113,8 +1322,11 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, hlay, dz_la kappa_mid(K) = 0.5*(kappa_out(K) + kappa_pred(K)) kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt tke_avg(K) = tke_avg(K) + dt_wt*0.5*(tke_pred(K) + tke(K)) + N2_mean(K) = N2_mean(K) + dt_wt*N2(K) + S2_mean(K) = S2_mean(K) + dt_wt*S2(K) kappa(K) = kappa_pred(K) ! First guess for the next iteration. enddo + ! call cpu_clock_end(id_clock_avg) endif @@ -1886,6 +2098,22 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) "If true, use a bug in vertex shear that zeros out viscosities at "//& "vertices on coastlines.", & default=.true., do_not_log=just_read.or.(.not.CS%KS_at_vertex)) + call get_param(param_file, mdl, "VERTEX_SHEAR_GEOMETRIC_MEAN", CS%VS_GeometricMean, & + "If true, use a geometric mean for moving diffusivity from "//& + "vertices to tracer points. False uses algebraic mean.", & + default=.false., do_not_log=just_read.or.(.not.CS%KS_at_vertex)) + call get_param(param_file, mdl, "VERTEX_SHEAR_THICKNESS_MEAN", CS%VS_ThicknessMean, & + "If true, apply thickness weighting to horizontal averagings of diffusivity "//& + "to tracer points in the kappa shear solver.", & + default=.false.) + call get_param(param_file, mdl, "VERTEX_SHEAR_KD_MEAN_ANSWER_DATE", CS%Kd_tracer_mean_answer_date, & + "Answer date for the algorithm for arithmetic horizontal mean of "//& + "Kd from vertices to tracer points. Less than 20250304 uses an older algorithm.", & + default=20250303, do_not_log=just_read.or.(.not.CS%KS_at_vertex)) + if ((CS%Kd_tracer_mean_answer_date<20250304).and.(CS%VS_ThicknessMean .or. CS%VS_GeometricMean)) then + call MOM_error(FATAL,"Cannot use VERTEX_SHEAR_THICKNESS_MEAN or VERTEX_SHEAR_GEOMETRIC_MEAN "//& + "with VERTEX_SHEAR_KD_MEAN_ANSWER_DATE < 20250304") + endif call get_param(param_file, mdl, "RINO_CRIT", CS%RiNo_crit, & "The critical Richardson number for shear mixing.", & units="nondim", default=0.25, do_not_log=just_read) @@ -2026,9 +2254,36 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) CS%diag => diag CS%id_Kd_shear = register_diag_field('ocean_model','Kd_shear', diag%axesTi, Time, & - 'Shear-driven Diapycnal Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) - CS%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesTi, Time, & - 'Shear-driven Turbulent Kinetic Energy', 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) + 'Shear-driven Diapycnal Diffusivity at horizontal tracer points', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) + if (CS%KS_at_vertex) then + CS%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesBi, Time, & + 'Shear-driven Turbulent Kinetic Energy at horizontal vertices', 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) + CS%id_Kd_vertex = register_diag_field('ocean_model','Kd_shear_vertex', diag%axesBi, Time, & + 'Shear-driven Diapycnal Diffusivity at horizontal vertices', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_S2_init = register_diag_field('ocean_model','S2_shear_in', diag%axesBi, Time, & + 'Interface shear squared at horizontal vertices, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) + CS%id_N2_init = register_diag_field('ocean_model','N2_shear_in', diag%axesBi, Time, & + 'Interface straitification at horizontal vertices, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) + CS%id_S2_mean = register_diag_field('ocean_model','S2_shear_mean', diag%axesBi, Time, & + 'Interface shear squared at horizontal vertices, averaged over timestep in kappa-shear', & + 's-2', conversion=US%s_to_T**2) + CS%id_N2_mean = register_diag_field('ocean_model','N2_shear_mean', diag%axesBi, Time, & + 'Interface straitification at horizontal vertices, averaged over timestep in kappa-shear', & + 's-2', conversion=US%s_to_T**2) + else + CS%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesTi, Time, & + 'Shear-driven Turbulent Kinetic Energy at horizontal tracer points', 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) + CS%id_S2_init = register_diag_field('ocean_model','S2_shear_in', diag%axesTi, Time, & + 'Interface shear squared at horizontal tracer points, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) + CS%id_N2_init = register_diag_field('ocean_model','N2_shear_in', diag%axesTi, Time, & + 'Interface straitification at horizontal tracer points, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) + CS%id_S2_mean = register_diag_field('ocean_model','S2_shear_mean', diag%axesTi, Time, & + 'Interface shear squared at horizontal tracer points, averaged over timestep in kappa-shear', & + 's-2', conversion=US%s_to_T**2) + CS%id_N2_mean = register_diag_field('ocean_model','N2_shear_mean', diag%axesTi, Time, & + 'Interface straitification at horizontal tracer points, averaged ove timestep in kappa-shear', & + 's-2', conversion=US%s_to_T**2) + endif end function kappa_shear_init From c79c5f67eb26edecb8a7e4571d735d7532d9e635 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Mon, 3 Mar 2025 15:36:18 -0500 Subject: [PATCH 2/8] Shorten two lines in MOM_kappa_shear to comply w/ MOM6 standard --- src/parameterizations/vertical/MOM_kappa_shear.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 8a92b3a525..f8b3dcd4ed 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -2272,11 +2272,13 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) 's-2', conversion=US%s_to_T**2) else CS%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesTi, Time, & - 'Shear-driven Turbulent Kinetic Energy at horizontal tracer points', 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) + 'Shear-driven Turbulent Kinetic Energy at horizontal tracer points', 'm2 s-2', & + conversion=US%Z_to_m**2*US%s_to_T**2) CS%id_S2_init = register_diag_field('ocean_model','S2_shear_in', diag%axesTi, Time, & 'Interface shear squared at horizontal tracer points, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) CS%id_N2_init = register_diag_field('ocean_model','N2_shear_in', diag%axesTi, Time, & - 'Interface straitification at horizontal tracer points, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) + 'Interface straitification at horizontal tracer points, as input to kappa-shear', 's-2', & + conversion=US%s_to_T**2) CS%id_S2_mean = register_diag_field('ocean_model','S2_shear_mean', diag%axesTi, Time, & 'Interface shear squared at horizontal tracer points, averaged over timestep in kappa-shear', & 's-2', conversion=US%s_to_T**2) From 867044fcc86582d8e01180d108bcf54e9756fc87 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Wed, 5 Mar 2025 08:57:58 -0500 Subject: [PATCH 3/8] More updates for Kappa_Shear horizontal averaging - Fix OMP directives for new diagnostics - Remove quasi-2d fields that broke OMP threading - Reworking horizontal averaging algorithm. --- .../vertical/MOM_kappa_shear.F90 | 253 ++++++++++-------- 1 file changed, 145 insertions(+), 108 deletions(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index f8b3dcd4ed..457e5d4590 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -217,7 +217,8 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 !$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,tv,G,GV,US, & - !$OMP CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io) + !$OMP CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io, & + !$OMP diag_N2_init,diag_S2_init,diag_N2_mean,diag_S2_mean) do j=js,je ! Convert layer thicknesses into geometric thickness in height units. @@ -453,19 +454,23 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ diag_N2_init, & ! Diagnostic of N2 as provided to this routine [T-2 ~> s-2] diag_S2_init, & ! Diagnostic of S2 as provided to this routine [T-2 ~> s-2] diag_N2_mean, & ! Diagnostic of N2 averaged over the timestep applied [T-2 ~> s-2] - diag_S2_mean, & ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2] - diag_Kd_vertex ! Diagnostic of Kd at vertex [H Z T-1 ~> m2 s-1 or Pa s] + diag_S2_mean ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & dz_3d ! Vertical distance between interface heights [Z ~> m]. real, dimension(SZIB_(G),SZK_(GV)) :: & + h_2d, & ! A 2-D version of h [H ~> m or kg m-2]. dz_2d, & ! Vertical distance between interface heights [Z ~> m]. u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1]. T_2d, S_2d, rho_2d ! 2-D versions of T [C ~> degC], S [S ~> ppt], and rho [R ~> kg m-3]. - real, dimension(SZIB_(G),SZK_(GV)+1,2) :: & - kappa_2d ! Quasi 2-D versions of kappa_io [H Z T-1 ~> m2 s-1 or Pa s] - real, dimension(SZIB_(G),SZK_(GV),2) :: & - h_2d ! Quasi 2-D version of h [H ~> m or kg m-2]. + real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: & + kappa_vertex, & ! kappa_io on vertices [H Z T-1 ~> m2 s-1 or Pa s] + h_int_mask ! masked h sumed over bounding layers on vertices [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & + hweight_ul, & ! thickness weighting at cell center from upper-left vertex + hweight_ur, & ! thickness weighting at cell center from upper-right vertex + hweight_ll, & ! thickness weighting at cell center from lower-left vertex + hweight_lr ! thickness weighting at cell center from lower-right vertex real, dimension(SZIB_(G),SZK_(GV)+1) :: & tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2]. real, dimension(SZK_(GV)) :: & @@ -494,6 +499,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ real :: dz_massless ! A layer thickness that is considered massless [H ~> m or kg m-2] real :: I_hwt ! The inverse of the masked thickness weights [H-1 ~> m-1 or m2 kg-1]. real :: I_Prandtl ! The inverse of the turbulent Prandtl number [nondim]. + real :: hweight_denom ! A weighting factor for thickness averaging [H ~> m or kg m-2] logical :: use_temperature ! If true, temperature and salinity have been ! allocated and are being used as state variables. @@ -502,8 +508,10 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! merged into nearby massive layers. real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for ! interpolating back to the original index space [nondim]. - real :: a1, a2, a3, a4, b1, b2, b3, b4, bT ! Terms in horizontal averaging of diffusivity to tracer location - integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc, J2, J2m1 + ! Terms in horizontal averaging of diffusivity to tracer location + real :: Kd_ll, Kd_lr, Kd_ul, Kd_ur, weight_ll, weight_lr, weight_ul, weight_ur + ! Other useful indices + integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc ! Diagnostics that should be deleted? isB = G%isc-1 ; ieB = G%iecB ; jsB = G%jsc-1 ; jeB = G%jecB ; nz = GV%ke @@ -512,7 +520,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_S2_init>0) diag_S2_init(:,:,:) = 0.0 if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0 if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 - if (CS%id_Kd_vertex>0) diag_Kd_vertex(:,:,:) = 0.0 + kappa_vertex(:,:,:) = 0.0 use_temperature = associated(tv%T) @@ -523,10 +531,58 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Convert layer thicknesses into geometric thickness in height units. call thickness_to_dz(h, tv, dz_3d, G, GV, US, halo_size=1) + if (CS%VS_ThicknessMean) then + !Loop over halo and pre compute a masked thickness around each interface + ! (excluding the top and bottom boundary) + do j=G%jsc-1,G%jec+1; do K=2,nz; do i=G%isc-1,G%iec+1 + h_int_mask(i,j,K) = G%mask2dT(i,j)*(h(i,j,k-1)+h(i,j,k))/2. + enddo; enddo ; enddo + !Compute the thickness weighting if requested for horizontal mean + do J=JsB,JeB ; do K=2,nz ; do I=IsB,IeB + hweight_denom = h_int_mask(i-1,j-1,k) + & !left,bot + 2.*h_int_mask(i-1,j,k) + & !left,cen + h_int_mask(i-1,j+1,k) + & !left,top + 2. * h_int_mask(i,j-1,k) + &!mid,bot + 4. * h_int_mask(i,j,k) + & !mid,cen + 2. * h_int_mask(i,j+1,k) + &!mid,top + h_int_mask(i+1,j-1,k) + & !right,bot + 2. * h_int_mask(i-1,j,k) + &!right,cen + h_int_mask(i+1,j+1,k) + & !right,top + GV%H_subroundoff + ! Computing the weighting of the 4 vertices to each cell center + hweight_ul(I,J,K) = ((h_int_mask(i-1,j,k) + & !left,cen + h_int_mask(i-1,j+1,k)) + &!left,top + (h_int_mask(i,j,k) + & !mid,cen + h_int_mask(i,j+1,k)) & !mid,top + ) / hweight_denom + hweight_ur(I,J,k) = ((h_int_mask(i,j,k) + & !mid,cen + h_int_mask(i,j+1,k)) + &!mid,top + (h_int_mask(i+1,j,k) + & !right,cen + h_int_mask(i+1,j+1,k)) &!right,top + ) / hweight_denom + hweight_ll(I,J,k) = ((h_int_mask(i-1,j-1,k) + &!left,bot + h_int_mask(i-1,j,k)) + & !left,cen + (h_int_mask(i+1,j-1,k) + &!mid,bot + h_int_mask(i+1,j,k)) & !mid,cen + ) / hweight_denom + hweight_lr(I,J,k) = ((h_int_mask(i,j-1,k) + & !mid,bot + h_int_mask(i,j,k)) + & !mid,cen + (h_int_mask(i+1,j-1,k) + &!right,bot + h_int_mask(i+1,j,k)) & !right,cen + ) / hweight_denom + enddo ; enddo ; enddo + else + ! Set some weights for horizontal averaging factors + weight_ll=0.25 + weight_lr=0.25 + weight_ul=0.25 + weight_ur=0.25 + endif + !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,tv,G,GV, & - !$OMP US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl) + !$OMP US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl, & + !$OMP diag_N2_init,diag_S2_init,diag_N2_mean,diag_S2_mean) do J=JsB,JeB - J2 = mod(J,2)+1 ; J2m1 = 3-J2 ! = mod(J-1,2)+1 ! Interpolate the various quantities to the corners, using masks. do k=1,nz ; do I=IsB,IeB @@ -551,7 +607,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ (G%mask2dT(i+1,j) * (h(i+1,j,k) * S_in(i+1,j,k)) + & G%mask2dT(i,j+1) * (h(i,j+1,k) * S_in(i,j+1,k))) ) * I_hwt endif - h_2d(I,k,J2) = ((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + & + h_2d(I,k) = ((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + & (G%mask2dT(i+1,j) * h(i+1,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) ) / & ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + & (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 ) @@ -585,23 +641,23 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Add a new layer if this one has mass. ! if ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless)) nzc = nzc+1 if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. & - (h_2d(I,k,J2) > dz_massless)) nzc = nzc+1 + (h_2d(I,k) > dz_massless)) nzc = nzc+1 ! Only merge clusters of massless layers. ! if ((h_lay(nzc) > dz_massless) .or. & ! ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless))) nzc = nzc+1 kc(k) = nzc - h_lay(nzc) = h_lay(nzc) + h_2d(I,k,J2) + h_lay(nzc) = h_lay(nzc) + h_2d(I,k) dz_lay(nzc) = dz_lay(nzc) + dz_2d(I,k) - u0xdz(nzc) = u0xdz(nzc) + u_2d(I,k)*h_2d(I,k,J2) - v0xdz(nzc) = v0xdz(nzc) + v_2d(I,k)*h_2d(I,k,J2) + u0xdz(nzc) = u0xdz(nzc) + u_2d(I,k)*h_2d(I,k) + v0xdz(nzc) = v0xdz(nzc) + v_2d(I,k)*h_2d(I,k) if (use_temperature) then - T0xdz(nzc) = T0xdz(nzc) + T_2d(I,k)*h_2d(I,k,J2) - S0xdz(nzc) = S0xdz(nzc) + S_2d(I,k)*h_2d(I,k,J2) + T0xdz(nzc) = T0xdz(nzc) + T_2d(I,k)*h_2d(I,k) + S0xdz(nzc) = S0xdz(nzc) + S_2d(I,k)*h_2d(I,k) else - T0xdz(nzc) = T0xdz(nzc) + rho_2d(I,k)*h_2d(I,k,J2) - S0xdz(nzc) = S0xdz(nzc) + rho_2d(I,k)*h_2d(I,k,J2) + T0xdz(nzc) = T0xdz(nzc) + rho_2d(I,k)*h_2d(I,k) + S0xdz(nzc) = S0xdz(nzc) + rho_2d(I,k)*h_2d(I,k) endif enddo kc(nz+1) = nzc+1 @@ -611,18 +667,18 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Now determine kf, the fractional weight of interface kc when ! interpolating between interfaces kc and kc+1. - kf(1) = 0.0 ; dz_in_lay = h_2d(I,1,J2) + kf(1) = 0.0 ; dz_in_lay = h_2d(I,1) do k=2,nz if (kc(k) > kc(k-1)) then - kf(k) = 0.0 ; dz_in_lay = h_2d(I,k,J2) + kf(k) = 0.0 ; dz_in_lay = h_2d(I,k) else - kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(I,k,J2) + kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(I,k) endif enddo kf(nz+1) = 0.0 else do k=1,nz - h_lay(k) = h_2d(I,k,J2) + h_lay(k) = h_2d(I,k) dz_lay(k) = dz_2d(I,k) u0xdz(k) = u_2d(I,k)*h_lay(k) ; v0xdz(k) = v_2d(I,k)*h_lay(k) enddo @@ -666,7 +722,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Extrapolate from the vertically reduced grid back to the original layers. if (nz == nzc) then do K=1,nz+1 - kappa_2d(I,K,J2) = kappa_avg(K) + kappa_vertex(I,J,K) = kappa_avg(K) if (CS%all_layer_TKE_bug) then tke_2d(i,K) = tke(K) else @@ -685,16 +741,13 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_S2_init>0) then ; do K=1,nz+1 diag_S2_init(i,j,K) = S2_init(K) enddo ; endif - if (CS%id_Kd_vertex>0) then ; do K=1,nz+1 - diag_Kd_vertex(i,j,K) = kappa_avg(K) - enddo ; endif else do K=1,nz+1 if (kf(K) == 0.0) then - kappa_2d(I,K,J2) = kappa_avg(kc(K)) + kappa_vertex(I,J,K) = kappa_avg(kc(K)) tke_2d(I,K) = tke_avg(kc(K)) else - kappa_2d(I,K,J2) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1) + kappa_vertex(I,J,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1) tke_2d(I,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1) endif enddo @@ -704,7 +757,6 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_S2_mean>0) diag_S2_mean(I,J,K) = S2_mean(kc(K)) if (CS%id_N2_init>0) diag_N2_init(I,J,K) = N2_init(kc(K)) if (CS%id_S2_init>0) diag_S2_init(I,J,K) = S2_init(kc(K)) - if (CS%id_Kd_vertex>0) diag_Kd_vertex(I,J,K) = kappa_avg(kc(K)) else if (CS%id_N2_mean>0) then diag_N2_mean(I,J,K) = (1.0-kf(K)) * N2_mean(kc(K)) & @@ -722,101 +774,86 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ diag_S2_init(I,J,K) = (1.0-kf(K)) * S2_init(kc(K)) & + kf(K) *S2_init(kc(K)+1) endif - if (CS%id_Kd_vertex>0) then - diag_Kd_vertex(I,J,K) = (1.0-kf(K)) * kappa_avg(kc(K)) & - + kf(K) *kappa_avg(kc(K)+1) - endif endif enddo endif ! call cpu_clock_end(Id_clock_setup) else ! Land points, still inside the i-loop. do K=1,nz+1 - kappa_2d(I,K,J2) = 0.0 ; tke_2d(I,K) = 0.0 + kappa_vertex(I,J,K) = 0.0 ; tke_2d(I,K) = 0.0 enddo endif ; enddo ! i-loop if (CS%VS_viscosity_bug) then do K=1,nz+1 ; do I=IsB,IeB tke_io(I,J,K) = G%mask2dBu(I,J) * tke_2d(I,K) - kv_io(I,J,K) = ( G%mask2dBu(I,J) * kappa_2d(I,K,J2) ) * CS%Prandtl_turb + kv_io(I,J,K) = ( G%mask2dBu(I,J) * kappa_vertex(I,J,K) ) * CS%Prandtl_turb enddo; enddo else do K=1,nz+1 ; do I=IsB,IeB tke_io(I,J,K) = tke_2d(I,K) - kv_io(I,J,K) = kappa_2d(I,K,J2) * CS%Prandtl_turb + kv_io(I,J,K) = kappa_vertex(I,J,K) * CS%Prandtl_turb enddo; enddo endif - if (J>=G%jsc) then - if (CS%Kd_tracer_mean_answer_date < 20250304) then - ! This is a non-thickness weighted arithmetic mean, but it is not bitwise identical to the more general - ! version of the mean below the "else" that includes a mathematically equivalent form. - do K=1,nz+1 ; do i=G%isc,G%iec + enddo ! end of J-loop + + + do j=G%jsc,G%jec + if (CS%Kd_tracer_mean_answer_date < 20250304) then + ! This is a non-thickness weighted arithmetic mean, but it is not bitwise identical to the more general + ! version of the mean below the "else" that includes a mathematically equivalent form. + do K=1,nz+1 ; do i=G%isc,G%iec + ! Set the diffusivities in tracer columns from the values at vertices. + kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * & + ((kappa_vertex(I-1,J-1,K)+kappa_vertex(I,J,K)) +& + (kappa_vertex(I-1,J,K)+kappa_vertex(I,J-1,K))) + enddo; enddo + else + do i=G%isc,G%iec + ! It doesn't matter what these values are so make them zero. + kappa_io(i,j,1) = 0.0 + kappa_io(i,j,nz+1) = 0.0 + enddo + if (CS%VS_GeometricMean) then + do K=2,nz ; do i=G%isc,G%iec ! Set the diffusivities in tracer columns from the values at vertices. - kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * & - ((kappa_2d(I-1,K,J2m1)+kappa_2d(I,K,J2)) +& - (kappa_2d(I-1,K,J2)+kappa_2d(I,K,J2m1))) + ! The geometric mean is zero if any component is zero, hence the need + ! and sensitivity to a value of kappa_trunc in regions on boundaries of shear zones. + Kd_ll = max(kappa_vertex(I-1,J-1,K),CS%kappa_trunc) ! Diffusivity, lower left + Kd_ur = max(kappa_vertex(I,J,K),CS%kappa_trunc) ! Diffusivity, upper right + Kd_ul = max(kappa_vertex(I-1,J,K),CS%kappa_trunc) ! Diffusivity, upper left + Kd_lr = max(kappa_vertex(I,J-1,K),CS%kappa_trunc) ! Difusivity, lower right + if (CS%VS_ThicknessMean) then + weight_ll = hweight_ll(i,j,k) + weight_ur = hweight_ur(i,j,k) + weight_ul = hweight_ul(i,j,k) + weight_lr = hweight_lr(i,j,k) + endif + + ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: + kappa_io(i,j,K) = G%mask2dT(i,j) * & + ( ((Kd_ll**weight_ll)*(Kd_ur**weight_ur)) * & + ((Kd_ul**weight_ul)*(Kd_lr**weight_lr)) ) enddo; enddo - else - do i=G%isc,G%iec - ! It doesn't matter what these values are so make them zero. - kappa_io(i,1,j) = 0.0 - kappa_io(i,nz+1,j) = 0.0 - enddo - if (CS%VS_GeometricMean) then - do K=2,nz ; do i=G%isc,G%iec - ! Set the diffusivities in tracer columns from the values at vertices. - ! The geometric mean is zero if any component is zero, hence the need - ! and sensitivity to a value of kappa_trunc in regions on boundaries of shear zones. - a1 = max(kappa_2d(I-1,K,J2m1),CS%kappa_trunc) ! Diffusivity, lower left - a2 = max(kappa_2d(I,K,J2),CS%kappa_trunc) ! Diffusivity, upper right - a3 = max(kappa_2d(I-1,K,J2),CS%kappa_trunc) ! Diffusivity, upper left - a4 = max(kappa_2d(I,K,J2m1),CS%kappa_trunc) ! Difusivity, lower right - if (CS%VS_ThicknessMean) then - b1 = h_2d(I-1,k,J2m1)+h_2d(I-1,k-1,J2m1) ! 2x thickness, lower left - b2 = h_2d(I,k,J2)+h_2d(I,k-1,J2) ! 2x thickness, upper right - b3 = h_2d(I-1,k,J2)+h_2d(I-1,k-1,J2) ! 2x thickness, upper left - b4 = h_2d(I,k,J2m1)+h_2d(I,k-1,J2m1) ! 2x thickness, upper right - else - b1 = 1. - b2 = 1. - b3 = 1. - b4 = 1. - endif - bT = (b1 + b2) + (b3 + b4) - ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: - if (bT>GV%H_subroundoff) then - kappa_io(i,j,K) = G%mask2dT(i,j) * & - ( (a1**(b1/bT)*a2**(b2/bT)) * (a3**(b3/bT)*a4**(b4/bT)) ) - endif - enddo; enddo - else !arithmetic mean - do K=2,nz ; do i=G%isc,G%iec - a1 = kappa_2d(I-1,K,J2m1) ! Diffusivity, lower left - a2 = kappa_2d(I,K,J2) ! Diffusivity, upper right - a3 = kappa_2d(I-1,K,J2) ! Diffusivity, upper left - a4 = kappa_2d(I,K,J2m1) ! Difusivity, lower right - if (CS%VS_ThicknessMean) then - b1 = h_2d(I-1,k,J2m1)+h_2d(I-1,k-1,J2m1) ! 2x thickness, lower left - b2 = h_2d(I,k,J2)+h_2d(I,k-1,J2) ! 2x thickness, upper right - b3 = h_2d(I-1,k,J2)+h_2d(I-1,k-1,J2) ! 2x thickness, upper left - b4 = h_2d(I,k,J2m1)+h_2d(I,k-1,J2m1) ! 2x thickness, upper right - else - b1 = 1. - b2 = 1. - b3 = 1. - b4 = 1. - endif - bT = (b1+b2) + (b3+b4) - ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: - if (bT>GV%H_subroundoff) then - kappa_io(i,j,K) = G%mask2dT(i,j)*((a1*b1+a2*b2)+(a3*b3+a4*b4))/bT - endif - enddo; enddo - endif - endif - endif - enddo ! end of J-loop + else !arithmetic mean + do K=2,nz ; do i=G%isc,G%iec + Kd_ll = kappa_vertex(I-1,J-1,K) ! Diffusivity, lower left + Kd_ur = kappa_vertex(I,J,K) ! Diffusivity, upper right + Kd_ul = kappa_vertex(I-1,J,K) ! Diffusivity, upper left + Kd_lr = kappa_vertex(I,J-1,K) ! Difusivity, lower right + if (CS%VS_ThicknessMean) then + weight_ll = hweight_ll(i,j,k) + weight_ur = hweight_ur(i,j,k) + weight_ul = hweight_ul(i,j,k) + weight_lr = hweight_lr(i,j,k) + endif + ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: + kappa_io(i,j,K) = G%mask2dT(i,j) * & + ((Kd_ll*weight_ll + Kd_ur*weight_ur) + (Kd_ul*weight_ul + Kd_lr*weight_lr)) + enddo; enddo + endif ! geometric/arithmetic mean + endif ! Answer date + enddo ! j loop if (CS%debug) then call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s) @@ -825,7 +862,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag) if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag) - if (CS%id_Kd_vertex > 0) call post_data(CS%id_Kd_vertex, diag_Kd_vertex, CS%diag) + if (CS%id_Kd_vertex > 0) call post_data(CS%id_Kd_vertex, kappa_vertex, CS%diag) if (CS%id_N2_init > 0) call post_data(CS%id_N2_init, diag_N2_init, CS%diag) if (CS%id_S2_init > 0) call post_data(CS%id_S2_init, diag_S2_init, CS%diag) if (CS%id_N2_mean > 0) call post_data(CS%id_N2_mean, diag_N2_mean, CS%diag) From 9fa3d3d184597c7ede546d02b20f395ee176459d Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Wed, 5 Mar 2025 09:07:04 -0500 Subject: [PATCH 4/8] Initialize new variables to zero in kappa_shear --- src/parameterizations/vertical/MOM_kappa_shear.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 457e5d4590..2edb76a5ab 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -521,6 +521,11 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0 if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 kappa_vertex(:,:,:) = 0.0 + hweight_ul(:,:,:) = 0.0 + hweight_ur(:,:,:) = 0.0 + hweight_ll(:,:,:) = 0.0 + hweight_lr(:,:,:) = 0.0 + h_int_mask(:,:,:) = 0.0 use_temperature = associated(tv%T) From 14455e6f6a7a4e91de5358b26312991e527c555c Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Wed, 5 Mar 2025 09:16:10 -0500 Subject: [PATCH 5/8] Only initialize some arrays in kappa shear if they are used --- src/parameterizations/vertical/MOM_kappa_shear.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 2edb76a5ab..02a20b19b5 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -521,11 +521,6 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0 if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 kappa_vertex(:,:,:) = 0.0 - hweight_ul(:,:,:) = 0.0 - hweight_ur(:,:,:) = 0.0 - hweight_ll(:,:,:) = 0.0 - hweight_lr(:,:,:) = 0.0 - h_int_mask(:,:,:) = 0.0 use_temperature = associated(tv%T) @@ -537,6 +532,11 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ call thickness_to_dz(h, tv, dz_3d, G, GV, US, halo_size=1) if (CS%VS_ThicknessMean) then + hweight_ul(:,:,:) = 0.0 + hweight_ur(:,:,:) = 0.0 + hweight_ll(:,:,:) = 0.0 + hweight_lr(:,:,:) = 0.0 + h_int_mask(:,:,:) = 0.0 !Loop over halo and pre compute a masked thickness around each interface ! (excluding the top and bottom boundary) do j=G%jsc-1,G%jec+1; do K=2,nz; do i=G%isc-1,G%iec+1 From 71cb1351984c7fbf9d096d19fcd6e490f8e1a24f Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Wed, 5 Mar 2025 15:12:09 -0500 Subject: [PATCH 6/8] Add VERTEX_SHEAR_GEOMETRIC_MEAN_KDMIN runtime option to Kappa Shear --- .../vertical/MOM_kappa_shear.F90 | 28 ++++++++++++++----- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 02a20b19b5..8b53933eaa 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -90,6 +90,10 @@ module MOM_kappa_shear !! greater than 1. The lower limit for the permitted fractional !! decrease is (1 - 0.5/kappa_src_max_chg). These limits could !! perhaps be made dynamic with an improved iterative solver. + real :: VS_GeometricMean_Kd_min !< A minimum diffusivity for computing the horizontal averages + !! when using the geometric mean with VERTEX_SHEAR=True. The model + !! is sensitive to this value, which is a drawback of using the + !! geometric average as currently implemented. logical :: psurf_bug !< If true, do a simple average of the cell surface pressures to get a !! surface pressure at the corner if VERTEX_SHEAR=True. Otherwise mask !! out any land points in the average. @@ -824,10 +828,10 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Set the diffusivities in tracer columns from the values at vertices. ! The geometric mean is zero if any component is zero, hence the need ! and sensitivity to a value of kappa_trunc in regions on boundaries of shear zones. - Kd_ll = max(kappa_vertex(I-1,J-1,K),CS%kappa_trunc) ! Diffusivity, lower left - Kd_ur = max(kappa_vertex(I,J,K),CS%kappa_trunc) ! Diffusivity, upper right - Kd_ul = max(kappa_vertex(I-1,J,K),CS%kappa_trunc) ! Diffusivity, upper left - Kd_lr = max(kappa_vertex(I,J-1,K),CS%kappa_trunc) ! Difusivity, lower right + Kd_ll = max(kappa_vertex(I-1,J-1,K),CS%VS_GeometricMean_Kd_min) ! Diffusivity, lower left + Kd_ur = max(kappa_vertex(I,J,K),CS%VS_GeometricMean_Kd_min) ! Diffusivity, upper right + Kd_ul = max(kappa_vertex(I-1,J,K),CS%VS_GeometricMean_Kd_min) ! Diffusivity, upper left + Kd_lr = max(kappa_vertex(I,J-1,K),CS%VS_GeometricMean_Kd_min) ! Difusivity, lower right if (CS%VS_ThicknessMean) then weight_ll = hweight_ll(i,j,k) weight_ur = hweight_ur(i,j,k) @@ -2152,9 +2156,19 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) "Answer date for the algorithm for arithmetic horizontal mean of "//& "Kd from vertices to tracer points. Less than 20250304 uses an older algorithm.", & default=20250303, do_not_log=just_read.or.(.not.CS%KS_at_vertex)) - if ((CS%Kd_tracer_mean_answer_date<20250304).and.(CS%VS_ThicknessMean .or. CS%VS_GeometricMean)) then - call MOM_error(FATAL,"Cannot use VERTEX_SHEAR_THICKNESS_MEAN or VERTEX_SHEAR_GEOMETRIC_MEAN "//& - "with VERTEX_SHEAR_KD_MEAN_ANSWER_DATE < 20250304") + if (CS%Kd_tracer_mean_answer_date<20250304) then + if (CS%VS_ThicknessMean .or. CS%VS_GeometricMean) then + call MOM_error(FATAL,"Cannot use VERTEX_SHEAR_THICKNESS_MEAN or VERTEX_SHEAR_GEOMETRIC_MEAN "//& + "with VERTEX_SHEAR_KD_MEAN_ANSWER_DATE < 20250304") + endif + else + if (CS%VS_GeometricMean) then + call get_param(param_file, mdl, "VERTEX_SHEAR_GEOMETRIC_MEAN_KDMIN", & + CS%VS_GeometricMean_Kd_min, "If using the geometric mean in vertex shear, "//& + "use this minimum value for Kd. This is an ad-hoc parameter, the "//& + "diffusivities on the edge of shear regions are sensitive to the choice.",& + units="m2 s-1",default=0.0, scale=GV%m2_s_to_HZ_T, do_not_log=just_read) + endif endif call get_param(param_file, mdl, "RINO_CRIT", CS%RiNo_crit, & "The critical Richardson number for shear mixing.", & From 8eaf06a543a198d0252bc26f00fc2fa3b9d5f3d1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 5 Mar 2025 18:06:49 -0500 Subject: [PATCH 7/8] *+Revise horizontal avg in Calc_kappa_shear_vertex Revised the Calc_kappa_shear_vertex code where it averages the diffusivities at the vertices back to the tracer points when there is either a geometric mean being used or thickness weighting. With these changes, the code should now pass dimensional consistency testing, work properly with openMP threading enabled, and give sensible values of the diffusivity for massless layers. Also corrected the indenting of several recently added lines to follow the MOM6 2-space indenting convention. This commit can change answers slightly in some cases where VERTEX_SHEAR = True, especially if VERTEX_SHEAR_GEOMETRIC_MEAN = True or VERTEX_SHEAR_THICKNESS_MEAN = True. The recently added runtime parameter VERTEX_SHEAR_KD_MEAN_ANSWER_DATE (which has not yet been added to the dev/gfdl branch of the code) was no longer needed and has been eliminated. --- .../vertical/MOM_kappa_shear.F90 | 295 +++++++----------- 1 file changed, 111 insertions(+), 184 deletions(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 02a20b19b5..307f0acbad 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -97,12 +97,11 @@ module MOM_kappa_shear !! time average TKE when there is mass in all layers. Otherwise always !! report the time-averaged TKE, as is currently done when there !! are some massless layers. - integer :: Kd_tracer_mean_answer_date !< Answer date for Kd tracer point horizontal averaging logical :: VS_viscosity_bug !< If true, use a bug in the calculation of the viscosity that sets !! it to zero for all vertices that are on a coastline. logical :: VS_GeometricMean !< If true use geometric averaging for Kd from vertices to tracer points logical :: VS_ThicknessMean !< If true use thickness weighting when averaging Kd from vertices to - !!tracer points + !! tracer points logical :: restrictive_tolerance_check !< If false, uses the less restrictive tolerance check to !! determine if a timestep is acceptable for the KS_it outer iteration !! loop, as the code was originally written. True uses the more @@ -349,10 +348,8 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & kappa_2d(i,K) = kappa_avg(kc(K)) tke_2d(i,K) = tke_avg(kc(K)) else - kappa_2d(i,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + & - kf(K) * kappa_avg(kc(K)+1) - tke_2d(i,K) = (1.0-kf(K)) * tke_avg(kc(K)) + & - kf(K) * tke_avg(kc(K)+1) + kappa_2d(i,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1) + tke_2d(i,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1) endif enddo do K=1,nz+1 @@ -362,22 +359,14 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & if (CS%id_N2_init>0) diag_N2_init(i,j,K) = N2_init(kc(K)) if (CS%id_S2_init>0) diag_S2_init(i,j,K) = S2_init(kc(K)) else - if (CS%id_N2_mean>0) then - diag_N2_mean(i,j,K) = (1.0-kf(K)) * N2_mean(kc(K)) & - + kf(K) *N2_mean(kc(K)+1) - endif - if (CS%id_S2_mean>0) then - diag_S2_mean(i,j,K) = (1.0-kf(K)) * S2_mean(kc(K)) & - + kf(K) *S2_mean(kc(K)+1) - endif - if (CS%id_N2_init>0) then - diag_N2_init(i,j,K) = (1.0-kf(K)) * N2_init(kc(K)) & - + kf(K) *N2_init(kc(K)+1) - endif - if (CS%id_S2_init>0) then - diag_S2_init(i,j,K) = (1.0-kf(K)) * S2_init(kc(K)) & - + kf(K) *S2_init(kc(K)+1) - endif + if (CS%id_N2_mean>0) & + diag_N2_mean(i,j,K) = (1.0-kf(K)) * N2_mean(kc(K)) + kf(K) * N2_mean(kc(K)+1) + if (CS%id_S2_mean>0) & + diag_S2_mean(i,j,K) = (1.0-kf(K)) * S2_mean(kc(K)) + kf(K) * S2_mean(kc(K)+1) + if (CS%id_N2_init>0) & + diag_N2_init(i,j,K) = (1.0-kf(K)) * N2_init(kc(K)) + kf(K) * N2_init(kc(K)+1) + if (CS%id_S2_init>0) & + diag_S2_init(i,j,K) = (1.0-kf(K)) * S2_init(kc(K)) + kf(K) * S2_init(kc(K)+1) endif enddo endif @@ -457,20 +446,18 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ diag_S2_mean ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & - dz_3d ! Vertical distance between interface heights [Z ~> m]. + dz_3d ! Vertical distance between interface heights [Z ~> m]. + real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: & + kappa_vertex ! Diffusivity at interfaces and vertices [H Z T-1 ~> m2 s-1 or Pa s] + real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: & + h_vert ! Thicknesses interpolated to vertices [H Z T-1 ~> m2 s-1 or Pa s] real, dimension(SZIB_(G),SZK_(GV)) :: & - h_2d, & ! A 2-D version of h [H ~> m or kg m-2]. + h_2d, & ! A 2-D version of h interpolated to vertices [H ~> m or kg m-2]. dz_2d, & ! Vertical distance between interface heights [Z ~> m]. u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1]. T_2d, S_2d, rho_2d ! 2-D versions of T [C ~> degC], S [S ~> ppt], and rho [R ~> kg m-3]. - real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: & - kappa_vertex, & ! kappa_io on vertices [H Z T-1 ~> m2 s-1 or Pa s] - h_int_mask ! masked h sumed over bounding layers on vertices [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & - hweight_ul, & ! thickness weighting at cell center from upper-left vertex - hweight_ur, & ! thickness weighting at cell center from upper-right vertex - hweight_ll, & ! thickness weighting at cell center from lower-left vertex - hweight_lr ! thickness weighting at cell center from lower-right vertex + real, dimension(SZIB_(G),SZK_(GV)+1) :: & + kappa_2d ! 2-D slice of kappa_vert [H Z T-1 ~> m2 s-1 or Pa s] real, dimension(SZIB_(G),SZK_(GV)+1) :: & tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2]. real, dimension(SZK_(GV)) :: & @@ -497,9 +484,9 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ real :: dz_in_lay ! The running sum of the thickness in a layer [H ~> m or kg m-2] real :: k0dt ! The background diffusivity times the timestep [H Z ~> m2 or kg m-1] real :: dz_massless ! A layer thickness that is considered massless [H ~> m or kg m-2] - real :: I_hwt ! The inverse of the masked thickness weights [H-1 ~> m-1 or m2 kg-1]. + real :: I_hwt ! The inverse of the sum of the adjacent masked thickness weights [H-1 ~> m-1 or m2 kg-1] + real :: I_htot ! The inverse of the sum of the thicknesses at adjacent vertices [H-1 ~> m-1 or m2 kg-1] real :: I_Prandtl ! The inverse of the turbulent Prandtl number [nondim]. - real :: hweight_denom ! A weighting factor for thickness averaging [H ~> m or kg m-2] logical :: use_temperature ! If true, temperature and salinity have been ! allocated and are being used as state variables. @@ -508,9 +495,10 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! merged into nearby massive layers. real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for ! interpolating back to the original index space [nondim]. - ! Terms in horizontal averaging of diffusivity to tracer location - real :: Kd_ll, Kd_lr, Kd_ul, Kd_ur, weight_ll, weight_lr, weight_ul, weight_ur - ! Other useful indices + real :: h_SW, h_SE, h_NW, h_NE ! Thicknesses at adjacent vertices [H ~> m or kg m-2] + real :: mks_to_HZ_T ! A factor used to restore dimensional scaling after the geomentric mean + ! diffusivity is taken using thickness weighted powers [H Z s m-2 T-1 ~> 1] + ! or [H Z m s kg-1 T-1 ~> 1] integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc ! Diagnostics that should be deleted? @@ -531,61 +519,8 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Convert layer thicknesses into geometric thickness in height units. call thickness_to_dz(h, tv, dz_3d, G, GV, US, halo_size=1) - if (CS%VS_ThicknessMean) then - hweight_ul(:,:,:) = 0.0 - hweight_ur(:,:,:) = 0.0 - hweight_ll(:,:,:) = 0.0 - hweight_lr(:,:,:) = 0.0 - h_int_mask(:,:,:) = 0.0 - !Loop over halo and pre compute a masked thickness around each interface - ! (excluding the top and bottom boundary) - do j=G%jsc-1,G%jec+1; do K=2,nz; do i=G%isc-1,G%iec+1 - h_int_mask(i,j,K) = G%mask2dT(i,j)*(h(i,j,k-1)+h(i,j,k))/2. - enddo; enddo ; enddo - !Compute the thickness weighting if requested for horizontal mean - do J=JsB,JeB ; do K=2,nz ; do I=IsB,IeB - hweight_denom = h_int_mask(i-1,j-1,k) + & !left,bot - 2.*h_int_mask(i-1,j,k) + & !left,cen - h_int_mask(i-1,j+1,k) + & !left,top - 2. * h_int_mask(i,j-1,k) + &!mid,bot - 4. * h_int_mask(i,j,k) + & !mid,cen - 2. * h_int_mask(i,j+1,k) + &!mid,top - h_int_mask(i+1,j-1,k) + & !right,bot - 2. * h_int_mask(i-1,j,k) + &!right,cen - h_int_mask(i+1,j+1,k) + & !right,top - GV%H_subroundoff - ! Computing the weighting of the 4 vertices to each cell center - hweight_ul(I,J,K) = ((h_int_mask(i-1,j,k) + & !left,cen - h_int_mask(i-1,j+1,k)) + &!left,top - (h_int_mask(i,j,k) + & !mid,cen - h_int_mask(i,j+1,k)) & !mid,top - ) / hweight_denom - hweight_ur(I,J,k) = ((h_int_mask(i,j,k) + & !mid,cen - h_int_mask(i,j+1,k)) + &!mid,top - (h_int_mask(i+1,j,k) + & !right,cen - h_int_mask(i+1,j+1,k)) &!right,top - ) / hweight_denom - hweight_ll(I,J,k) = ((h_int_mask(i-1,j-1,k) + &!left,bot - h_int_mask(i-1,j,k)) + & !left,cen - (h_int_mask(i+1,j-1,k) + &!mid,bot - h_int_mask(i+1,j,k)) & !mid,cen - ) / hweight_denom - hweight_lr(I,J,k) = ((h_int_mask(i,j-1,k) + & !mid,bot - h_int_mask(i,j,k)) + & !mid,cen - (h_int_mask(i+1,j-1,k) + &!right,bot - h_int_mask(i+1,j,k)) & !right,cen - ) / hweight_denom - enddo ; enddo ; enddo - else - ! Set some weights for horizontal averaging factors - weight_ll=0.25 - weight_lr=0.25 - weight_ul=0.25 - weight_ur=0.25 - endif - - !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,tv,G,GV, & - !$OMP US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl, & + !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,tv,G,GV,US,CS,kappa_io, & + !$OMP dz_massless,k0dt,p_surf,dt,tke_io,kv_io,kappa_vertex,h_vert,I_Prandtl, & !$OMP diag_N2_init,diag_S2_init,diag_N2_mean,diag_S2_mean) do J=JsB,JeB @@ -727,7 +662,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Extrapolate from the vertically reduced grid back to the original layers. if (nz == nzc) then do K=1,nz+1 - kappa_vertex(I,J,K) = kappa_avg(K) + kappa_2d(I,K) = kappa_avg(K) if (CS%all_layer_TKE_bug) then tke_2d(i,K) = tke(K) else @@ -749,10 +684,10 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ else do K=1,nz+1 if (kf(K) == 0.0) then - kappa_vertex(I,J,K) = kappa_avg(kc(K)) + kappa_2d(I,K) = kappa_avg(kc(K)) tke_2d(I,K) = tke_avg(kc(K)) else - kappa_vertex(I,J,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1) + kappa_2d(I,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1) tke_2d(I,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1) endif enddo @@ -763,102 +698,102 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%id_N2_init>0) diag_N2_init(I,J,K) = N2_init(kc(K)) if (CS%id_S2_init>0) diag_S2_init(I,J,K) = S2_init(kc(K)) else - if (CS%id_N2_mean>0) then - diag_N2_mean(I,J,K) = (1.0-kf(K)) * N2_mean(kc(K)) & - + kf(K) *N2_mean(kc(K)+1) - endif - if (CS%id_S2_mean>0) then - diag_S2_mean(I,J,K) = (1.0-kf(K)) * S2_mean(kc(K)) & - + kf(K) *S2_mean(kc(K)+1) - endif - if (CS%id_N2_init>0) then - diag_N2_init(I,J,K) = (1.0-kf(K)) * N2_init(kc(K)) & - + kf(K) *N2_init(kc(K)+1) - endif - if (CS%id_S2_init>0) then - diag_S2_init(I,J,K) = (1.0-kf(K)) * S2_init(kc(K)) & - + kf(K) *S2_init(kc(K)+1) - endif + if (CS%id_N2_mean>0) & + diag_N2_mean(I,J,K) = (1.0-kf(K)) * N2_mean(kc(K)) + kf(K) * N2_mean(kc(K)+1) + if (CS%id_S2_mean>0) & + diag_S2_mean(I,J,K) = (1.0-kf(K)) * S2_mean(kc(K)) + kf(K) * S2_mean(kc(K)+1) + if (CS%id_N2_init>0) & + diag_N2_init(I,J,K) = (1.0-kf(K)) * N2_init(kc(K)) + kf(K) * N2_init(kc(K)+1) + if (CS%id_S2_init>0) & + diag_S2_init(I,J,K) = (1.0-kf(K)) * S2_init(kc(K)) + kf(K) * S2_init(kc(K)+1) endif enddo endif ! call cpu_clock_end(Id_clock_setup) else ! Land points, still inside the i-loop. do K=1,nz+1 - kappa_vertex(I,J,K) = 0.0 ; tke_2d(I,K) = 0.0 + kappa_2d(I,K) = 0.0 ; tke_2d(I,K) = 0.0 enddo endif ; enddo ! i-loop + ! Store the 2-d slices back in the 3-d arrays for restarts or interpolation back to tracer points. + if (CS%VS_ThicknessMean) then + do K=1,nz+1 ; do I=IsB,IeB + h_vert(I,J,k) = h_2d(I,k) + enddo ; enddo + endif if (CS%VS_viscosity_bug) then do K=1,nz+1 ; do I=IsB,IeB + kappa_vertex(I,J,K) = kappa_2d(I,K) tke_io(I,J,K) = G%mask2dBu(I,J) * tke_2d(I,K) kv_io(I,J,K) = ( G%mask2dBu(I,J) * kappa_vertex(I,J,K) ) * CS%Prandtl_turb - enddo; enddo + enddo ; enddo else do K=1,nz+1 ; do I=IsB,IeB + kappa_vertex(I,J,K) = kappa_2d(I,K) tke_io(I,J,K) = tke_2d(I,K) kv_io(I,J,K) = kappa_vertex(I,J,K) * CS%Prandtl_turb - enddo; enddo + enddo ; enddo endif enddo ! end of J-loop + ! Set the diffusivities in tracer columns from the values at vertices. - do j=G%jsc,G%jec - if (CS%Kd_tracer_mean_answer_date < 20250304) then - ! This is a non-thickness weighted arithmetic mean, but it is not bitwise identical to the more general - ! version of the mean below the "else" that includes a mathematically equivalent form. - do K=1,nz+1 ; do i=G%isc,G%iec - ! Set the diffusivities in tracer columns from the values at vertices. - kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * & - ((kappa_vertex(I-1,J-1,K)+kappa_vertex(I,J,K)) +& - (kappa_vertex(I-1,J,K)+kappa_vertex(I,J-1,K))) - enddo; enddo - else - do i=G%isc,G%iec - ! It doesn't matter what these values are so make them zero. - kappa_io(i,j,1) = 0.0 - kappa_io(i,j,nz+1) = 0.0 - enddo + !$OMP parallel do default(private) shared(G,kappa_io) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! The turbulent length scales (and hence turbulent diffusivity) should always go to 0 at the top and bottom. + kappa_io(i,j,1) = 0.0 + kappa_io(i,j,nz+1) = 0.0 + enddo ; enddo + if (CS%VS_ThicknessMean) then + ! This conversion factor is required to allow for aribtrary fracional powers of the diffusivities. + if (CS%VS_GeometricMean) mks_to_HZ_T = 1.0 / GV%HZ_T_to_MKS + !$OMP parallel do default(private) shared(nz,G,GV,CS,kappa_io,kappa_vertex,h_vert) + do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + h_SW = 0.5 * (h_vert(I-1,J-1,k) + h_vert(I-1,J-1,k-1)) + h_NE = 0.5 * (h_vert(I,J,k) + h_vert(I,J,k-1)) + h_NW = 0.5 * (h_vert(I-1,J,k) + h_vert(I-1,J,k-1)) + h_SE = 0.5 * (h_vert(I,J-1,k) + h_vert(I,J-1,k-1)) if (CS%VS_GeometricMean) then - do K=2,nz ; do i=G%isc,G%iec - ! Set the diffusivities in tracer columns from the values at vertices. - ! The geometric mean is zero if any component is zero, hence the need - ! and sensitivity to a value of kappa_trunc in regions on boundaries of shear zones. - Kd_ll = max(kappa_vertex(I-1,J-1,K),CS%kappa_trunc) ! Diffusivity, lower left - Kd_ur = max(kappa_vertex(I,J,K),CS%kappa_trunc) ! Diffusivity, upper right - Kd_ul = max(kappa_vertex(I-1,J,K),CS%kappa_trunc) ! Diffusivity, upper left - Kd_lr = max(kappa_vertex(I,J-1,K),CS%kappa_trunc) ! Difusivity, lower right - if (CS%VS_ThicknessMean) then - weight_ll = hweight_ll(i,j,k) - weight_ur = hweight_ur(i,j,k) - weight_ul = hweight_ul(i,j,k) - weight_lr = hweight_lr(i,j,k) - endif - - ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: - kappa_io(i,j,K) = G%mask2dT(i,j) * & - ( ((Kd_ll**weight_ll)*(Kd_ur**weight_ur)) * & - ((Kd_ul**weight_ul)*(Kd_lr**weight_lr)) ) - enddo; enddo - else !arithmetic mean - do K=2,nz ; do i=G%isc,G%iec - Kd_ll = kappa_vertex(I-1,J-1,K) ! Diffusivity, lower left - Kd_ur = kappa_vertex(I,J,K) ! Diffusivity, upper right - Kd_ul = kappa_vertex(I-1,J,K) ! Diffusivity, upper left - Kd_lr = kappa_vertex(I,J-1,K) ! Difusivity, lower right - if (CS%VS_ThicknessMean) then - weight_ll = hweight_ll(i,j,k) - weight_ur = hweight_ur(i,j,k) - weight_ul = hweight_ul(i,j,k) - weight_lr = hweight_lr(i,j,k) - endif - ! The following expression generalizes the geometric mean at tracer points, w/ or w/o thickness weighting: - kappa_io(i,j,K) = G%mask2dT(i,j) * & - ((Kd_ll*weight_ll + Kd_ur*weight_ur) + (Kd_ul*weight_ul + Kd_lr*weight_lr)) - enddo; enddo - endif ! geometric/arithmetic mean - endif ! Answer date - enddo ! j loop + if ((h_SW + h_NE) + (h_NW + h_SE) > 0.0) then + ! The geometric mean is zero if any component is zero, hence the need to use a floor + ! on the value of kappa_trunc in regions on boundaries of shear zones. + I_htot = 1.0 / ((h_SW + h_NE) + (h_NW + h_SE)) + kappa_io(i,j,K) = G%mask2dT(i,j) * mks_to_HZ_T * & + ( ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J-1,K),CS%kappa_trunc))**(h_SW*I_htot) * & + (GV%HZ_T_to_MKS * max(kappa_vertex(I,J,K),CS%kappa_trunc))**(h_NE*I_htot)) * & + ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J,K),CS%kappa_trunc))**(h_NW*I_htot) * & + (GV%HZ_T_to_MKS * max(kappa_vertex(I,J-1,K),CS%kappa_trunc))**(h_SE*I_htot)) ) + else + ! If all points have zero thickness, the thikncess-weighted geometric mean is undefined, so use + ! the non-thickness weighted geometric mean instead. + kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( & + (max(kappa_vertex(I-1,J-1,K),CS%kappa_trunc) * max(kappa_vertex(I,J,K),CS%kappa_trunc)) * & + (max(kappa_vertex(I-1,J,K),CS%kappa_trunc) * max(kappa_vertex(I,J-1,K),CS%kappa_trunc)) )) + endif + else + ! The following expression is a thickness weighted arithmetic mean at tracer points: + I_htot = 1.0 / (((h_SW + h_NE) + (h_NW + h_SE)) + GV%H_subroundoff) + kappa_io(i,j,K) = G%mask2dT(i,j) * & + (((kappa_vertex(I-1,J-1,K)*h_SW) + (kappa_vertex(I,J,K)*h_NE)) + & + ((kappa_vertex(I-1,J,K)*h_NW) + (kappa_vertex(I,J-1,K)*h_SE))) * I_htot + endif + enddo ; enddo ; enddo + elseif (CS%VS_GeometricMean) then ! The geometic mean diffusivities are not thickness weighted. + !$OMP parallel do default(private) shared(nz,G,CS,kappa_io,kappa_vertex) + do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( & + (max(kappa_vertex(I-1,J-1,K),CS%kappa_trunc) * max(kappa_vertex(I,J,K),CS%kappa_trunc)) * & + (max(kappa_vertex(I-1,J,K),CS%kappa_trunc) * max(kappa_vertex(I,J-1,K),CS%kappa_trunc)) )) + enddo ; enddo ; enddo + else ! Use a non-thickness weighted arithmetic mean. + !$OMP parallel do default(private) shared(nz,G,CS,kappa_io,kappa_vertex) + do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * & + ((kappa_vertex(I-1,J-1,K) + kappa_vertex(I,J,K)) +& + (kappa_vertex(I-1,J,K) + kappa_vertex(I,J-1,K))) + enddo ; enddo ; enddo + endif if (CS%debug) then call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s) @@ -2148,14 +2083,6 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) "If true, apply thickness weighting to horizontal averagings of diffusivity "//& "to tracer points in the kappa shear solver.", & default=.false.) - call get_param(param_file, mdl, "VERTEX_SHEAR_KD_MEAN_ANSWER_DATE", CS%Kd_tracer_mean_answer_date, & - "Answer date for the algorithm for arithmetic horizontal mean of "//& - "Kd from vertices to tracer points. Less than 20250304 uses an older algorithm.", & - default=20250303, do_not_log=just_read.or.(.not.CS%KS_at_vertex)) - if ((CS%Kd_tracer_mean_answer_date<20250304).and.(CS%VS_ThicknessMean .or. CS%VS_GeometricMean)) then - call MOM_error(FATAL,"Cannot use VERTEX_SHEAR_THICKNESS_MEAN or VERTEX_SHEAR_GEOMETRIC_MEAN "//& - "with VERTEX_SHEAR_KD_MEAN_ANSWER_DATE < 20250304") - endif call get_param(param_file, mdl, "RINO_CRIT", CS%RiNo_crit, & "The critical Richardson number for shear mixing.", & units="nondim", default=0.25, do_not_log=just_read) @@ -2305,27 +2232,27 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) CS%id_S2_init = register_diag_field('ocean_model','S2_shear_in', diag%axesBi, Time, & 'Interface shear squared at horizontal vertices, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) CS%id_N2_init = register_diag_field('ocean_model','N2_shear_in', diag%axesBi, Time, & - 'Interface straitification at horizontal vertices, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) + 'Interface stratification at horizontal vertices, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) CS%id_S2_mean = register_diag_field('ocean_model','S2_shear_mean', diag%axesBi, Time, & 'Interface shear squared at horizontal vertices, averaged over timestep in kappa-shear', & 's-2', conversion=US%s_to_T**2) CS%id_N2_mean = register_diag_field('ocean_model','N2_shear_mean', diag%axesBi, Time, & - 'Interface straitification at horizontal vertices, averaged over timestep in kappa-shear', & + 'Interface stratification at horizontal vertices, averaged over timestep in kappa-shear', & 's-2', conversion=US%s_to_T**2) else CS%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesTi, Time, & - 'Shear-driven Turbulent Kinetic Energy at horizontal tracer points', 'm2 s-2', & - conversion=US%Z_to_m**2*US%s_to_T**2) + 'Shear-driven Turbulent Kinetic Energy at horizontal tracer points', & + 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) CS%id_S2_init = register_diag_field('ocean_model','S2_shear_in', diag%axesTi, Time, & 'Interface shear squared at horizontal tracer points, as input to kappa-shear', 's-2', conversion=US%s_to_T**2) CS%id_N2_init = register_diag_field('ocean_model','N2_shear_in', diag%axesTi, Time, & - 'Interface straitification at horizontal tracer points, as input to kappa-shear', 's-2', & - conversion=US%s_to_T**2) + 'Interface stratification at horizontal tracer points, as input to kappa-shear', & + 's-2', conversion=US%s_to_T**2) CS%id_S2_mean = register_diag_field('ocean_model','S2_shear_mean', diag%axesTi, Time, & 'Interface shear squared at horizontal tracer points, averaged over timestep in kappa-shear', & 's-2', conversion=US%s_to_T**2) CS%id_N2_mean = register_diag_field('ocean_model','N2_shear_mean', diag%axesTi, Time, & - 'Interface straitification at horizontal tracer points, averaged ove timestep in kappa-shear', & + 'Interface stratification at horizontal tracer points, averaged ove timestep in kappa-shear', & 's-2', conversion=US%s_to_T**2) endif From 1e3742157ff69be42af8d926178aeadd568b19b8 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Thu, 6 Mar 2025 07:57:38 -0500 Subject: [PATCH 8/8] Shorten variable name VS_GeometricMean_Kd_min to VS_GeoMean_Kdmin to help with line length issues. --- .../vertical/MOM_kappa_shear.F90 | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 57111324f0..53d6b36e4a 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -90,7 +90,7 @@ module MOM_kappa_shear !! greater than 1. The lower limit for the permitted fractional !! decrease is (1 - 0.5/kappa_src_max_chg). These limits could !! perhaps be made dynamic with an improved iterative solver. - real :: VS_GeometricMean_Kd_min !< A minimum diffusivity for computing the horizontal averages + real :: VS_GeoMean_Kdmin !< A minimum diffusivity for computing the horizontal averages !! when using the geometric mean with VERTEX_SHEAR=True. The model !! is sensitive to this value, which is a drawback of using the !! geometric average as currently implemented. @@ -764,16 +764,16 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! on the value of kappa_trunc in regions on boundaries of shear zones. I_htot = 1.0 / ((h_SW + h_NE) + (h_NW + h_SE)) kappa_io(i,j,K) = G%mask2dT(i,j) * mks_to_HZ_T * & - ( ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J-1,K),CS%VS_GeometricMean_Kd_min))**(h_SW*I_htot) * & - (GV%HZ_T_to_MKS * max(kappa_vertex(I,J,K),CS%VS_GeometricMean_Kd_min))**(h_NE*I_htot)) * & - ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J,K),CS%VS_GeometricMean_Kd_min))**(h_NW*I_htot) * & - (GV%HZ_T_to_MKS * max(kappa_vertex(I,J-1,K),CS%VS_GeometricMean_Kd_min))**(h_SE*I_htot)) ) + ( ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin))**(h_SW*I_htot) * & + (GV%HZ_T_to_MKS * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin))**(h_NE*I_htot)) * & + ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin))**(h_NW*I_htot) * & + (GV%HZ_T_to_MKS * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin))**(h_SE*I_htot)) ) else ! If all points have zero thickness, the thikncess-weighted geometric mean is undefined, so use ! the non-thickness weighted geometric mean instead. kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( & - (max(kappa_vertex(I-1,J-1,K),CS%VS_GeometricMean_Kd_min) * max(kappa_vertex(I,J,K),CS%VS_GeometricMean_Kd_min)) * & - (max(kappa_vertex(I-1,J,K),CS%VS_GeometricMean_Kd_min) * max(kappa_vertex(I,J-1,K),CS%VS_GeometricMean_Kd_min)) )) + (max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin)) * & + (max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin)) )) endif else ! The following expression is a thickness weighted arithmetic mean at tracer points: @@ -787,8 +787,8 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ !$OMP parallel do default(private) shared(nz,G,CS,kappa_io,kappa_vertex) do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( & - (max(kappa_vertex(I-1,J-1,K),CS%VS_GeometricMean_Kd_min) * max(kappa_vertex(I,J,K),CS%VS_GeometricMean_Kd_min)) * & - (max(kappa_vertex(I-1,J,K),CS%VS_GeometricMean_Kd_min) * max(kappa_vertex(I,J-1,K),CS%VS_GeometricMean_Kd_min)) )) + (max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin)) * & + (max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin)) )) enddo ; enddo ; enddo else ! Use a non-thickness weighted arithmetic mean. !$OMP parallel do default(private) shared(nz,G,CS,kappa_io,kappa_vertex) @@ -2089,7 +2089,7 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) default=.false.) if (CS%VS_GeometricMean) then call get_param(param_file, mdl, "VERTEX_SHEAR_GEOMETRIC_MEAN_KDMIN", & - CS%VS_GeometricMean_Kd_min, "If using the geometric mean in vertex shear, "//& + CS%VS_GeoMean_Kdmin, "If using the geometric mean in vertex shear, "//& "use this minimum value for Kd. This is an ad-hoc parameter, the "//& "diffusivities on the edge of shear regions are sensitive to the choice.",& units="m2 s-1",default=0.0, scale=GV%m2_s_to_HZ_T, do_not_log=just_read)