added ML_diffusivity shape function #111
Open
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
aakashsane
force-pushed
the
ML_diffusivity_3
branch
from
b213d2c to
7b8f824
Compare
mnlevy1981
reviewed
Jan 29, 2026
src/shared/cvmix_kpp.F90
Outdated
Comment on lines
3750
to
3881
|
|
||
| !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,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] | ||
| 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 | ||
|
|
||
| ! Local variables | ||
| integer :: iz, nlev , kbl , kinv | ||
| 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, BEnt | ||
| real(kind=cvmix_r8), dimension(4) :: coeffs | ||
| type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in | ||
|
|
||
| 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 | ||
| Fdelrho = cvmix_one | ||
| 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 | ||
| sigma(:) = cvmix_zero | ||
| Bflux(1) = Bsfc_ns ! non-solar surface buoyancy boundary condition for all kbl | ||
| ! Set OBL_depth(1)=top cell center values | ||
| 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) .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)) | ||
| ! comment by Aakash: the above is the original line, | ||
| ! it gives me errors so I changed it to below call. not sure if it is correct | ||
| call cvmix_math_poly_interp(coeffs, CVMIX_MATH_INTERP_QUAD, & | ||
| sigma(iz:iz+1), Bflux(iz:iz+1), & | ||
| sigma(iz+2), Bflux(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 | ||
| 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_ER_depth | ||
|
|
Contributor
There was a problem hiding this comment.
Your merge of the latest master introduced a second copy of cvmix_kpp_compute_ER_depth(). Did you change this function in a previous commit? I'd suggest deleting the first instance of this function, in case there are differences in this second copy...
mnlevy1981
requested changes
Jan 30, 2026
| call cvmix_math_poly_interp(coeffs, sigma(iz:iz+2), Bflux(iz:iz+2)) | ||
| ! call cvmix_math_poly_interp(coeffs, sigma(iz:iz+2), Bflux(iz:iz+2)) | ||
| ! comment by Aakash: the above is the original line, | ||
| ! it gives me errors so I changed it to below call. not sure if it is correct |
Contributor
There was a problem hiding this comment.
What error did you get? These two calls to cvmix_math_poly_interp are not equivalent, and will change answers regardless of the value of ML_diffusivity:
call cvmix_math_poly_interp(coeffs, sigma(iz:iz+2), Bflux(iz:iz+2))fits a quadratic such that P(sigma(i)) = Bflux(i) for i = iz, iz+1, iz+2 (standard quadratic interpolation)
call cvmix_math_poly_interp(coeffs, CVMIX_MATH_INTERP_QUAD, &
sigma(iz:iz+1), Bflux(iz:iz+1), &
sigma(iz+2), Bflux(iz+2))fits a quadratic such that P(sigma(i)) = Bflux(i) for i = iz, iz+1, and P'(sigma(iz)) = (Bflux(iz) - Bflux(iz+2)) / (sigma(iz) - sigma(iz+2)). If we were okay with this interpolation change, you'd still need to change your arguments to
call cvmix_math_poly_interp(coeffs, CVMIX_MATH_INTERP_QUAD, &
sigma(iz+1:iz+2), Bflux(iz+1:iz+2), &
sigma(iz), Bflux(iz))And this would interpolate at iz+1 and iz+2 while having the derivative at iz+1 match the finite difference computed between iz and iz+1...
Anyway, if you could print out the three sigma values and three Bflux values that cause problems in your runs then I'll look into what's going wrong in the interpolation routine we want to use.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Adds option to activate shape function obtained from machine learning (genetic programming and least squares fitting).
The relevant flag is 'CVmix_kpp_params_in%ML_diffusivity', when set to True uses the mixing coefficients.
The mixing coefficients are obtained using equations similar to the ones in the GFDL's ePBL scheme.
The corresponding article preprint is: Sane et al. (2025) https://doi.org/10.31219/osf.io/uab7v_v2