From 61734e6ce9e32af5e327393658a3166e4cdb07c6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 Jul 2025 10:58:05 -0600 Subject: [PATCH 1/9] Update stokes_most --- src/shared/cvmix_kpp.F90 | 289 ++++++++++++++++++++++++++++++--------- 1 file changed, 228 insertions(+), 61 deletions(-) diff --git a/src/shared/cvmix_kpp.F90 b/src/shared/cvmix_kpp.F90 index 7b2b9ad5..3d770dd0 100644 --- a/src/shared/cvmix_kpp.F90 +++ b/src/shared/cvmix_kpp.F90 @@ -79,7 +79,7 @@ module cvmix_kpp public :: cvmix_kpp_compute_bulk_Richardson public :: cvmix_kpp_compute_turbulent_scales public :: cvmix_kpp_compute_unresolved_shear - public :: cvmix_kpp_composite_Gshape !STOKES_MOST + public :: cvmix_kpp_composite_Gshape !STOKES_MOST public :: cvmix_kpp_compute_StokesXi !STOKES_MOST ! These are public for testing, may end up private later public :: cvmix_kpp_compute_shape_function_coeffs @@ -88,7 +88,8 @@ module cvmix_kpp public :: cvmix_kpp_compute_nu_at_OBL_depth_LMD94 public :: cvmix_kpp_EFactor_model public :: cvmix_kpp_ustokes_SL_model - + public :: cvmix_kpp_compute_ER_depth + public :: cvmix_kpp_quad_fit interface cvmix_coeffs_kpp module procedure cvmix_coeffs_kpp_low @@ -3354,7 +3355,7 @@ subroutine cvmix_kpp_composite_Gshape(sigma , Gat1, Gsig, dGdsig) G_1 = MAX( cvmix_zero , MIN( Gat1 , G_m ) ) - if (sigma .lT. sig_m) then + if (sigma .lt. sig_m) then sig = MAX( sigma , cvmix_zero) Gsig = sig + sig * sig * (a2Gsig + a3Gsig * sig) dGdsig = cvmix_one + sig * (2.0_cvmix_r8 * a2Gsig + 3.0_cvmix_r8 * a3Gsig * sig) @@ -3374,66 +3375,60 @@ end subroutine cvmix_kpp_composite_Gshape ! !IROUTINE: cvmix_coeffs_bkgnd_wrap ! !INTERFACE: - subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, & - surf_buoy_force, surf_fric_vel, omega_w2x, uE, vE, uS, vS, & - uSbar, vSbar, uS_SLD, vS_SLD, uSbar_SLD, vSbar_SLD, & - StokesXI, CVmix_kpp_params_user) - ! !DESCRIPTION: -! Compute the Stokes similarity parameter, StokesXI, and Entrainment Rule, BEdE\_ER, from -! surface layer integrated TKE production terms as parameterized in -! Large et al., 2020 (doi:10.1175/JPO-D-20-0308.1) +! Compute the Stokes similarity parameter StokesXI, and Entrainment Rule BEdE, +! and SL integrated tke production terms as +! parameterized in Large et al., 2021 (doi:10.1175/JPO-D-20-0308.1) !\\ !\\ + subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, surf_buoy_force, & + surfBuoy_NS,surf_fric_vel, omega_w2x, uE, vE, uS, vS, uSbar, vSbar, & + uS_SL, vS_SL, uSb_SL, vSb_SL, & + StokesXI,BEdE_ER,PU_TKE,PS_TKE,PB_TKE,CVmix_kpp_params_user) + ! !INPUT PARAMETERS: - real(cvmix_r8), dimension(:), intent(in) :: zi, zk !< Cell interface and center heights < 0 + real(cvmix_r8), dimension(:), intent(in) :: zi, zk !< Cell interface and center heights <= 0 integer, intent(in) :: kSL !< cell index of Surface Layer Depth - real(cvmix_r8), intent(in) :: SLDepth !< Surface Layer Depth Integration limit - real(cvmix_r8), intent(in) :: surf_buoy_force !< Surface buoyancy flux forcing + real(cvmix_r8), intent(in) :: SLDepth !< Surface Layer Depth > 0 + real(cvmix_r8), intent(in) :: surf_buoy_force,surfBuoy_NS !< Surface buoyancy flux forcing, non-solar real(cvmix_r8), intent(in) :: surf_fric_vel, omega_w2x !< Surface wind forcing from x-axis - real(cvmix_r8), dimension(:), intent(in) :: uE, vE !< Eulerian velocity at centers - real(cvmix_r8), intent(in) :: uS_SLD, vS_SLD !< Stokes drift at SLDepth - real(cvmix_r8), intent(in) :: uSbar_SLD, vSbar_SLD !< Average Stokes drift cell kSL to SLDepth - + real(cvmix_r8), dimension(:), intent(in) :: uE, vE !< Eulerian velocity at cell centers + real(cvmix_r8), dimension(:), intent(in) :: uS, vS !< Stokes drift at interfaces + real(cvmix_r8), dimension(:), intent(in) :: uSbar, vSbar !< Cell average Stokes drift + real(cvmix_r8), intent(in) :: uS_SL, vS_SL !< Stokes drift at SLDepth + real(cvmix_r8), intent(in) :: uSb_SL, vSb_SL !< Average Stokes drift cell top to SLDepth type(cvmix_kpp_params_type), intent(in), optional, target :: CVmix_kpp_params_user -! !INPUT/OUTPUT PARAMETERS: - real(cvmix_r8), dimension(:), intent(inout) :: uS, vS !< Stokes drift at interfaces - real(cvmix_r8), dimension(:), intent(inout) :: uSbar, vSbar !< Cell average Stokes drift +! !OUTPUT PARAMETERS: real(cvmix_r8), intent(inout) :: StokesXI !< Stokes similarity parameter + real(cvmix_r8), intent(inout) :: BEdE_ER ! Entrainment rule product [m3 s-3] + real(cvmix_r8), intent(inout) :: PU_TKE,PS_TKE,PB_TKE ! Shear, Stkes, Buoyancy SL TKE Production [m3 s-3] !EOP !BOC +! !LOCAL PARAMETERS: type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in - real(cvmix_r8) :: PU, PS , PB ! surface layer TKE production terms - real(cvmix_r8) :: uS_TMP, vS_TMP, uSbar_TMP, vSbar_TMP ! Temporary Store + real(cvmix_r8) :: PU, PS , PB ! surface layer TKE production terms real(cvmix_r8) :: ustar, delH, delU, delV, omega_E2x, cosOmega, sinOmega real(cvmix_r8) :: BLDepth, TauMAG, TauCG, TauDG, taux0, tauy0, Stk0 , Pinc - real(cvmix_r8) :: PBfact , CempCGm ! Empirical constants - real(cvmix_r8) :: dtop, tauEtop, tauxtop, tauytop ! Cell top values - real(cvmix_r8) :: dbot, tauEbot, tauxbot, tauybot, sigbot, Gbot ! Cell bottom values - integer :: ktmp ! vertical loop + real(cvmix_r8) :: PBfact , CempCGm, CempCGs ! Empirical constants + real(cvmix_r8) :: Cu, Cs, Cb ! Entrainment Rule coefficients + real(cvmix_r8) :: dtop, tauEtop, tauxtop, tauytop ! Cell top values + real(cvmix_r8) :: dbot, tauEbot, tauxbot, tauybot, sigbot, Gbot ! Cell bottom values + integer :: ktmp ! vertical loop CVmix_kpp_params_in => CVmix_kpp_params_saved if (present(CVmix_kpp_params_user)) then CVmix_kpp_params_in => CVmix_kpp_params_user end if - if ( CVmix_kpp_params_in%lStokesMOST ) then - - ! Move bottom of cell kSL up to Surface Layer Extent = SLDepth - uS_TMP = uS(kSL+1) - vS_TMP = vS(kSL+1) - uSbar_TMP = uSbar(kSL) - vSbar_TMP = vSbar(kSL) - uS(kSL+1) = uS_SLD - vS(kSL+1) = vS_SLD - uSbar(kSL)= uSbar_SLD - vSbar(kSL)= vSbar_SLD - - CempCGm= 3.5_cvmix_r8 + Cu = 0.023_cvmix_r8 + Cs = 0.038_cvmix_r8 + Cb = 0.96_cvmix_r8 + CempCGm = 3.5_cvmix_r8 + CempCGs = 4.7_cvmix_r8 ustar = MAX( surf_fric_vel , 1.e-4_cvmix_r8 ) ! > 0 taux0 = ustar**2 * cos(omega_w2x) @@ -3442,10 +3437,13 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, & BLDepth = SLDepth / CVmix_kpp_params_in%surf_layer_ext ! Parameterized Buoyancy production of TKE - PBfact = 0.110_cvmix_r8 - PB = PBfact * MAX( -surf_buoy_force * BLdepth , cvmix_zero ) + PBfact = 0.0893759_cvmix_r8 + PB = MAX( -surfBuoy_NS * BLdepth * PBfact , cvmix_zero ) + PBfact = 0.00215_cvmix_r8 * CempCGs + PB = PB + PBfact * BLdepth * ( abs(surf_buoy_force) - surf_buoy_force ) - ! Compute Both Shear Production Terms down from Surface = initial top values + cosOmega = 0.0 + sinOmega = 0.0 PU = 0.0 PS = 0.0 dtop = 0.0 @@ -3455,18 +3453,19 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, & tauxtop = taux0 tauytop = tauy0 - do ktmp = 1, kSL - ! SLdepth can be between cell interfaces kSL and kSL+1 - delH = min( max(cvmix_zero, SLdepth - dtop), (zi(ktmp) - zi(ktmp+1) ) ) - dbot = MIN( dtop + delH , SLdepth) - sigbot = dbot / BLdepth - Gbot = cvmix_kpp_composite_shape(sigbot) - TauMAG = ustar * ustar * Gbot / sigbot +! Integrate Both Shear Production Terms from Surface to top of surface layer cell + do ktmp = 1, kSL-1 delU = uE(ktmp) - uE(ktmp+1) delV = vE(ktmp) - vE(ktmp+1) Omega_E2x= atan2( delV , delU ) cosOmega = cos(Omega_E2x) sinOmega = sin(Omega_E2x) + + delH = zi(ktmp) - zi(ktmp+1) + dbot = dtop + delH + sigbot = dbot / BLdepth + Gbot = cvmix_kpp_composite_shape(sigbot) + TauMAG = ustar * ustar * Gbot / sigbot tauCG = CempCGm * Gbot * (taux0 * cosOmega - tauy0 * sinOmega) ! tauDG = sqrt( TauMAG**2 - tauCG**2 ) ! G tauDG = TauMAG ! E @@ -3490,22 +3489,190 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, & tauEtop = tauEbot enddo - ! Compute Stokes similarity parameter - StokesXI = PS / MAX( PU + PS + PB , 1.e-12_cvmix_r8 ) +! Integrate from top of surface layer cell to Surface layer Depth + delH = SLDepth + zi(ktmp) + sigbot = CVmix_kpp_params_in%surf_layer_ext + Gbot = cvmix_kpp_composite_shape(sigbot) + TauMAG = ustar * ustar * Gbot / sigbot + tauCG = CempCGm * Gbot * (taux0 * cosOmega - tauy0 * sinOmega) + ! tauDG = sqrt( TauMAG**2 - tauCG**2 ) ! G + tauDG = TauMAG ! E + tauxbot = tauDG * cosOmega - tauCG * sinOmega + tauybot = tauDG * sinOmega + tauCG * cosOmega + tauEbot = (tauxbot * delU + tauybot * delV) / (zk(ktmp) - zk(ktmp+1) ) + ! Increment Eulerian Shear Production + Pinc = 0.5_cvmix_r8 * (tauEbot + tauEtop) * delH + PU = PU + MAX( Pinc , cvmix_zero ) + + ! Increment Stokes Shear Production + Pinc = tauxtop*uS(ktmp) - tauxbot*uS_SL + tauytop*vS(ktmp) - tauybot*vS_SL + Pinc = Pinc - (tauxtop-tauxbot) * uSb_SL - (tauytop-tauybot) * vSb_SL + PS = PS + MAX( Pinc , cvmix_zero ) + + ! Compute Stokes similarity parameter, Entrainment Rule, and TKE prodction Rates + if ( (PU + PB + PS) > cvmix_zero ) then + StokesXI = PS / (PU + PS + PB) + else + StokesXI = cvmix_zero + endif + BEdE_ER = MAX( ( Cu*PU + Cs*PS + Cb*PB ) , cvmix_zero ) + PU_TKE = PU + PS_TKE = PS + PB_TKE = PB + +!EOC + + end subroutine cvmix_kpp_compute_StokesXi + +!BOP + +! !DESCRIPTION: +! Use Entrainment Rule, BEdE_ER, To Find Entrainment Flux and Depth +! for each assumed OBL_depth = cell centers, +! until the boundary layer depth, ERdepth > 0 kER_depth are determined, OR +! if there is no viable solution ERdepth = -1 , kER_depth=-1 + + subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,deltaRho,OBL_depth, & + uStar,Bsfc_ns,surfBuoy,StokesXI,BEdE_ER,ERdepth, & + CVMix_kpp_params_user) + +! !INPUT PARAMETERS: + real(cvmix_r8), dimension(:), intent(in) :: & + z_inter, & ! Interface heights <= 0 [m] + Nsq ! Column of Buoyancy Gradients at interfaces + real(cvmix_r8), dimension(:), intent(in) :: & + OBL_depth, & ! Array of assumed OBL depths >0 at cell centers [m] + deltaRho, & ! Column of Local density difference from surface at centers [kg m-3] + surfBuoy, & ! surface Buoyancy flux surface to OBL_depth + StokesXI, & ! Stokes similarity parameter given OBL_depth + BEdE_ER ! Parameterized Entrainment Rule given OBL_depth + real(cvmix_r8), intent(in) :: uStar ! surface friction velocity [m s-1] + real(cvmix_r8), intent(in) :: Bsfc_ns ! non-solar surface buoyancy flux boundary condition + + type(cvmix_kpp_params_type), optional, target, intent(in) :: CVmix_kpp_params_user + +! !OUTPUT PARAMETERS: + real(cvmix_r8), intent(out) :: ERdepth ! ER Boundary Layer Depth [m] +!EOP +!BOC - ! Restore bottom of cell kSL at zi(kSL+1) with stored Stokes Drift ; ditto average over cell kSL - uS(kSL+1) = uS_TMP - vS(kSL+1) = vS_TMP - uSbar(kSL) = uSbar_TMP - vSbar(kSL) = vSbar_TMP + ! Local variables + integer :: iz, nlev , kbl , kinv + real(cvmix_r8), dimension(size(OBL_depth)+1) :: zdepth, BEdE, BEnt ! surface then cell-centers<0 + real(cvmix_r8), dimension(size(z_inter)+1) :: sigma, Bflux ! interface values + real(cvmix_r8) :: ws_i, Cemp_Rs, Gsig_i, Fdelrho, Bnonlocal, sigE, maxNsq + real(kind=cvmix_r8), dimension(4) :: coeffs + type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in - else ! not lStokesMOST - StokesXI = cvmix_zero - end if + CVmix_kpp_params_in => CVmix_kpp_params_saved + if (present(CVmix_kpp_params_user)) then + CVmix_kpp_params_in => CVmix_kpp_params_user + end if + + nlev = size(OBL_depth) + Cemp_Rs = 4.7_cvmix_r8 +! Gat1 = cvmix_zero + Fdelrho = cvmix_one +! kinv = MAXLOC( Nsq ) ! interface index of maximum stratification, N2>0 (inversion) +! maxNsq = Nsq( kinv ) + maxNsq = 0.0 + do kbl = 2, nlev+1 + if ( Nsq(kbl) > maxNsq ) then + kinv = kbl + maxNsq = Nsq(kbl) + endif + enddo + +! Set default output values (no solution) + ERdepth = -cvmix_one +! Set surface values + zdepth(1) = cvmix_zero + BEdE(1) = cvmix_zero + BEnt(1) = cvmix_zero + sigma(:) = cvmix_zero + Bflux(1) = Bsfc_ns ! non-solar surface buoyancy boundary condition for all kbl +! Set OBL_depth(1)=top cell center values + BEnt(1) = cvmix_zero + zdepth(2) = -OBL_depth(1) + BEdE(2) = cvmix_zero + + do kbl = 2, nlev + zdepth(kbl+1)= -OBL_depth(kbl) + BEdE(kbl+1) = cvmix_zero + + sigma(kbl+1) = cvmix_one + Bflux(kbl+1) = cvmix_zero + sigma(kbl+2) = -z_inter(kbl+1)/OBL_depth(kbl) ! > 1.0 + Bflux(kbl+2) = cvmix_zero + Bnonlocal = 0.5_cvmix_r8 * Cemp_Rs * (abs(surfBuoy(kbl)) - surfBuoy(kbl)) + + do iz = kbl,1,-1 + if (iz .gt. 1) then + sigma(iz) = -z_inter(iz)/OBL_depth(kbl) ! < 1.0 + call cvmix_kpp_compute_turbulent_scales(sigma(iz), OBL_depth(kbl), surfBuoy(kbl), uStar, StokesXI(kbl), & !0d + w_s=ws_i , CVmix_kpp_params_user=CVmix_kpp_params_user) + Gsig_i = cvmix_kpp_composite_shape(sigma(iz)) + Bflux(iz) = Gsig_i * (OBL_depth(kbl) * ws_i * Nsq(iz) - Bnonlocal) + endif + + ! find the peak + if ( (Bflux(iz+1) > Bflux(iz+2)) .and. (Bflux(iz+1) .ge. Bflux(iz)) .and. (Bflux(iz+1) > cvmix_zero) ) then + call cvmix_kpp_quad_fit(iz, sigma, Bflux, sigE, BEnt(kbl+1) ) + Fdelrho = cvmix_one + BEdE(kbl+1) = Fdelrho*BEnt(kbl+1)*sigE*OBL_depth(kbl) + exit ! No exit leaves BEdE(kbl+1) = cvmix_zero + endif + enddo + + if ( (BEdE(kbl+1) > BEdE_ER(kbl)) .and. (Bsfc_ns < cvmix_zero) ) then + call cvmix_math_poly_interp(coeffs, CVmix_kpp_params_in%interp_type, & + zdepth(kbl:kbl+1) , BEdE(kbl:kbl+1), & + zdepth(kbl-1) , BEdE(kbl-1) ) ! surface values if kbl=2; + coeffs(1) = coeffs(1) - BEdE_ER(kbl) + ERdepth = -cvmix_math_cubic_root_find(coeffs, 0.5_cvmix_r8 * & + (zdepth(kbl)+zdepth(kbl+1))) + + exit + endif + + enddo !EOC - end subroutine cvmix_kpp_compute_StokesXi + end subroutine cvmix_kpp_compute_ER_depth + + +! !DESCRIPTION: +! Find maximum y0, at x0; given values of x and y at i, i+1 and i+2 + + subroutine cvmix_kpp_quad_fit( i, x, y, x0, y0 ) + +! !INPUT PARAMETERS: + real(cvmix_r8), dimension(:), intent(in) :: x, y !< + integer, intent(in) :: i !< + +! !OUTPUT PARAMETERS: + real(cvmix_r8), intent(out) :: x0, y0 !< quadratic fit extreme + +! !LOCAL PARAMETERS: + real(cvmix_r8) :: a , b, c !< y = a x^2 + b x + c + + a = (x(i+1)-x(i)) * (x(i+2)*x(i+2) - x(i+1)*x(i+1)) - (x(i+2)-x(i+1)) * (x(i+1)*x(i+1) - x(i)*x(i)) + if (a == 0.0) then + x0 = x(i+1) + y0 = y(i+1) + else + a = ( (x(i+1)-x(i)) * (y(i+2)-y(i+1)) - (x(i+2)-x(i+1)) * (y(i+1)-y(i)) ) / a + + b = ( (y(i+1) - y(i)) - a * ( x(i+1)*x(i+1) - x(i)*x(i) ) ) / (x(i+1)-x(i)) + + c = y(i) - a * x(i)*x(i) - b * x(i) + + x0 = real(-0.5,cvmix_r8) * b / a + y0 = a * x0*x0 + b * x0 + c + endif + + end subroutine cvmix_kpp_quad_fit end module cvmix_kpp From 268829e943df11b5e5144dd5e60b07361aa0741b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 Jul 2025 11:10:44 -0600 Subject: [PATCH 2/9] Delete unsed argument, deltaRho --- src/shared/cvmix_kpp.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/shared/cvmix_kpp.F90 b/src/shared/cvmix_kpp.F90 index 3d770dd0..e5762e0c 100644 --- a/src/shared/cvmix_kpp.F90 +++ b/src/shared/cvmix_kpp.F90 @@ -3532,7 +3532,7 @@ end subroutine cvmix_kpp_compute_StokesXi ! until the boundary layer depth, ERdepth > 0 kER_depth are determined, OR ! if there is no viable solution ERdepth = -1 , kER_depth=-1 - subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,deltaRho,OBL_depth, & + subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & uStar,Bsfc_ns,surfBuoy,StokesXI,BEdE_ER,ERdepth, & CVMix_kpp_params_user) @@ -3542,7 +3542,6 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,deltaRho,OBL_depth, & Nsq ! Column of Buoyancy Gradients at interfaces real(cvmix_r8), dimension(:), intent(in) :: & OBL_depth, & ! Array of assumed OBL depths >0 at cell centers [m] - deltaRho, & ! Column of Local density difference from surface at centers [kg m-3] surfBuoy, & ! surface Buoyancy flux surface to OBL_depth StokesXI, & ! Stokes similarity parameter given OBL_depth BEdE_ER ! Parameterized Entrainment Rule given OBL_depth From 81c42655c5bc771760c2b787764f453712b21d07 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 Jul 2025 11:15:27 -0600 Subject: [PATCH 3/9] Replace *SLD with *SL --- src/drivers/cvmix_kpp_drv.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/drivers/cvmix_kpp_drv.F90 b/src/drivers/cvmix_kpp_drv.F90 index 37be82bb..1afdafbc 100644 --- a/src/drivers/cvmix_kpp_drv.F90 +++ b/src/drivers/cvmix_kpp_drv.F90 @@ -703,10 +703,10 @@ Subroutine cvmix_kpp_driver() vS=vS, & uSbar=uSbar, & vSbar=vSbar, & - uS_SLD=cvmix_one, & - vS_SLD=cvmix_one, & - uSbar_SLD=cvmix_one, & - vSbar_SLD=cvmix_one, & + uS_SL=cvmix_one, & + vS_SL=cvmix_one, & + uSbar_SL=cvmix_one, & + vSbar_SL=cvmix_one, & StokesXI=StokesXI) call cvmix_kpp_compute_OBL_depth(Ri_bulk, zw_iface, OBL_depth7, & kOBL_depth, zt, surf_fric=cvmix_one, & From 09be4f4e2dce3f742e3e589b4660705c91af5694 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 Jul 2025 11:23:26 -0600 Subject: [PATCH 4/9] Replace uSbar_SL => uSb_SL; vSbar_SL => vSb_SL --- src/drivers/cvmix_kpp_drv.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/drivers/cvmix_kpp_drv.F90 b/src/drivers/cvmix_kpp_drv.F90 index 1afdafbc..1ad8d8ac 100644 --- a/src/drivers/cvmix_kpp_drv.F90 +++ b/src/drivers/cvmix_kpp_drv.F90 @@ -705,8 +705,8 @@ Subroutine cvmix_kpp_driver() vSbar=vSbar, & uS_SL=cvmix_one, & vS_SL=cvmix_one, & - uSbar_SL=cvmix_one, & - vSbar_SL=cvmix_one, & + uSb_SL=cvmix_one, & + vSb_SL=cvmix_one, & StokesXI=StokesXI) call cvmix_kpp_compute_OBL_depth(Ri_bulk, zw_iface, OBL_depth7, & kOBL_depth, zt, surf_fric=cvmix_one, & From bb136050bfa808184968523dfcaa4662473764a7 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 Jul 2025 11:47:55 -0600 Subject: [PATCH 5/9] Fix args passed to cvmix_kpp_driver --- src/drivers/cvmix_kpp_drv.F90 | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/drivers/cvmix_kpp_drv.F90 b/src/drivers/cvmix_kpp_drv.F90 index 1ad8d8ac..bf574616 100644 --- a/src/drivers/cvmix_kpp_drv.F90 +++ b/src/drivers/cvmix_kpp_drv.F90 @@ -62,7 +62,7 @@ Subroutine cvmix_kpp_driver() real(cvmix_r8) :: hmix1, hmix5, hmix7, ri_crit, layer_thick1, & layer_thick4, layer_thick5, layer_thick7, OBL_depth4, & OBL_depth5, OBL_depth7, N, Nsqr - real(cvmix_r8) :: StokesXI + real(cvmix_r8) :: StokesXI, BEdE_ER, PU_TKE, PS_TKE, PB_TKE real(cvmix_r8) :: kOBL_depth, Bslope, Vslope real(cvmix_r8) :: sigma6, OBL_depth6, surf_buoy_force6, surf_fric_vel6, & vonkarman6, tao, rho0, grav, alpha, Qnonpen, Cp0, & @@ -690,11 +690,16 @@ Subroutine cvmix_kpp_driver() vSbar(nlev7), & source=cvmix_one) StokesXI = cvmix_one + BEdE_ER = cvmix_one + PU_TKE = cvmix_one + PS_TKE = cvmix_one + PB_TKE = cvmix_one call cvmix_kpp_compute_StokesXi(zi=zt, & zk=zw_iface, & kSL=(nlev7+1)/2, & SLDepth=cvmix_zero, & surf_buoy_force=cvmix_one, & + surfBuoy_NS=cvmix_one, & surf_fric_vel=cvmix_one, & omega_w2x=cvmix_one, & uE=ones, & @@ -703,11 +708,15 @@ Subroutine cvmix_kpp_driver() vS=vS, & uSbar=uSbar, & vSbar=vSbar, & - uS_SL=cvmix_one, & - vS_SL=cvmix_one, & - uSb_SL=cvmix_one, & - vSb_SL=cvmix_one, & - StokesXI=StokesXI) + uS_SL=cvmix_one, & + vS_SL=cvmix_one, & + uSb_SL=cvmix_one, & + vSb_SL=cvmix_one, & + StokesXI=StokesXI, & + BEdE_ER=BEdE_ER, & + PU_TKE=PU_TKE, & + PS_TKE=PS_TKE, & + PB_TKE=PB_TKE) call cvmix_kpp_compute_OBL_depth(Ri_bulk, zw_iface, OBL_depth7, & kOBL_depth, zt, surf_fric=cvmix_one, & surf_buoy=ones, Xi=ones) From 9f4da9c6adada153435a1c38144f837c9d54e71c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 Jul 2025 13:11:37 -0600 Subject: [PATCH 6/9] Set SLDepth=cvmix_one and change BLDepth -> BLdepth --- src/drivers/cvmix_kpp_drv.F90 | 2 +- src/shared/cvmix_kpp.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/drivers/cvmix_kpp_drv.F90 b/src/drivers/cvmix_kpp_drv.F90 index bf574616..57d6df80 100644 --- a/src/drivers/cvmix_kpp_drv.F90 +++ b/src/drivers/cvmix_kpp_drv.F90 @@ -697,7 +697,7 @@ Subroutine cvmix_kpp_driver() call cvmix_kpp_compute_StokesXi(zi=zt, & zk=zw_iface, & kSL=(nlev7+1)/2, & - SLDepth=cvmix_zero, & + SLDepth=cvmix_one, & surf_buoy_force=cvmix_one, & surfBuoy_NS=cvmix_one, & surf_fric_vel=cvmix_one, & diff --git a/src/shared/cvmix_kpp.F90 b/src/shared/cvmix_kpp.F90 index e5762e0c..0e8aaa98 100644 --- a/src/shared/cvmix_kpp.F90 +++ b/src/shared/cvmix_kpp.F90 @@ -3434,7 +3434,7 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, surf_buoy_force, taux0 = ustar**2 * cos(omega_w2x) tauy0 = ustar**2 * sin(omega_w2x) Stk0 = sqrt( uS(1)**2 + vS(1)**2 ) - BLDepth = SLDepth / CVmix_kpp_params_in%surf_layer_ext + BLdepth = SLDepth / CVmix_kpp_params_in%surf_layer_ext ! Parameterized Buoyancy production of TKE PBfact = 0.0893759_cvmix_r8 From 20eee020983a84582fea913c6c634e52a0d336d0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 20 Aug 2025 13:22:58 -0600 Subject: [PATCH 7/9] Addresses Mike Levy comments --- src/shared/cvmix_kpp.F90 | 37 +++++++++++++++---------------------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/src/shared/cvmix_kpp.F90 b/src/shared/cvmix_kpp.F90 index 0e8aaa98..6364eae1 100644 --- a/src/shared/cvmix_kpp.F90 +++ b/src/shared/cvmix_kpp.F90 @@ -3410,11 +3410,15 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, surf_buoy_force, ! !LOCAL PARAMETERS: type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in - real(cvmix_r8) :: PU, PS , PB ! surface layer TKE production terms + real(cvmix_r8), parameter :: CempCGm = 3.5_cvmix_r8 ! Coeff. relating cross-shear mtm flux to sfc stress [nondim] + real(cvmix_r8), parameter :: CempCGs = 4.7_cvmix_r8 ! Coeff. relating non-local scalar flux to sfc flux [nondim] + real(cvmix_r8), parameter :: Cu = 0.023_cvmix_r8 ! TKE shear production weight [nondim] + real(cvmix_r8), parameter :: Cs = 0.038_cvmix_r8 ! TKE Stokes production weight [nondim] + real(cvmix_r8), parameter :: Cb = 0.96_cvmix_r8 ! TKE buoyancy production weight [nondim] + real(cvmix_r8) :: PBfact ! Ratio of TKE surface layer production to w*^3 [nondim] + real(cvmix_r8) :: PU, PS , PB ! Surface layer TKE production terms increments [m3 s-3] real(cvmix_r8) :: ustar, delH, delU, delV, omega_E2x, cosOmega, sinOmega real(cvmix_r8) :: BLDepth, TauMAG, TauCG, TauDG, taux0, tauy0, Stk0 , Pinc - real(cvmix_r8) :: PBfact , CempCGm, CempCGs ! Empirical constants - real(cvmix_r8) :: Cu, Cs, Cb ! Entrainment Rule coefficients real(cvmix_r8) :: dtop, tauEtop, tauxtop, tauytop ! Cell top values real(cvmix_r8) :: dbot, tauEbot, tauxbot, tauybot, sigbot, Gbot ! Cell bottom values integer :: ktmp ! vertical loop @@ -3424,12 +3428,6 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, surf_buoy_force, CVmix_kpp_params_in => CVmix_kpp_params_user end if - Cu = 0.023_cvmix_r8 - Cs = 0.038_cvmix_r8 - Cb = 0.96_cvmix_r8 - CempCGm = 3.5_cvmix_r8 - CempCGs = 4.7_cvmix_r8 - ustar = MAX( surf_fric_vel , 1.e-4_cvmix_r8 ) ! > 0 taux0 = ustar**2 * cos(omega_w2x) tauy0 = ustar**2 * sin(omega_w2x) @@ -3515,10 +3513,10 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, surf_buoy_force, else StokesXI = cvmix_zero endif - BEdE_ER = MAX( ( Cu*PU + Cs*PS + Cb*PB ) , cvmix_zero ) - PU_TKE = PU - PS_TKE = PS - PB_TKE = PB + BEdE_ER = MAX( ( Cu*PU + Cs*PS + Cb*PB ) , cvmix_zero ) + PU_TKE = PU + PS_TKE = PS + PB_TKE = PB !EOC @@ -3558,9 +3556,9 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & ! Local variables integer :: iz, nlev , kbl , kinv - real(cvmix_r8), dimension(size(OBL_depth)+1) :: zdepth, BEdE, BEnt ! surface then cell-centers<0 + real(cvmix_r8), dimension(size(OBL_depth)+1) :: zdepth, BEdE ! surface then cell-centers<0 real(cvmix_r8), dimension(size(z_inter)+1) :: sigma, Bflux ! interface values - real(cvmix_r8) :: ws_i, Cemp_Rs, Gsig_i, Fdelrho, Bnonlocal, sigE, maxNsq + real(cvmix_r8) :: ws_i, Cemp_Rs, Gsig_i, Fdelrho, Bnonlocal, sigE, maxNsq, BEnt real(kind=cvmix_r8), dimension(4) :: coeffs type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in @@ -3571,10 +3569,7 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & nlev = size(OBL_depth) Cemp_Rs = 4.7_cvmix_r8 -! Gat1 = cvmix_zero Fdelrho = cvmix_one -! kinv = MAXLOC( Nsq ) ! interface index of maximum stratification, N2>0 (inversion) -! maxNsq = Nsq( kinv ) maxNsq = 0.0 do kbl = 2, nlev+1 if ( Nsq(kbl) > maxNsq ) then @@ -3588,11 +3583,9 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & ! Set surface values zdepth(1) = cvmix_zero BEdE(1) = cvmix_zero - BEnt(1) = cvmix_zero sigma(:) = cvmix_zero Bflux(1) = Bsfc_ns ! non-solar surface buoyancy boundary condition for all kbl ! Set OBL_depth(1)=top cell center values - BEnt(1) = cvmix_zero zdepth(2) = -OBL_depth(1) BEdE(2) = cvmix_zero @@ -3617,9 +3610,9 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & ! find the peak if ( (Bflux(iz+1) > Bflux(iz+2)) .and. (Bflux(iz+1) .ge. Bflux(iz)) .and. (Bflux(iz+1) > cvmix_zero) ) then - call cvmix_kpp_quad_fit(iz, sigma, Bflux, sigE, BEnt(kbl+1) ) + call cvmix_kpp_quad_fit(iz, sigma, Bflux, sigE, BEnt ) Fdelrho = cvmix_one - BEdE(kbl+1) = Fdelrho*BEnt(kbl+1)*sigE*OBL_depth(kbl) + BEdE(kbl+1) = Fdelrho*BEnt*sigE*OBL_depth(kbl) exit ! No exit leaves BEdE(kbl+1) = cvmix_zero endif enddo From 66a067cac14fa1dd50f8a1c0355e15768b9cd691 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Wed, 27 Aug 2025 14:17:04 -0600 Subject: [PATCH 8/9] Remove cvmix_kpp_quad_fit The function cvmix_kpp_quad_fit did two things: 1. Given three values of (sigma, Bflux), find quadratic polynomial that interpolates function 2. Find the root of the derivative of the quadratic, and evaluate the polynomial at that point I created an interface for cvmix_math_poly_interp() to do the interpolation (that routine could do quadratic interpolation, but based on two points and a derivative instead of three points), and then find the derivative and evaluate the polynomial in cvmix_kpp_compute_ER_depth() instead of a separate function --- src/shared/cvmix_kpp.F90 | 51 +++++++------------- src/shared/cvmix_math.F90 | 99 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 112 insertions(+), 38 deletions(-) diff --git a/src/shared/cvmix_kpp.F90 b/src/shared/cvmix_kpp.F90 index 6364eae1..ff3ced08 100644 --- a/src/shared/cvmix_kpp.F90 +++ b/src/shared/cvmix_kpp.F90 @@ -89,7 +89,6 @@ module cvmix_kpp public :: cvmix_kpp_EFactor_model public :: cvmix_kpp_ustokes_SL_model public :: cvmix_kpp_compute_ER_depth - public :: cvmix_kpp_quad_fit interface cvmix_coeffs_kpp module procedure cvmix_coeffs_kpp_low @@ -3609,8 +3608,22 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & endif ! find the peak - if ( (Bflux(iz+1) > Bflux(iz+2)) .and. (Bflux(iz+1) .ge. Bflux(iz)) .and. (Bflux(iz+1) > cvmix_zero) ) then - call cvmix_kpp_quad_fit(iz, sigma, Bflux, sigE, BEnt ) + if ( (Bflux(iz+1) .gt. Bflux(iz+2)) .and. (Bflux(iz+1) .ge. Bflux(iz)) .and. & + (Bflux(iz+1).gt. cvmix_zero) ) then + ! Find sigE (the root of the derivative of the quadratic polynomial + ! interpolating (sigma(i), Bflux(i)) for i in [iz, iz+1, iz+2]) + ! Also find BEnt (value of quadratic at sigE) + call cvmix_math_poly_interp(coeffs, sigma(iz:iz+2), Bflux(iz:iz+2)) + ! coeffs(3) = 0 => sigma(iz), sigma(iz+1), and sigma(iz+2) are not unique values + ! so there interpolation returned a linear equation. In this case we select + ! (sigma(iz+1), Bflux(iz+1)) as the maximum. + if (coeffs(3) .eq. cvmix_zero) then + sigE = sigma(iz+1) + Bent = Bflux(iz+1) + else + sigE = -0.5_cvmix_r8 * (coeffs(2) / coeffs(3)) + Bent = cvmix_math_evaluate_cubic(coeffs, sigE) + endif Fdelrho = cvmix_one BEdE(kbl+1) = Fdelrho*BEnt*sigE*OBL_depth(kbl) exit ! No exit leaves BEdE(kbl+1) = cvmix_zero @@ -3635,36 +3648,4 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & end subroutine cvmix_kpp_compute_ER_depth -! !DESCRIPTION: -! Find maximum y0, at x0; given values of x and y at i, i+1 and i+2 - - subroutine cvmix_kpp_quad_fit( i, x, y, x0, y0 ) - -! !INPUT PARAMETERS: - real(cvmix_r8), dimension(:), intent(in) :: x, y !< - integer, intent(in) :: i !< - -! !OUTPUT PARAMETERS: - real(cvmix_r8), intent(out) :: x0, y0 !< quadratic fit extreme - -! !LOCAL PARAMETERS: - real(cvmix_r8) :: a , b, c !< y = a x^2 + b x + c - - a = (x(i+1)-x(i)) * (x(i+2)*x(i+2) - x(i+1)*x(i+1)) - (x(i+2)-x(i+1)) * (x(i+1)*x(i+1) - x(i)*x(i)) - if (a == 0.0) then - x0 = x(i+1) - y0 = y(i+1) - else - a = ( (x(i+1)-x(i)) * (y(i+2)-y(i+1)) - (x(i+2)-x(i+1)) * (y(i+1)-y(i)) ) / a - - b = ( (y(i+1) - y(i)) - a * ( x(i+1)*x(i+1) - x(i)*x(i) ) ) / (x(i+1)-x(i)) - - c = y(i) - a * x(i)*x(i) - b * x(i) - - x0 = real(-0.5,cvmix_r8) * b / a - y0 = a * x0*x0 + b * x0 + c - endif - - end subroutine cvmix_kpp_quad_fit - end module cvmix_kpp diff --git a/src/shared/cvmix_math.F90 b/src/shared/cvmix_math.F90 index f9c0df7b..dee7585d 100644 --- a/src/shared/cvmix_math.F90 +++ b/src/shared/cvmix_math.F90 @@ -23,6 +23,7 @@ module cvmix_math ! !USES: use cvmix_kinds_and_types, only : cvmix_r8, & + cvmix_zero, & cvmix_one !EOP @@ -47,16 +48,21 @@ module cvmix_math public :: cvmix_math_cubic_root_find public :: cvmix_math_evaluate_cubic + interface cvmix_math_poly_interp + module procedure linear_through_cubic_poly_interp + module procedure traditional_quad_poly_interp + end interface cvmix_math_poly_interp + !EOP contains !BOP -! !IROUTINE: cvmix_math_poly_interp +! !IROUTINE: linear_through_cubic_poly_interp ! !INTERFACE: - subroutine cvmix_math_poly_interp(coeffs, interp_type, x, y, x0, y0) + subroutine linear_through_cubic_poly_interp(coeffs, interp_type, x, y, x0, y0) ! !DESCRIPTION: ! Given (x(1), y(1)), (x(2), y(2)), and possibly (x0, y0), compute coeffs = @@ -178,14 +184,97 @@ subroutine cvmix_math_poly_interp(coeffs, interp_type, x, y, x0, y0) !EOC - end subroutine cvmix_math_poly_interp + end subroutine linear_through_cubic_poly_interp + +!BOP + +! !IROUTINE: traditional_quad_poly_interp +! !INTERFACE: + + subroutine traditional_quad_poly_interp(coeffs, x, y) + +! !DESCRIPTION: +! Given (x(1), y(1)), (x(2), y(2)), and (x(3), y(3)), compute coeffs = +! $(/a_0, a_1, a_2, a_3=0/)$ such that the quadratic polynomial $f(x) = \sum a_nx^n$ +! interpolates those three points. I.e. $f(x(1)) = y(1)$, $f(x(2)) = y(2)$, and +! $f(x(3)) = y(3)$. +! \\ +! \\ + +! !INPUT PARAMETERS: + real(cvmix_r8), dimension(3), intent(in) :: x, y +! !OUTPUT PARAMETERS: + real(cvmix_r8), dimension(4), intent(inout) :: coeffs + +!EOP +!BOC + + ! Local variables + real(cvmix_r8) :: det, inv_det + + coeffs(:) = cvmix_zero + det = -(x(3) - x(2)) * (x(3) - x(1)) * (x(2) - x(1)) + ! Return all 0s if any two points are the same + if (det .eq. cvmix_zero) then + return + endif + ! Given a matrix + ! + ! [x(1)*x(1) x(1) 1] + ! A = [x(2)*x(2) x(2) 1] + ! [x(3)*x(3) x(3) 1] + ! + ! Want coeff(1), coeff(2), and coeff(3) such that + ! + ! [coeff(3)] [y(1)] + ! A * [coeff(2)] = [y(2)] + ! [coeff(1)] [y(3)] + ! + ! The inverse is + ! 1 [ (x(2)-x(3)) -(x(1)-x(3)) (x(1)-x(2)) ] + ! Ainv = --- * [-(x(2)*x(2)-x(3)*x(3)) (x(1)*x(1)-x(3)*x(3)) -(x(1)*x(1)-x(2)*x(2))] + ! det [ x(2)*x(3)*(x(2)-x(3)) -x(1)*x(3)*(x(1)-x(3)) x(1)*x(2)*(x(1)-x(2))] + ! + ! And the coefficients of interpolating polynomial is the above is + ! + ! [coeff(3)] [y(1)] + ! [coeff(2)] = Ainv * [y(2)] + ! [coeff(1)] [y(3)] + inv_det = cvmix_one / det + coeffs(3) = inv_det * ((x(2)-x(3))*y(1) - (x(1)-x(3))*y(2) + (x(1)-x(2))*y(3)) + coeffs(2) = inv_det * (-(x(2)*x(2)-x(3)*x(3))*y(1) + (x(1)*x(1)-x(3)*x(3))*y(2) -(x(1)*x(1)-x(2)*x(2))*y(3)) + ! Now that we know coeffs(3) and coeffs(2), we can determine coeffs(1) by evaluation: + ! coeffs(1) + coeffs(2)*x(1) + coeffs(3)*x(1)*x(1) = y(1) + ! coeffs(1) = y(1) - coeffs(2)*x(1) - coeffs(3)*x(1)*x(1) + coeffs(1) = y(1) - coeffs(2)*x(1) - coeffs(3)*x(1)*x(1) + +!EOC + + end subroutine traditional_quad_poly_interp + +!BOP + +! !IROUTINE: cvmix_math_cubic_root_find +! !INTERFACE: function cvmix_math_cubic_root_find(coeffs, x0) +! !DESCRIPTION: +! Use Newton's Method to find a root of $f(x) = a_0 + a_1x + a_2x^2 + a_3x^3$ +!\\ +!\\ + +! !INPUT PARAMETERS: real(cvmix_r8), dimension(4), intent(in) :: coeffs real(cvmix_r8), intent(in) :: x0 +! !OUTPUT PARAMETERS: real(cvmix_r8) :: cvmix_math_cubic_root_find + +!EOP +!BOC + + ! Local Variables real(cvmix_r8) :: fun_val, root, slope integer :: it_cnt @@ -200,6 +289,8 @@ function cvmix_math_cubic_root_find(coeffs, x0) end do cvmix_math_cubic_root_find = root +!EOC + end function cvmix_math_cubic_root_find !BOP @@ -244,6 +335,8 @@ function cvmix_math_evaluate_cubic(coeffs, x_in, fprime) fprime = fprime + coeffs(i)*real(i-1,cvmix_r8)*(x_in**(i-2)) end do +!EOC + end function cvmix_math_evaluate_cubic end module cvmix_math From 971b3663685211798ad1a962e027bc79c68a3323 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Wed, 27 Aug 2025 14:26:46 -0600 Subject: [PATCH 9/9] Fix whitespace --- src/shared/cvmix_kpp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shared/cvmix_kpp.F90 b/src/shared/cvmix_kpp.F90 index ff3ced08..82570b02 100644 --- a/src/shared/cvmix_kpp.F90 +++ b/src/shared/cvmix_kpp.F90 @@ -3609,7 +3609,7 @@ subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, & ! find the peak if ( (Bflux(iz+1) .gt. Bflux(iz+2)) .and. (Bflux(iz+1) .ge. Bflux(iz)) .and. & - (Bflux(iz+1).gt. cvmix_zero) ) then + (Bflux(iz+1) .gt. cvmix_zero) ) then ! Find sigE (the root of the derivative of the quadratic polynomial ! interpolating (sigma(i), Bflux(i)) for i in [iz, iz+1, iz+2]) ! Also find BEnt (value of quadratic at sigE)