Update stokes most #106

gustavo-marques · 2025-07-18T17:49:25Z

Add new subroutine cvmix_kpp_compute_ER_depth to calculate the entrainment flux and corresponding depth using the entrainment rule. For each assumed boundary layer depth (OBL_depth) at cell centers, the subroutine iteratively searches for the first depth where the entrainment depth (ERdepth) becomes positive. The corresponding index (kER_depth) is also returned. If no viable solution is found, the subroutine returns ERdepth = -1 and kER_depth = -1.
Add new subroutine cvmix_kpp_quad_fit to estimate the maximum of a quadratic fit.
Update cvmix_kpp_compute_StokesXi
Compute Shear, Stokes, Buoyancy, surface-layer TKE Production

This commit introduces an alternative method for computing the boundary layer depth (BLD) in KPP, based on the entrainment rule. This method is activated when the non-solar surface buoyancy flux (surfBuoy_NS) is negative and STOKES_MOST = True. In all other cases, the standard Richardson Number-based method is used. The non-solar surface buoyancy flux is estimated internally using an exponential attenuation function, since only the total surface buoyancy flux (solar + non-solar) is currently available to the module. Ideally, the non-solar component should be provided explicitly. Several new diagnostics have been added, and some of the existing ones have been updated. Specifically, the parameterized shear, buoyancy, and Stokes TKE production terms are calculated. We conducted multiple tests where these terms are passed to the MLE module and used in the denominator of the Bodner streamfunction. However, because of the small Cr value we adopted (0.01), this had minimal effect on the solutions, so we decided to keep the original simpler implementation and not pursue this option for now. This commit should be evaluated together with CVMix/CVMix-src#106.

mnlevy1981

Sorry to sit on this for a month, but there are some minor updates needed. I'm happy to help you with them next week.

mnlevy1981 · 2025-08-15T14:54:10Z

src/shared/cvmix_kpp.F90

  public :: cvmix_kpp_ustokes_SL_model
-
+  public :: cvmix_kpp_compute_ER_depth
+  public :: cvmix_kpp_quad_fit


There's a polynomial fit in cvmix_math.F90, though it needs a little updating because the wants two points and a derivative instead of the three points used here. Still, it makes more sense to add this functionality to the existing routine rather than create a new subroutine here.

mnlevy1981 · 2025-08-15T14:56:48Z

src/shared/cvmix_kpp.F90

+    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


These should be declared as real(cvmix_r8), parameter (and it would be nice to see a comment describing units / where these values come from)

mnlevy1981 · 2025-08-15T14:59:09Z

src/shared/cvmix_kpp.F90

+    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 )


I think PBfact could also be a Fortran parameter (0.0893759_cvmix_r8) and then

PB = MAX( -surfBuoy_NS * BLdepth * PBfact , cvmix_zero ) & + PBfact * CempCGs * BLdepth * ( abs(surf_buoy_force) - surf_buoy_force )

PBfact cannot be used as a parameter because it is used with different values.

mnlevy1981 · 2025-08-15T15:00:12Z

src/shared/cvmix_kpp.F90

+
+  nlev      = size(OBL_depth)
+  Cemp_Rs   = 4.7_cvmix_r8
+!  Gat1      = cvmix_zero


We can just remove lines instead of commenting them out

mnlevy1981 · 2025-08-15T15:00:18Z

src/shared/cvmix_kpp.F90

+! kinv      = MAXLOC( Nsq )      ! interface index of maximum stratification, N2>0 (inversion)
+! maxNsq    = Nsq( kinv )


We can just remove lines instead of commenting them out

mnlevy1981 · 2025-08-15T15:01:04Z

src/shared/cvmix_kpp.F90

+  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


This is a duplicate of line 3591; should it be removed, or should it be BEnt(2) instead of BEnt(1)? Looking at how BEnt is used, I see a couple of BEnt(kbl+1) but no other index is evaluated. Maybe this should be a scalar instead of an array?

BEnt is now a scalar.

gustavo-marques · 2025-08-20T19:27:35Z

Thank you, @mnlevy1981, for all the suggestions. We have addressed all comments except for switching to the polynomial fit in cvmix_math.F90.

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

gustavo-marques · 2025-08-28T23:12:43Z

@mnlevy1981, here is the difference in SST when adding your polynomial fit after 1 month ( B = baseline):

This is more than I expected, so I will run for longer to check the effects on climate-related metrics.

gustavo-marques · 2025-09-04T16:10:16Z

This PR has minimal impact on climate metrics.
I approve it.

* Add entrainment rule BLD calculation to KPP This commit introduces an alternative method for computing the boundary layer depth (BLD) in KPP, based on the entrainment rule. This method is activated when the non-solar surface buoyancy flux (surfBuoy_NS) is negative and STOKES_MOST = True. In all other cases, the standard Richardson Number-based method is used. The non-solar surface buoyancy flux is estimated internally using an exponential attenuation function, since only the total surface buoyancy flux (solar + non-solar) is currently available to the module. Ideally, the non-solar component should be provided explicitly. Several new diagnostics have been added, and some of the existing ones have been updated. Specifically, the parameterized shear, buoyancy, and Stokes TKE production terms are calculated. We conducted multiple tests where these terms are passed to the MLE module and used in the denominator of the Bodner streamfunction. However, because of the small Cr value we adopted (0.01), this had minimal effect on the solutions, so we decided to keep the original simpler implementation and not pursue this option for now. This commit should be evaluated together with CVMix/CVMix-src#106. * Fix typos * Update CVMix to latest commit Include new entreinment rule depth changes. * Revert unit conversion for Vt2 and GoRho * Update description and remove empty line * Split description of variables * Refactor KPP StokesMOST diagnostics and allocations - Removed unused `ustar` array from KPP_CS type. - Restricted registration, allocation, and posting of StokesMOST diagnostics (`StokesXI`, `Lam2`, `BEdE_ER`, `ERdepth`, `RNdepth`, `PU_TKE`, `PS_TKE`, `PB_TKE`) to cases when `CS%StokesMOST` is enabled. - Moved allocations of related arrays inside StokesMOST conditional block. - Updated condition for calculating `surfBuoy_NS` to avoid division by zero in alog(buoyFlux(i,j,2)/buoyFlux(i,j,3)). - Wrapped diagnostic and debug checks in `if (CS%StokesMOST)` conditions. - Minor cleanup of comments and diagnostic labeling.

* Add entrainment rule BLD calculation to KPP This commit introduces an alternative method for computing the boundary layer depth (BLD) in KPP, based on the entrainment rule. This method is activated when the non-solar surface buoyancy flux (surfBuoy_NS) is negative and STOKES_MOST = True. In all other cases, the standard Richardson Number-based method is used. The non-solar surface buoyancy flux is estimated internally using an exponential attenuation function, since only the total surface buoyancy flux (solar + non-solar) is currently available to the module. Ideally, the non-solar component should be provided explicitly. Several new diagnostics have been added, and some of the existing ones have been updated. Specifically, the parameterized shear, buoyancy, and Stokes TKE production terms are calculated. We conducted multiple tests where these terms are passed to the MLE module and used in the denominator of the Bodner streamfunction. However, because of the small Cr value we adopted (0.01), this had minimal effect on the solutions, so we decided to keep the original simpler implementation and not pursue this option for now. This commit should be evaluated together with CVMix/CVMix-src#106. * Fix typos * Update CVMix to latest commit Include new entreinment rule depth changes. * Revert unit conversion for Vt2 and GoRho * Update description and remove empty line * Split description of variables * Refactor KPP StokesMOST diagnostics and allocations - Removed unused `ustar` array from KPP_CS type. - Restricted registration, allocation, and posting of StokesMOST diagnostics (`StokesXI`, `Lam2`, `BEdE_ER`, `ERdepth`, `RNdepth`, `PU_TKE`, `PS_TKE`, `PB_TKE`) to cases when `CS%StokesMOST` is enabled. - Moved allocations of related arrays inside StokesMOST conditional block. - Updated condition for calculating `surfBuoy_NS` to avoid division by zero in alog(buoyFlux(i,j,2)/buoyFlux(i,j,3)). - Wrapped diagnostic and debug checks in `if (CS%StokesMOST)` conditions. - Minor cleanup of comments and diagnostic labeling. * Add support for Lam2 in MLE and KPP - Introduced `Lam2` as an optional pointer field in `vertvisc_type` and restart registry. - Extended `mixedlayer_restrat` and `mixedlayer_restrat_Bodner` to accept/use `Lam2` when available, with fallback to equilibrium value. - Added `KPP_get_Lam2` routine to expose Lam2 from KPP control structure. - Connected `Lam2` into `diabatic_ALE_legacy` and `diabatic_ALE` and copied into `visc%Lam2` when allocated. - Enabled allocation and restart registration of `Lam2` when STOKES_MOST and Bodner MLE are used. - Updated MOM.F90 call sites to conditionally pass `Lam2`. This integrates Langmuir turbulence effects via Lam2 into MLE restratification. * Fix duplicate declaration of Lam2 * Change paramFile to param_file * Fixed OMP directive in KPP_get_Lam2 * Defined Lam2 in diabatic_ALE * Add ParameterBlock for KPP so STOKES_MOST can be read * Fix segfault by checking CS%visc before accessing Lam2 Wrapped the mixedlayer_restrat call with an outer check on `associated(CS%visc)` to ensure the derived type is valid before testing `CS%visc%Lam2`. This prevents invalid memory references when `CS%visc` is unassociated. * Replace associated with allocated * Replace associated with allocated * Replace allocated with associated for checking CS%visc%Lam2 * Only call KPP_get_Lam2 if visc%Lam2 is associated * Implement conditional logic for passing Lam2 to MLE. Read `StokesMOST` and `wave_enhanced_ustar` to determine whether Lam2 should be passed to the `mixedlayer_restrat` subroutine. * Safely handle optional Lam2 pointer in mixedlayer_restrat Added local logical flag `haveLam2` to ensure Lam2 is only used when both present and associated. This avoids potential segfaults from evaluating `associated(Lam2)` when Lam2 is not passed. * Improve description of WAVE_ENHANCED_USTAR

gustavo-marques added 6 commits July 18, 2025 10:58

Update stokes_most

61734e6

Delete unsed argument, deltaRho

268829e

Replace *SLD with *SL

81c4265

Replace uSbar_SL => uSb_SL; vSbar_SL => vSb_SL

09be4f4

Fix args passed to cvmix_kpp_driver

bb13605

Set SLDepth=cvmix_one and change BLDepth -> BLdepth

9f4da9c

gustavo-marques mentioned this pull request Aug 4, 2025

Add entrainment rule BLD calculation to KPP NCAR/MOM6#379

Merged

mnlevy1981 requested changes Aug 15, 2025

View reviewed changes

Addresses Mike Levy comments

20eee02

mnlevy1981 added 2 commits August 27, 2025 14:17

Fix whitespace

971b366

mnlevy1981 approved these changes Sep 9, 2025

View reviewed changes

mnlevy1981 merged commit d20b989 into CVMix:master Sep 9, 2025
1 check passed

		! kinv = MAXLOC( Nsq ) ! interface index of maximum stratification, N2>0 (inversion)
		! maxNsq = Nsq( kinv )

Update stokes most #106

Update stokes most #106

Uh oh!

Conversation

gustavo-marques commented Jul 18, 2025

Uh oh!

mnlevy1981 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gustavo-marques commented Aug 20, 2025

Uh oh!

gustavo-marques commented Aug 28, 2025

Uh oh!

gustavo-marques commented Sep 4, 2025

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants