+*Convert MLD to ML_thickness in diabatic

  Moved the call to convert mixed layer depths into mixed layer thicknesses
from right before the call to mixedlayer_restrat into the various diabatic
routines just after they are calculated.  This code rearrangement changes
answers in non-Boussinesq mode because the layer specific volumes will have
evolved between these two calls, but it is consistent with the mixed
layer thicknesses as used in the boundary layer parameterizations.

  A new argument was added to bulkmixedlayer to return the mixed layer
thickness that was already being calculated.  The previous argument Hml
was renamed to BLD and changed from a pointer into a simple array.  Both
are intent(inout) rather than intent(out) to preserve values in halos.

  This commit also revises the logic around the allocation of visc%h_ML and
its registration as a restart variable, including handling cases where an older
restart file is being read that includes MLD but not h_ML.  Because the mixed
layer depth is used in almost cases, CS%Hml in the control structure for
MOM.F90 was changed from a pointer to an allocatable, with values of 0 in those
cases (e.g., when ADIABATIC is true) where it is not used.

  In addition, the argument Hml to the various diabatic routines was changed to
MLD to more clearly reflect that it is a depth and not a thickness.  Outside of
diabatic, BML is only used to provide the tracer point boundary layer depths
to the calling routines, so it does not need a halo update.  Instead, all of
the halo updates for the various elements of visc that do need halo updates
after the diabatic calls are collected into one place for efficiency and more
accurate timings.

  The scale arguments of 34 checksum calls were revised from GV%H_to_m to
GV%H_to_MKS for more accurate checksums in non-Boussinesq configurations by
avoiding multiplication by an reference specific volume, and instead only
rescaling by an integer power of 2.

  All answers are bitwise identical in Boussinesq mode, but in non-Boussinesq
mode there are changes in answers in cases that use the boundary layer
thicknesses obtained from the boundary layer parameterizations in subsequent
calculations, such as mixed layer restratification with some options.  There are
also changes to the arguments of several publicly visible routines.

src/core/MOM.F90

@@ -97,7 +97,6 @@ module MOM
 use MOM_hor_index,             only : hor_index_type, hor_index_init
 use MOM_hor_index,             only : rotate_hor_index
 use MOM_interface_heights,     only : find_eta, calc_derived_thermo, thickness_to_dz
-use MOM_interface_heights,     only : convert_MLD_to_ML_thickness
 use MOM_interface_filter,      only : interface_filter, interface_filter_init, interface_filter_end
 use MOM_interface_filter,      only : interface_filter_CS
 use MOM_internal_tides,        only : int_tide_CS
@@ -211,8 +210,8 @@ module MOM
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta_av_bc
                     !< free surface height or column mass time averaged over the last
                     !! baroclinic dynamics time step [H ~> m or kg m-2]
-  real, dimension(:,:), pointer :: &
-    Hml => NULL()   !< active mixed layer depth [Z ~> m]
+  real, dimension(:,:), pointer :: Hml => NULL()
+                    !< active mixed layer depth, or 0 if there is no boundary layer scheme [Z ~> m]
   real :: time_in_cycle !< The running time of the current time-stepping cycle
                     !! in calls that step the dynamics, and also the length of
                     !! the time integral of ssh_rint [T ~> s].
@@ -1327,10 +1326,6 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
                     CS%uhtr, CS%vhtr, G%HI, haloshift=0, scale=GV%H_to_MKS*US%L_to_m**2)
     endif
     call cpu_clock_begin(id_clock_ml_restrat)
-    if (associated(CS%visc%MLD)) then
-      call safe_alloc_ptr(CS%visc%h_ML, G%isd, G%ied, G%jsd, G%jed)
-      call convert_MLD_to_ML_thickness(CS%visc%MLD, h, CS%visc%h_ML, CS%tv, G, GV, halo=1)
-    endif
     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)
     call cpu_clock_end(id_clock_ml_restrat)
@@ -2128,6 +2123,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   logical :: Boussinesq        ! If true, this run is fully Boussinesq
   logical :: semi_Boussinesq   ! If true, this run is partially non-Boussinesq
   logical :: use_KPP           ! If true, diabatic is using KPP vertical mixing
+  logical :: MLE_use_PBL_MLD   ! If true, use stored boundary layer depths for submesoscale restratification.
   integer :: nkml, nkbl, verbosity, write_geom
   integer :: dynamics_stencil  ! The computational stencil for the calculations
                                ! in the dynamic core.
@@ -2732,7 +2728,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   if (use_frazil) allocate(CS%tv%frazil(isd:ied,jsd:jed), source=0.0)
   if (bound_salinity) allocate(CS%tv%salt_deficit(isd:ied,jsd:jed), source=0.0)
 
-  if (bulkmixedlayer .or. use_temperature) allocate(CS%Hml(isd:ied,jsd:jed), source=0.0)
+  allocate(CS%Hml(isd:ied,jsd:jed), source=0.0)
 
   if (bulkmixedlayer) then
     GV%nkml = nkml ; GV%nk_rho_varies = nkml + nkbl
@@ -3271,11 +3267,23 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   CS%mixedlayer_restrat = mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, &
                                                   CS%mixedlayer_restrat_CSp, restart_CSp)
   if (CS%mixedlayer_restrat) then
+    if (GV%Boussinesq .and. associated(CS%visc%h_ML)) then
+      ! This is here to allow for a transition of restart files between model versions.
+      call get_param(param_file, "MOM", "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
+                     default=.false., do_not_log=.true.)
+      if (MLE_use_PBL_MLD .and. .not.query_initialized(CS%visc%h_ML, "h_ML", restart_CSp) .and. &
+          associated(CS%visc%MLD)) then
+        do j=js,je ; do i=is,ie ; CS%visc%h_ML(i,j) = GV%Z_to_H * CS%visc%MLD(i,j) ; enddo ; enddo
+      endif
+    endif
+
     if (.not.(bulkmixedlayer .or. CS%use_ALE_algorithm)) &
       call MOM_error(FATAL, "MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
     ! When DIABATIC_FIRST=False and using CS%visc%ML in mixedlayer_restrat we need to update after a restart
     if (.not. CS%diabatic_first .and. associated(CS%visc%MLD)) &
       call pass_var(CS%visc%MLD, G%domain, halo=1)
+    if (.not. CS%diabatic_first .and. associated(CS%visc%h_ML)) &
+      call pass_var(CS%visc%h_ML, G%domain, halo=1)
   endif
 
   call MOM_diagnostics_init(MOM_internal_state, CS%ADp, CS%CDp, Time, G, GV, US, &
@@ -3567,7 +3575,7 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp)
   ! hML is needed when using the ice shelf module
   call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., &
                  do_not_log=.true.)
-  if (use_ice_shelf .and. associated(CS%Hml)) then
+  if (use_ice_shelf) then
     call register_restart_field(CS%Hml, "hML", .false., restart_CSp, &
                                 "Mixed layer thickness", "m", conversion=US%Z_to_m)
   endif
@@ -3707,11 +3715,9 @@ subroutine extract_surface_state(CS, sfc_state_in)
   enddo ; enddo ; endif
 
   ! copy Hml into sfc_state, so that caps can access it
-  if (associated(CS%Hml)) then
-    do j=js,je ; do i=is,ie
-      sfc_state%Hml(i,j) = CS%Hml(i,j)
-    enddo ; enddo
-  endif
+  do j=js,je ; do i=is,ie
+    sfc_state%Hml(i,j) = CS%Hml(i,j)
+  enddo ; enddo
 
   if (CS%Hmix < 0.0) then  ! A bulk mixed layer is in use, so layer 1 has the properties
     if (use_temperature) then ; do j=js,je ; do i=is,ie
@@ -3876,7 +3882,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
       do k=1,nz
         call calculate_TFreeze(CS%tv%S(is:ie,j,k), pres(is:ie), T_freeze(is:ie), CS%tv%eqn_of_state)
         do i=is,ie
-          depth_ml = min(CS%HFrz, (US%Z_to_m*GV%m_to_H)*CS%visc%MLD(i,j))
+          depth_ml = min(CS%HFrz, CS%visc%h_ML(i,j))
           if (depth(i) + h(i,j,k) < depth_ml) then
             dh = h(i,j,k)
           elseif (depth(i) < depth_ml) then

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -164,7 +164,7 @@ module MOM_bulk_mixed_layer
 !> This subroutine partially steps the bulk mixed layer model.
 !! See \ref BML for more details.
 subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, CS, &
-                          optics, Hml, aggregate_FW_forcing, dt_diag, last_call)
+                          optics, BLD, H_ml, aggregate_FW_forcing, dt_diag, last_call)
   type(ocean_grid_type),      intent(inout) :: G      !< The ocean's grid structure.
   type(verticalGrid_type),    intent(in)    :: GV     !< The ocean's vertical grid structure.
   type(unit_scale_type),      intent(in)    :: US     !< A dimensional unit scaling type
@@ -195,7 +195,10 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
   type(optics_type),          pointer       :: optics !< The structure that can be queried for the
                                                       !! inverse of the vertical absorption decay
                                                       !! scale for penetrating shortwave radiation.
-  real, dimension(:,:),       pointer       :: Hml    !< Active mixed layer depth [Z ~> m]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                              intent(inout) :: BLD    !< Active mixed layer depth [Z ~> m]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                              intent(inout) :: H_ml   !< Active mixed layer thickness [H ~> m or kg m-2].
   logical,                    intent(in)    :: aggregate_FW_forcing !< If true, the net incoming and
                                                      !! outgoing surface freshwater fluxes are
                                                      !! combined before being applied, instead of
@@ -605,25 +608,27 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
       CS%ML_depth(i,j) = h(i,0)  ! Store the diagnostic.
     enddo ; endif
 
-    if (associated(Hml)) then
-      ! Return the mixed layerd depth in [Z ~> m].
-      if (GV%Boussinesq .or. GV%semi_Boussinesq) then
-        do i=is,ie
-          Hml(i,j) = G%mask2dT(i,j) * GV%H_to_Z*h(i,0)
-        enddo
+    ! Return the mixed layer depth in [Z ~> m].
+    if (GV%Boussinesq .or. GV%semi_Boussinesq) then
+      do i=is,ie
+        BLD(i,j) = G%mask2dT(i,j) * GV%H_to_Z*h(i,0)
+      enddo
+    else
+      do i=is,ie ; dp_ml(i) = GV%g_Earth * GV%H_to_RZ * h(i,0) ; enddo
+      if (associated(tv%p_surf)) then
+        do i=is,ie ; p_sfc(i) = tv%p_surf(i,j) ; enddo
       else
-        do i=is,ie ; dp_ml(i) = GV%g_Earth * GV%H_to_RZ * h(i,0) ; enddo
-        if (associated(tv%p_surf)) then
-          do i=is,ie ; p_sfc(i) = tv%p_surf(i,j) ; enddo
-        else
-          do i=is,ie ; p_sfc(i) = 0.0 ; enddo
-        endif
-        call average_specific_vol(T(:,0), S(:,0), p_sfc, dp_ml, SpV_ml, tv%eqn_of_state)
-        do i=is,ie
-          Hml(i,j) = G%mask2dT(i,j) * GV%H_to_RZ * SpV_ml(i) *  h(i,0)
-        enddo
+        do i=is,ie ; p_sfc(i) = 0.0 ; enddo
       endif
+      call average_specific_vol(T(:,0), S(:,0), p_sfc, dp_ml, SpV_ml, tv%eqn_of_state)
+      do i=is,ie
+        BLD(i,j) = G%mask2dT(i,j) * GV%H_to_RZ * SpV_ml(i) *  h(i,0)
+      enddo
     endif
+    ! Return the mixed layer thickness in [H ~> m or kg m-2].
+    do i=is,ie
+      H_ml(i,j) = G%mask2dT(i,j) * h(i,0)
+    enddo
 
 ! At this point, return water to the original layers, but constrained to
 ! still be sorted.  After this point, all the water that is in massive