Adds the Bodner et al. 2023 version of MLE

.testing/tc2.a/MOM_tc_variant

@@ -1,3 +1,9 @@
 #override TOPO_CONFIG = "spoon"
 #override REMAPPING_SCHEME = "PPM_H4"
 #override REGRIDDING_COORDINATE_MODE = "SIGMA"
+MLE_USE_PBL_MLD = True
+MLE%USE_BODNER23 = True
+MLE%BLD_DECAYING_TFILTER = 86400.
+MLE%MLD_DECAYING_TFILTER = 259200.
+MLE%BLD_GROWING_TFILTER = 300.
+MLE%MLD_GROWING_TFILTER = 3600.

src/core/MOM.F90

@@ -1263,7 +1263,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
     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%VarMix, G, GV, US, CS%mixedlayer_restrat_CSp)
+                            CS%visc%sfc_buoy_flx, CS%VarMix, G, GV, US, CS%mixedlayer_restrat_CSp)
     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

src/core/MOM_unit_tests.F90

@@ -11,6 +11,7 @@ module MOM_unit_tests
 use MOM_random,                     only : random_unit_tests
 use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests
 use MOM_CFC_cap,                    only : CFC_cap_unit_tests
+use MOM_mixed_layer_restrat,        only : mixedlayer_restrat_unit_tests
 implicit none ; private
 
 public unit_tests
@@ -40,6 +41,8 @@ subroutine unit_tests(verbosity)
        "MOM_unit_tests: near_boundary_unit_tests FAILED")
     if (CFC_cap_unit_tests(verbose)) call MOM_error(FATAL, &
        "MOM_unit_tests: CFC_cap_unit_tests FAILED")
+    if (mixedlayer_restrat_unit_tests(verbose)) call MOM_error(FATAL, &
+       "MOM_unit_tests: mixedlayer_restrat_unit_tests FAILED")
   endif
 
 end subroutine unit_tests

src/core/MOM_variables.F90

@@ -257,8 +257,8 @@ module MOM_variables
     Ray_v       !< The Rayleigh drag velocity to be applied to each layer at v-points [Z T-1 ~> m s-1].
 
   ! 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(:,:) :: MLD => NULL() !< Instantaneous active mixing layer depth [Z ~> m].
+  real, pointer, dimension(:,:) :: sfc_buoy_flx !< Surface buoyancy flux (derived) [Z2 T-3 ~> m2 s-3].
   real, pointer, dimension(:,:,:) :: Kd_shear => NULL()
                 !< The shear-driven turbulent diapycnal diffusivity at the interfaces between layers
                 !! in tracer columns [Z2 T-1 ~> m2 s-1].

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -712,6 +712,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
       ! If visc%MLD exists, copy KPP's BLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
     endif
+    if (associated(visc%sfc_buoy_flx)) then
+      visc%sfc_buoy_flx(:,:) = CS%KPP_buoy_flux(:,:,1)
+      call pass_var(visc%sfc_buoy_flx, G%domain, halo=1)
+    endif
 
     if (.not.CS%KPPisPassive) then
       !$OMP parallel do default(shared)
@@ -854,6 +858,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
       call energetic_PBL_get_MLD(CS%ePBL, visc%MLD, G, US)
       call pass_var(visc%MLD, G%domain, halo=1)
     endif
+    if (associated(visc%sfc_buoy_flx)) then
+      visc%sfc_buoy_flx(:,:) = SkinBuoyFlux(:,:)
+      call pass_var(visc%sfc_buoy_flx, G%domain, halo=1)
+    endif
 
     ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
     do K=2,nz ; do j=js,je ; do i=is,ie
@@ -1306,6 +1314,10 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
       ! If visc%MLD exists, copy KPP's BLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
     endif
+    if (associated(visc%sfc_buoy_flx)) then
+      visc%sfc_buoy_flx(:,:) = CS%KPP_buoy_flux(:,:,1)
+      call pass_var(visc%sfc_buoy_flx, G%domain, halo=1)
+    endif
 
     if (showCallTree) call callTree_waypoint("done with KPP_calculate (diabatic)")
     if (CS%debug) then
@@ -1391,6 +1403,10 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
       call energetic_PBL_get_MLD(CS%ePBL, visc%MLD, G, US)
       call pass_var(visc%MLD, G%domain, halo=1)
     endif
+    if (associated(visc%sfc_buoy_flx)) then
+      visc%sfc_buoy_flx(:,:) = SkinBuoyFlux(:,:)
+      call pass_var(visc%sfc_buoy_flx, G%domain, halo=1)
+    endif
 
     ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
     do K=2,nz ; do j=js,je ; do i=is,ie
@@ -1900,6 +1916,10 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
       ! If visc%MLD exists, copy KPP's BLD into it
       if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
     endif
+    if (associated(visc%sfc_buoy_flx)) then
+      visc%sfc_buoy_flx(:,:) = CS%KPP_buoy_flux(:,:,1)
+      call pass_var(visc%sfc_buoy_flx, G%domain, halo=1)
+    endif
 
     if (.not. CS%KPPisPassive) then
       !$OMP parallel do default(shared)

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -1870,7 +1870,7 @@ subroutine set_visc_register_restarts(HI, GV, US, param_file, visc, restart_CS)
   ! Local variables
   logical :: use_kappa_shear, KS_at_vertex
   logical :: adiabatic, useKPP, useEPBL
-  logical :: use_CVMix_shear, MLE_use_PBL_MLD, use_CVMix_conv
+  logical :: use_CVMix_shear, MLE_use_PBL_MLD, MLE_use_Bodner, use_CVMix_conv
   integer :: isd, ied, jsd, jed, nz
   real :: hfreeze !< If hfreeze > 0 [Z ~> m], melt potential will be computed.
   character(len=40)  :: mdl = "MOM_set_visc"  ! This module's name.
@@ -1942,6 +1942,15 @@ subroutine set_visc_register_restarts(HI, GV, US, param_file, visc, restart_CS)
     call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
   endif
 
+  ! visc%sfc_buoy_flx is used to communicate the state of the (e)PBL or KPP to the rest of the model
+  call get_param(param_file, mdl, "MLE%USE_BODNER23", MLE_use_Bodner, &
+                 default=.false., do_not_log=.true.)
+  if (MLE_use_PBL_MLD .or. MLE_use_Bodner) then
+    call safe_alloc_ptr(visc%sfc_buoy_flx, isd, ied, jsd, jed)
+    call register_restart_field(visc%sfc_buoy_flx, "SFC_BFLX", .false., restart_CS, &
+                                "Instantaneous surface buoyancy flux", "m2 s-3", &
+                                conversion=US%Z_to_m**2*US%s_to_T**3)
+  endif
 
 end subroutine set_visc_register_restarts