Skip to content

Commit b213d2c

Browse files
aakashsaneaakashsane
committed
added ML_diffusivity shape function
1 parent 9187de9 commit b213d2c

1 file changed

Lines changed: 70 additions & 6 deletions

File tree

‎src/shared/cvmix_kpp.F90‎

Lines changed: 70 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,7 @@ module cvmix_kpp
6565
integer, parameter :: LANGMUIR_ENTRAINMENT_LWF16 = 1
6666
integer, parameter :: LANGMUIR_ENTRAINMENT_LF17 = 2
6767
integer, parameter :: LANGMUIR_ENTRAINMENT_RWHGK16 = 3
68+
integer, parameter :: ML_DIFFUSIVITY_SHAPE = 5 ! ML_Diffusivity
6869

6970
! !PUBLIC MEMBER FUNCTIONS:
7071

@@ -214,6 +215,9 @@ module cvmix_kpp
214215
real(cvmix_r8) :: ER_Cs ! Entrainment Rule TKE Stokes production weight [nondim]
215216
real(cvmix_r8) :: ER_Cu ! Entrainment Rule TKE shear production weight [nondim]
216217

218+
logical :: ML_diffusivity ! uses the Machine Learned equations to set diffusivity
219+
! from Sane et al. (2025) https://doi.org/10.31219/osf.io/uab7v_v2
220+
217221
end type cvmix_kpp_params_type
218222

219223
!EOP
@@ -232,7 +236,8 @@ subroutine cvmix_init_kpp(ri_crit, minOBLdepth, maxOBLdepth, minVtsqr, &
232236
old_vals, lEkman, lStokesMOST, lMonOb, lnoDGat1, &
233237
lenhanced_diff, lnonzero_surf_nonlocal, &
234238
Langmuir_mixing_str, Langmuir_entrainment_str, &
235-
l_LMD_ws, ER_Cb, ER_Cs, ER_Cu, CVmix_kpp_params_user)
239+
l_LMD_ws, ML_diffusivity, &
240+
ER_Cb, ER_Cs, ER_Cu, CVmix_kpp_params_user)
236241

237242
! !DESCRIPTION:
238243
! Initialization routine for KPP mixing.
@@ -267,6 +272,7 @@ subroutine cvmix_init_kpp(ri_crit, minOBLdepth, maxOBLdepth, minVtsqr, &
267272
lenhanced_diff, &
268273
lnonzero_surf_nonlocal, &
269274
l_LMD_ws
275+
logical, optional, intent(in) :: ML_diffusivity
270276

271277
! !OUTPUT PARAMETERS:
272278
type(cvmix_kpp_params_type), intent(inout), target, optional :: &
@@ -459,6 +465,8 @@ subroutine cvmix_init_kpp(ri_crit, minOBLdepth, maxOBLdepth, minVtsqr, &
459465
case ('ParabolicNonLocal')
460466
call cvmix_put_kpp('MatchTechnique', CVMIX_KPP_PARABOLIC_NONLOCAL, &
461467
CVmix_kpp_params_user)
468+
case ('ML_Diffusivity_Shape') ! ML_Diffusivity shape function
469+
call cvmix_put_kpp('MatchTechnique', ML_DIFFUSIVITY_SHAPE, CVmix_kpp_params_user)
462470
case DEFAULT
463471
print*, "ERROR: ", trim(MatchTechnique), " is not a valid choice ", &
464472
"for MatchTechnique!"
@@ -469,6 +477,16 @@ subroutine cvmix_init_kpp(ri_crit, minOBLdepth, maxOBLdepth, minVtsqr, &
469477
CVmix_kpp_params_user)
470478
end if
471479

480+
! --- Add this block to for ML_diffusivity ---
481+
if (present(ML_diffusivity)) then ! ML_diffusivity
482+
call cvmix_put_kpp('ML_diffusivity', ML_diffusivity, CVmix_kpp_params_user)
483+
if (ML_diffusivity) then
484+
call cvmix_put_kpp('MatchTechnique', ML_DIFFUSIVITY_SHAPE, CVmix_kpp_params_user)
485+
end if
486+
else
487+
call cvmix_put_kpp('ML_diffusivity', .false., CVmix_kpp_params_user)
488+
end if
489+
472490
if (present(old_vals)) then
473491
select case (trim(old_vals))
474492
case ("overwrite")
@@ -678,7 +696,7 @@ subroutine cvmix_coeffs_kpp_wrap(CVmix_vars, CVmix_kpp_params_user)
678696
CVmix_vars%SurfaceBuoyancyForcing, &
679697
nlev, max_nlev, &
680698
CVmix_vars%LangmuirEnhancementFactor, &
681-
CVmix_vars%StokesMostXi, &
699+
CVmix_vars%StokesMostXi, CVmix_vars%Coriolis, &
682700
CVmix_kpp_params_user)
683701

684702
call cvmix_update_wrap(CVmix_kpp_params_in%handle_old_vals, max_nlev, &
@@ -702,7 +720,7 @@ subroutine cvmix_coeffs_kpp_low(Mdiff_out, Tdiff_out, Sdiff_out, zw, zt, &
702720
old_Mdiff, old_Tdiff, old_Sdiff, OBL_depth, &
703721
kOBL_depth, Tnonlocal, Snonlocal, surf_fric,&
704722
surf_buoy, nlev, max_nlev, Langmuir_EFactor,&
705-
StokesXI,CVmix_kpp_params_user)
723+
StokesXI,Coriolis,CVmix_kpp_params_user)
706724

707725
! !DESCRIPTION:
708726
! Computes vertical diffusion coefficients for the KPP boundary layer mixing
@@ -727,6 +745,7 @@ subroutine cvmix_coeffs_kpp_low(Mdiff_out, Tdiff_out, Sdiff_out, zw, zt, &
727745
! Langmuir enhancement factor
728746
real(cvmix_r8), intent(in), optional :: Langmuir_EFactor
729747
real(cvmix_r8), intent(in), optional :: StokesXI
748+
real(cvmix_r8), intent(in), optional :: Coriolis ! required for ML_diffusivity
730749
! !INPUT/OUTPUT PARAMETERS:
731750
real(cvmix_r8), dimension(max_nlev+1), intent(inout) :: Mdiff_out, &
732751
Tdiff_out, &
@@ -794,6 +813,14 @@ subroutine cvmix_coeffs_kpp_low(Mdiff_out, Tdiff_out, Sdiff_out, zw, zt, &
794813
! Parameters for Stokes_MOST
795814
real(cvmix_r8) :: Gcomposite, Hsigma, sigh, T_NLenhance , S_NLenhance , XIone
796815

816+
! Parameters for the machine learning part, ML_diffusivity
817+
real(cvmix_r8) :: sigma_max ! sigma coordinate of location of maximum diffusivity
818+
real(cvmix_r8) :: g_sigma ! shape function at a sigma coordinate.
819+
real(cvmix_r8) :: L_h ! Non-dimensional L_h = B*OBL_depth/u_*^3
820+
real(cvmix_r8) :: E_h ! Non-dimensional E_h = OBL_depth * Coriolis /u_*
821+
real(cvmix_r8) :: F_inter_func ! Stands for F_intermediate_function,
822+
! Non-dimensional intermediate function used to calculate sigma_max
823+
797824
! Constant from params
798825
integer :: interp_type2, MatchTechnique
799826

@@ -1208,11 +1235,46 @@ subroutine cvmix_coeffs_kpp_low(Mdiff_out, Tdiff_out, Sdiff_out, zw, zt, &
12081235
call cvmix_kpp_compute_turbulent_scales(sigma, OBL_depth, surf_buoy, &
12091236
surf_fric, XIone, w_m, w_s, &
12101237
CVmix_kpp_params_user)
1238+
1239+
if (CVmix_kpp_params_in%ML_diffusivity) then ! ML_diffusivity
1240+
!!! Evaluating L_h and E_h
1241+
L_h = -surf_buoy * OBL_depth / (surf_fric ** 3.0) ! OBL / Monin-Obukhov-Depth
1242+
E_h = OBL_depth * Coriolis / surf_fric ! OBL / Ekman Depth i.e. hf/u*
1243+
F_inter_func = ( cvmix_one / ( 0.0712 + 0.4380 * exp(-1.0*(2.6821 * L_h)) ) ) + 1.5845
1244+
sigma_max = (F_inter_func * E_h) / ( 1.7908*(F_inter_func * E_h) + 0.6904)
1245+
!!! capping sigma_max between 0.1 and 0.7
1246+
sigma_max = min( max(sigma_max, 0.1), 0.7)
1247+
endif
12111248
do kw=2,kwup
12121249
! (3b) Evaluate G(sigma) at each cell interface
1213-
MshapeAtS = cvmix_math_evaluate_cubic(Mshape, sigma(kw))
1214-
TshapeAtS = cvmix_math_evaluate_cubic(Tshape, sigma(kw))
1215-
SshapeAtS = cvmix_math_evaluate_cubic(Sshape, sigma(kw))
1250+
1251+
if (CVmix_kpp_params_in%ML_diffusivity) then ! ML_diffusivity
1252+
1253+
! ML-diffusivity modification:
1254+
if (sigma(kw) .le. sigma_max ) then ! ML based shape function
1255+
! quadratic part above sigma_max
1256+
g_sigma = (2.0*sigma(kw)/sigma_max) - (sigma(kw)/sigma_max)**2.0
1257+
g_sigma = g_sigma * (4.0/27.0) ! reduces the amplitude from 1 to 4/27 to match G(\sigma) of KPP
1258+
MshapeAtS = g_sigma
1259+
TshapeAtS = g_sigma
1260+
SshapeAtS = g_sigma
1261+
else
1262+
! cubic part below sigma_max
1263+
g_sigma = 2.0*((sigma(kw) - sigma_max)/(cvmix_one - sigma_max))**3.0 &
1264+
- 3.0*((sigma(kw) - sigma_max)/(cvmix_one - sigma_max))**2.0 &
1265+
+ cvmix_one
1266+
g_sigma = g_sigma * (4.0/27.0) ! reduces the amplitude from 1 to 4/27 to match G(\sigma) of KPP
1267+
MshapeAtS = g_sigma
1268+
TshapeAtS = g_sigma
1269+
SshapeAtS = g_sigma
1270+
end if
1271+
1272+
else
1273+
1274+
MshapeAtS = cvmix_math_evaluate_cubic(Mshape, sigma(kw))
1275+
TshapeAtS = cvmix_math_evaluate_cubic(Tshape, sigma(kw))
1276+
SshapeAtS = cvmix_math_evaluate_cubic(Sshape, sigma(kw))
1277+
end if
12161278
! The RWHGK16 Langmuir uses the shape function to shape the
12171279
! enhancement to the mixing coefficient.
12181280
ShapeNoMatchAtS = cvmix_math_evaluate_cubic(NMshape, sigma(kw))
@@ -1483,6 +1545,8 @@ subroutine cvmix_put_kpp_logical(varname, val, CVmix_kpp_params_user)
14831545
CVmix_kpp_params_out%lenhanced_diff = val
14841546
case ('l_LMD_ws')
14851547
CVmix_kpp_params_out%l_LMD_ws = val
1548+
case ('ML_diffusivity')
1549+
CVmix_kpp_params_out%ML_diffusivity = val
14861550
case DEFAULT
14871551
print*, "ERROR: ", trim(varname), " is not a boolean variable!"
14881552
stop 1

0 commit comments

Comments
 (0)