Option to add effects of Lam2 into MLE restratification.

src/core/MOM.F90

@@ -315,6 +315,10 @@ module MOM
   logical :: useMEKE                 !< If true, call the MEKE parameterization.
   logical :: use_stochastic_EOS      !< If true, use the stochastic EOS parameterizations.
   logical :: useWaves                !< If true, update Stokes drift
+  logical :: StokesMOST              !< If true, use Stokes Similarity package. Needed to decide if Lam2 should
+                                     !! be passed to mixedlayer_restrat.
+  logical :: wave_enhanced_ustar     !< If true, enhance ustar in Bodner23. Needed to decide if Lam2 should
+                                     !! be passed to mixedlayer_restrat.
   real :: dtbt_reset_period          !< The time interval between dynamic recalculation of the
                                      !! barotropic time step [T ~> s]. If this is negative dtbt is never
                                      !! calculated, and if it is 0, dtbt is calculated every step.
@@ -1427,8 +1431,17 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, &
                     CS%uhtr, CS%vhtr, G%HI, haloshift=0, unscale=GV%H_to_MKS*US%L_to_m**2)
     endif
     call cpu_clock_begin(id_clock_ml_restrat)
-    call mixedlayer_restrat(h, CS%uhtr, CS%vhtr, CS%tv, forces, dt, CS%visc%MLD, CS%visc%h_ML, &
-                            CS%visc%sfc_buoy_flx, CS%VarMix, G, GV, US, CS%mixedlayer_restrat_CSp)
+    if (CS%wave_enhanced_ustar .and. CS%StokesMOST) then
+      if (associated(CS%visc%Lam2)) then
+        call mixedlayer_restrat(h, CS%uhtr, CS%vhtr, CS%tv, forces, dt, CS%visc%MLD, CS%visc%h_ML, &
+                    CS%visc%sfc_buoy_flx, CS%VarMix, G, GV, US, CS%mixedlayer_restrat_CSp, CS%visc%Lam2)
+      else
+        call MOM_error(FATAL,'step_MOM_dynamics:CS%visc%Lam2 not associated')
+      endif
+    else
+      call mixedlayer_restrat(h, CS%uhtr, CS%vhtr, CS%tv, forces, dt, CS%visc%MLD, CS%visc%h_ML, &
+                  CS%visc%sfc_buoy_flx, CS%VarMix, G, GV, US, CS%mixedlayer_restrat_CSp)
+    endif
     call cpu_clock_end(id_clock_ml_restrat)
     call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil))
     if (CS%debug) then
@@ -2449,12 +2462,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   call get_param(param_file, '', "FPMIX", fpmix, &
                  "If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", &
                  default=.false., do_not_log=.true.)
-
   if (fpmix .and. .not. CS%split)  then
     call MOM_error(FATAL, "initialize_MOM: "//&
        "FPMIX=True only works when SPLIT=True.")
   endif
-
+  ! STOKES_MOST and  needed to
+  call get_param(param_file, '', 'STOKES_MOST', CS%StokesMOST, &
+                 'If True, use Stokes Similarity package.', &
+                 default=.False., do_not_log=.true.)
+  call get_param(param_file, '', "WAVE_ENHANCED_USTAR", CS%wave_enhanced_ustar, &
+             "If true, enhance ustar in Bodner23.", default=.false., do_not_log=.true.)
   call get_param(param_file, "MOM", "BOUSSINESQ", Boussinesq, &
                  "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.)
   call get_param(param_file, "MOM", "SEMI_BOUSSINESQ", semi_Boussinesq, &

src/core/MOM_variables.F90

@@ -291,6 +291,7 @@ module MOM_variables
 
   ! The following elements are pointers so they can be used as targets for pointers in the restart registry.
   real, pointer, dimension(:,:) :: MLD => NULL()  !< Instantaneous active mixing layer depth [Z ~> m].
+  real, pointer, dimension(:,:) :: Lam2 => NULL() !< (Langmuir Number)^-2  [nondim].
   real, pointer, dimension(:,:) :: h_ML => NULL() !< Instantaneous active mixing layer thickness [H ~> m or kg m-2].
   real, pointer, dimension(:,:) :: sfc_buoy_flx => NULL() !< Surface buoyancy flux (derived) [Z2 T-3 ~> m2 s-3].
   real, pointer, dimension(:,:,:) :: Kd_shear => NULL()

src/parameterizations/lateral/MOM_mixed_layer_restrat.F90

@@ -102,8 +102,9 @@ module MOM_mixed_layer_restrat
                                    !! front-length scales read from a file.
   type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
   logical :: use_Stanley_ML        !< If true, use the Stanley parameterization of SGS T variance
-  logical :: wave_enhanced_ustar   !< If true, enhance ustar for equilibrium surface waves (La-2=11),
-                                   !! following Eq. 28 in Bodner23.
+  logical :: wave_enhanced_ustar   !< If true, enhance ustar using surface waves, following Eq. 28 in Bodner23.
+                                   !! Use a Langmuir number if provided. Otherwise, assumes equilibrium
+                                   !! surface waves (La-2=11.).
   real    :: ustar_min             !< A minimum value of ustar in thickness units to avoid numerical
                                    !! problems [H T-1 ~> m s-1 or kg m-2 s-1]
 
@@ -149,7 +150,7 @@ module MOM_mixed_layer_restrat
 !> Driver for the mixed-layer restratification parameterization.
 !! The code branches between two different implementations depending
 !! on whether the bulk-mixed layer or a general coordinate are in use.
-subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux, VarMix, G, GV, US, CS)
+subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux, VarMix, G, GV, US, CS, Lam2)
   type(ocean_grid_type),                      intent(inout) :: G      !< Ocean grid structure
   type(verticalGrid_type),                    intent(in)    :: GV     !< Ocean vertical grid structure
   type(unit_scale_type),                      intent(in)    :: US     !< A dimensional unit scaling type
@@ -170,17 +171,29 @@ subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux,
                                                                       !! PBL scheme [Z2 T-3 ~> m2 s-3]
   type(VarMix_CS),                            intent(in)    :: VarMix !< Variable mixing control structure
   type(mixedlayer_restrat_CS),                intent(inout) :: CS     !< Module control structure
+  real, dimension(:,:),             optional, pointer       :: Lam2   !< (Langmuir Number)^-2  [nondim]
 
 
+  ! local variables
+  logical :: haveLam2 !< True if optional Lam2 argument is both present and associated
+
   if (.not. CS%initialized) call MOM_error(FATAL, "mixedlayer_restrat: "// &
          "Module must be initialized before it is used.")
 
+  ! Determine if Lam2 should be used
+  haveLam2 = .false.
+  if (present(Lam2)) haveLam2 = associated(Lam2)
+
   if (GV%nkml>0) then
     ! Original form, written for the isopycnal model with a bulk mixed layer
     call mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
   elseif (CS%use_Bodner) then
     ! Implementation of Bodner et al., 2023
-    call mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux)
+    if (haveLam2) then
+      call mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux, Lam2)
+    else
+      call mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux)
+    endif
   else
     ! Implementation of Fox-Kemper et al., 2008, to work in general coordinates
     call mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, G, GV, US, CS)
@@ -755,7 +768,7 @@ end function mu
 
 !> Calculates a restratifying flow in the mixed layer, following the formulation
 !! used in Bodner et al., 2023 (B22)
-subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, dt, BLD, h_MLD, bflux)
+subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, dt, BLD, h_MLD, bflux, Lam2)
   ! Arguments
   type(mixedlayer_restrat_CS),                intent(inout) :: CS     !< Module control structure
   type(ocean_grid_type),                      intent(inout) :: G      !< Ocean grid structure
@@ -776,6 +789,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
                                                                       !! the PBL scheme [H ~> m or kg m-2]
   real, dimension(:,:),                       pointer       :: bflux  !< Surface buoyancy flux provided by the
                                                                       !! PBL scheme [Z2 T-3 ~> m2 s-3]
+  real, dimension(:,:),             optional, pointer       :: Lam2   !< (Langmuir Number)^-2, which is defined as
+                                                                      !! Surface Stokes/ustar [nondim]
+
   ! Local variables
   real :: uhml(SZIB_(G),SZJ_(G),SZK_(GV)) ! zonal mixed layer transport [H L2 T-1 ~> m3 s-1 or kg s-1]
   real :: vhml(SZI_(G),SZJB_(G),SZK_(GV)) ! merid mixed layer transport [H L2 T-1 ~> m3 s-1 or kg s-1]
@@ -809,7 +825,6 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
   real :: w_star3         ! Cube of turbulent convective velocity [Z3 T-3 ~> m3 s-3]
   real :: u_star3         ! Cube of surface friction velocity [Z3 T-3 ~> m3 s-3]
   real :: E_ustar         ! Surface wave ustar enhancement factor   [nondim]
-  real :: Lam2            ! Reciprocal of the squrared turbulent Langmuir number [nondim]
   real :: r_wpup          ! reciprocal of vertical momentum flux [T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1]
   real :: absf            ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
   real :: f_h             ! Coriolis parameter at h-points [T-1 ~> s-1]
@@ -833,6 +848,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
                           ! fractional power [T3 m3 Z-3 s-3 ~> 1]
   real :: m2_s2_to_Z2_T2  ! Conversion factors to restore scaling after a term is raised to a
                           ! fractional power [Z2 s2 T-2 m-2 ~> 1]
+  real, parameter :: Lam2_eq = 11.       ! (Langmuir Number)^-2 assuming wind wave equilibrium [nondim]
   real, parameter :: two_thirds = 2./3.  ! [nondim]
   logical :: line_is_empty, keep_going
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
@@ -873,15 +889,20 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
   ! Extract the friction velocity from the forcing type.
   call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1)
 
-  ! wave enhancement of ustar following Eq. 28 in Bodner23
+  ! Wave Enhanced of ustar following Eq. 28 in Bodner23
   if (CS%wave_enhanced_ustar) then
-    ! Assuming wind wave equilibrium (Lam2=11) until Lam2 becomes available
-    Lam2 = 11.
-    E_ustar   =  sqrt( 1.0 + (Lam2 * 0.104) + (Lam2 * Lam2 * 0.00118))
-    ! Wave Enhanced
-    do j=js-1,je+1 ; do i=is-1,ie+1
-      U_star_2d(i,j) = E_ustar * U_star_2d(i,j)
-    enddo ; enddo
+    if (present(Lam2) .and. associated(Lam2)) then
+      do j=js-1,je+1 ; do i=is-1,ie+1
+        E_ustar   =  sqrt( 1.0 + (Lam2(i,j) * 0.104) + (Lam2(i,j) * Lam2(i,j) * 0.00118))
+        U_star_2d(i,j) = E_ustar * U_star_2d(i,j)
+      enddo ; enddo
+    else
+      ! Assuming wind wave equilibrium (Lam2=11)
+      E_ustar   =  sqrt( 1.0 + (Lam2_eq * 0.104) + (Lam2_eq * Lam2_eq * 0.00118))
+      do j=js-1,je+1 ; do i=is-1,ie+1
+        U_star_2d(i,j) = E_ustar * U_star_2d(i,j)
+      enddo ; enddo
+    endif
   endif
 
   if (CS%debug) then
@@ -1755,8 +1776,9 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
              "parameter a micron away from the equator.", &
              units="m2 s-2", default=1.0e-24, scale=US%m_to_Z**2*US%T_to_s**2)
     call get_param(param_file, mdl, "WAVE_ENHANCED_USTAR", CS%wave_enhanced_ustar, &
-             "If true, enhance ustar for equilibrium surface waves (La-2=11.), "// &
-             "following Eq. 28 in Bodner23.", default=.false.)
+             "If true, enhance ustar using surface waves, following Eq. 28 in Bodner23. " //&
+             "Use a Langmuir number if provided. Otherwise, assumes equilibrium "// &
+             "surface waves (La-2=11.).", default=.false.)
     call get_param(param_file, mdl, "TAIL_DH", CS%MLE_tail_dh, &
              "Fraction by which to extend the mixed-layer restratification "//&
              "depth used for a smoother stream function at the base of "//&

src/parameterizations/vertical/MOM_CVMix_KPP.F90

@@ -47,6 +47,8 @@ module MOM_CVMix_KPP
 public :: KPP_NonLocalTransport_saln
 public :: KPP_NonLocalTransport
 public :: KPP_get_BLD
+public :: KPP_get_Lam2
+
 
 ! Enumerated constants
 integer, private, parameter :: NLT_SHAPE_CVMix     = 0 !< Use the CVMix profile
@@ -161,6 +163,7 @@ module MOM_CVMix_KPP
 
   ! Diagnostics arrays
   real, pointer,     dimension(:,:)   :: OBLdepth  !< Depth (positive) of ocean boundary layer (OBL) [Z ~> m]
+  real, pointer,     dimension(:,:)   :: Lam2      !< La^(-2) = Ustk0/u* [nondim]
   real, allocatable, dimension(:,:,:) :: BulkRi    !< Bulk Richardson number for each layer [nondim]
   real, allocatable, dimension(:,:,:) :: N         !< Brunt-Vaisala frequency [T-1 ~> s-1]
   real, allocatable, dimension(:,:,:) :: N2        !< Squared Brunt-Vaisala frequency [T-2 ~> s-2]
@@ -184,7 +187,6 @@ module MOM_CVMix_KPP
   real, allocatable, dimension(:,:)   :: RNdepth   !< Percent use Ri Number boundary layer depth [nondim]
   real, allocatable, dimension(:,:)   :: StokesXI  !< Stokes similarity parameter  [nondim]
   real, allocatable, dimension(:,:)   :: BEdE_ER   !< Enrtainment Rule's Parameterized BEdE [ m3 s-3 ]
-  real, allocatable, dimension(:,:)   :: Lam2      !< La^(-2) = Ustk0/u*
   ! Other arrays
   real, allocatable, dimension(:,:)   :: kOBL      !< Level (+fraction) of OBL extent [nondim]
   real, allocatable, dimension(:,:)   :: OBLdepthprev !< previous Depth (positive) of OBL [Z ~> m]
@@ -1740,6 +1742,24 @@ subroutine KPP_get_BLD(CS, BLD, G, US, m_to_BLD_units)
 
 end subroutine KPP_get_BLD
 
+!> Copies CS%Lam2 into Lam2.
+subroutine KPP_get_Lam2(CS, Lam2, G, US)
+  type(KPP_CS),                     pointer     :: CS  !< Control structure for
+                                                       !! this module
+  type(ocean_grid_type),            intent(in)  :: G   !< Grid structure
+  type(unit_scale_type),            intent(in)  :: US  !< A dimensional unit scaling type
+  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: Lam2 !< (Langmuir Number)^-2 [nondim]
+
+  ! Local variables
+  integer :: i,j ! Horizontal indices
+
+  !$OMP parallel do default(none) shared(Lam2, CS, G)
+  do j = G%jsc, G%jec ; do i = G%isc, G%iec
+    Lam2(i,j) = CS%Lam2(i,j)
+  enddo ; enddo
+
+end subroutine KPP_get_Lam2
+
 !> Apply KPP non-local transport of surface fluxes for a given tracer
 subroutine KPP_NonLocalTransport(CS, G, GV, h, nonLocalTrans, surfFlux, &
                                  dt, diag, tr_ptr, scalar, flux_scale)

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -53,7 +53,7 @@ module MOM_diabatic_driver
 use MOM_internal_tides,      only : internal_tides_init, internal_tides_end, int_tide_CS
 use MOM_kappa_shear,         only : kappa_shear_is_used
 use MOM_CVMix_KPP,           only : KPP_CS, KPP_init, KPP_compute_BLD, KPP_calculate
-use MOM_CVMix_KPP,           only : KPP_end, KPP_get_BLD, register_KPP_restarts
+use MOM_CVMix_KPP,           only : KPP_end, KPP_get_BLD, register_KPP_restarts, KPP_get_Lam2
 use MOM_CVMix_KPP,           only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln
 use MOM_oda_incupd,          only : apply_oda_incupd, oda_incupd_CS
 use MOM_opacity,             only : opacity_init, opacity_end, opacity_CS
@@ -591,6 +591,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
 
   real, dimension(SZI_(G),SZJ_(G)) :: &
     U_star, &    ! The friction velocity [Z T-1 ~> m s-1].
+    Lam2,   &    !  (Langmuir Number)^-2 [nondim]
     KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
     SkinBuoyFlux, &  ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
@@ -784,10 +785,14 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     endif
 
     call KPP_get_BLD(CS%KPP_CSp, BLD(:,:), G, US)
+    if (associated(visc%Lam2)) then
+      call KPP_get_Lam2(CS%KPP_CSp, Lam2(:,:), G, US)
+    endif
     ! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.
     if (associated(visc%h_ML)) call convert_MLD_to_ML_thickness(BLD, h, visc%h_ML, tv, G, GV)
     if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
     if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1)
+    if (associated(visc%Lam2)) visc%Lam2(:,:) = Lam2(:,:)
 
     if (.not.CS%KPPisPassive) then
       !$OMP parallel do default(shared)
@@ -1321,6 +1326,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
 
   real, dimension(SZI_(G),SZJ_(G)) :: &
     U_star, &    ! The friction velocity [Z T-1 ~> m s-1].
+    Lam2,   &    !  (Langmuir Number)^-2 [nondim]
     KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
     KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
     SkinBuoyFlux, &  ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
@@ -1517,11 +1523,16 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
                        Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves)
     endif
 
+
     call KPP_get_BLD(CS%KPP_CSp, BLD(:,:), G, US)
+    if (associated(visc%Lam2)) then
+      call KPP_get_Lam2(CS%KPP_CSp, Lam2(:,:), G, US)
+    endif
     ! If visc%MLD or visc%h_ML exist, copy KPP's BLD into them with appropriate conversions.
     if (associated(visc%h_ML)) call convert_MLD_to_ML_thickness(BLD, h, visc%h_ML, tv, G, GV)
     if (associated(visc%MLD)) visc%MLD(:,:) = BLD(:,:)
     if (associated(visc%sfc_buoy_flx)) visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1) * US%L_to_Z**2
+    if (associated(visc%Lam2)) visc%Lam2(:,:) = Lam2(:,:)
 
     if (showCallTree) call callTree_waypoint("done with KPP_calculate (diabatic)")
     if (CS%debug) then