From 0daa5457001ce76b756e6c62bf38e3716e4c768e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 2 Oct 2024 15:09:50 -0600 Subject: [PATCH 1/3] Stokes MOST modifications to unresolved shear * Add a control of the Stokes enhancement to the turbulent scales contributing to the unresolved shear (Stokes_Vt = 0 for no enhancement; Stokes_Vt = Stokes_Xi for full enhancement); * Make a depth restriction to the turbulent velocity scale contribution to the unresolved shear, consistent with the contribution to diffusivity and viscosity (Vt_layer = 1 use full boundary layer depth; vt_layer = surf_layer_ext only includes the surface layer). Stokes MOST parameters are based on Vt_layer = 1; * Replace divisions with reciprocal multiplications. --- .../vertical/MOM_CVMix_KPP.F90 | 39 ++++++++++--------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 816d5d7498..c58178a977 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1028,6 +1028,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m] ! For Langmuir Calculations + real :: Vt_layer ! non-dimensional extent contribution to unresolved shear real :: LangEnhVt2 ! Langmuir enhancement for unresolved shear [nondim] real, dimension(GV%ke) :: U_H, V_H ! Velocities at tracer points [L T-1 ~> m s-1] real :: MLD_guess ! A guess at the mixed layer depth for calculating the Langmuir number [Z ~> m] @@ -1092,13 +1093,9 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl do k=1,GV%ke U_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)) V_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)) + uE_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)-Waves%US_x(I,j,k)-Waves%US_x(I-1,j,k)) + vE_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)-Waves%US_y(i,J,k)-Waves%US_y(i,J-1,k)) enddo - if (CS%StokesMOST) then - do k=1,GV%ke - uE_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)-Waves%US_x(I,j,k)-Waves%US_x(I-1,j,k)) - vE_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)-Waves%US_y(i,J,k)-Waves%US_y(i,J-1,k)) - enddo - endif ! things independent of position within the column Coriolis = 0.25*US%s_to_T*( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) + & @@ -1149,7 +1146,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl surfFricVel,waves%omega_w2x(i,j), uE_H, vE_H, uS_Hi, vS_Hi, uSbar_H, vSbar_H, uS_SLD,& vS_SLD, uSbar_SLD, vSbar_SLD, StokesXI, CVMix_kpp_params_user=CS%KPP_params ) StokesXI_1d(k) = StokesXI - StokesVt_1d(k) = StokesXI + StokesVt_1d(k) = 0.0 ! StokesXI ! average temperature, salinity, u and v over surface layer starting at ksfc delH = SLdepth_0d + iFaceHeight(ksfc-1) @@ -1175,6 +1172,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl Vk = vE_H(k) + vS_H(k) - surfV else !not StokesMOST + uS_H(k) = 0.5*(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) + vS_H(k) = 0.5*(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) StokesXI_1d(k) = 0.0 ! average temperature, salinity, u and v over surface layer @@ -1205,12 +1204,13 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl endif enddo - surfTemp = surfHtemp / hTot - surfSalt = surfHsalt / hTot - surfU = surfHu / hTot - surfV = surfHv / hTot - surfUs = surfHus / hTot - surfVs = surfHvs / hTot + I_hTot = 1./hTot + surfTemp = surfHtemp * I_hTot + surfSalt = surfHsalt * I_hTot + surfU = surfHu * I_hTot + surfV = surfHv * I_hTot + surfUs = surfHus * I_hTot + surfVs = surfHvs * I_hTot ! vertical shear between present layer and surface layer averaged surfU and surfV. ! C-grid average to get Uk and Vk on T-points. @@ -1296,12 +1296,13 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl z_inter(K) = US%Z_to_m*iFaceHeight(K) enddo - ! turbulent velocity scales w_s and w_m computed at the cell centers. - ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales - ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass - ! sigma=CS%surf_layer_ext for this calculation. + ! CVMix_kpp_compute_turbulent_scales_1d_OBL computes w_s velocity scale at cell centers for + ! CVmix_kpp_compute_bulk_Richardson call to CVmix_kpp_compute_unresolved_shear + ! at sigma=Vt_layer (CS%surf_layer_ext or 1.0) for this calculation. + ! StokesVt_1d controls Stokes enhancement (= 0 for none) + Vt_layer = 1.0 ! CS%surf_layer_ext call CVMix_kpp_compute_turbulent_scales( & ! 1d_OBL - CS%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext + Vt_layer, & ! (in) Boundary layer extent contributing to unresolved shear OBL_depth, & ! (in) OBL depth [m] surfBuoyFlux2, & ! (in) Buoyancy flux at surface [m2 s-3] surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] @@ -1341,7 +1342,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl N_iface=N_col, & ! Buoyancy frequency [s-1] EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim] LaSL=CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim] - bfsfc=surfBuoyFlux2, & ! surface buoyancy flux [m2 s-3] + bfsfc=surfBuoyFlux2, & ! surface buoyancy flux [m2 s-3] uStar=surfFricVel, & ! surface friction velocity [m s-1] CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters From d61e4e7fe706b95354c50986513a708ab59d6d9d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 11 Oct 2024 08:36:30 -0600 Subject: [PATCH 2/3] Include Vt_layer in the OMP shared clause --- src/parameterizations/vertical/MOM_CVMix_KPP.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index c58178a977..0c1d40e7ec 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1086,7 +1086,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl !$OMP uS_SLD, vS_SLD, uS_SLC, vS_SLC, uSbar_SLD, vSbar_SLD, & !$OMP StokesXI, StokesXI_1d, StokesVt_1d, kbl) & !$OMP shared(G, GV, CS, US, uStar, h, dz, buoy_scale, buoyFlux, & - !$OMP Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult) + !$OMP Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult, Vt_layer) do j = G%jsc, G%jec do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then From 574e6e58618c0106413656112d979634e1a9fa92 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 16 Oct 2024 16:06:46 -0600 Subject: [PATCH 3/3] Revert * I_hTot and protect uE_H/vE_H --- .../vertical/MOM_CVMix_KPP.F90 | 37 +++++++++++-------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 0c1d40e7ec..f0ef355374 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1093,10 +1093,13 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl do k=1,GV%ke U_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)) V_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)) - uE_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)-Waves%US_x(I,j,k)-Waves%US_x(I-1,j,k)) - vE_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)-Waves%US_y(i,J,k)-Waves%US_y(i,J-1,k)) enddo - + if (CS%StokesMOST) then + do k=1,GV%ke + uE_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)-Waves%US_x(I,j,k)-Waves%US_x(I-1,j,k)) + vE_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)-Waves%US_y(i,J,k)-Waves%US_y(i,J-1,k)) + enddo + endif ! things independent of position within the column Coriolis = 0.25*US%s_to_T*( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) + & (G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) ) @@ -1134,7 +1137,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl endif enddo - if (CS%StokesMOST) then + if (CS%StokesMOST) then surfBuoyFlux = buoy_scale * & (buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,k)+buoyFlux(i,j,k+1)) ) surfBuoyFlux2(k) = surfBuoyFlux @@ -1168,14 +1171,12 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl surfSalt = surfHsalt * I_hTot surfU = surfHu * I_hTot surfV = surfHv * I_hTot + Uk = uE_H(k) + uS_H(k) - surfU Vk = vE_H(k) + vS_H(k) - surfV else !not StokesMOST - uS_H(k) = 0.5*(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) - vS_H(k) = 0.5*(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) StokesXI_1d(k) = 0.0 - ! average temperature, salinity, u and v over surface layer ! use C-grid average to get u and v on T-points. surfHtemp = 0.0 @@ -1204,14 +1205,20 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl endif enddo - I_hTot = 1./hTot - surfTemp = surfHtemp * I_hTot - surfSalt = surfHsalt * I_hTot - surfU = surfHu * I_hTot - surfV = surfHv * I_hTot - surfUs = surfHus * I_hTot - surfVs = surfHvs * I_hTot - + !I_hTot = 1./hTot + !surfTemp = surfHtemp * I_hTot + !surfSalt = surfHsalt * I_hTot + !surfU = surfHu * I_hTot + !surfV = surfHv * I_hTot + !surfUs = surfHus * I_hTot + !surfVs = surfHvs * I_hTot + + surfTemp = surfHtemp / hTot + surfSalt = surfHsalt / hTot + surfU = surfHu / hTot + surfV = surfHv / hTot + surfUs = surfHus / hTot + surfVs = surfHvs / hTot ! vertical shear between present layer and surface layer averaged surfU and surfV. ! C-grid average to get Uk and Vk on T-points. Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU