diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90
index 5af0b774b0..b62f479354 100644
--- a/config_src/coupled_driver/ocean_model_MOM.F90
+++ b/config_src/coupled_driver/ocean_model_MOM.F90
@@ -442,8 +442,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
   real,    optional, intent(in)    :: cycle_length !< The duration of a coupled time stepping cycle [s].
 
   ! Local variables
-  type(time_type) :: Time_seg_start ! Stores the ocean model time at the start of this call to allow
-                            ! step_MOM to temporarily change the time as seen by internal modules.
+  type(time_type) :: Time_seg_start ! Stores the dynamic or thermodynamic ocean model time at the
+                            ! start of this call to allow step_MOM to temporarily change the time
+                            ! as seen by internal modules.
+  type(time_type) :: Time_thermo_start ! Stores the ocean model thermodynamics time at the start of
+                            ! this call to allow step_MOM to temporarily change the time as seen by
+                            ! internal modules.
   type(time_type) :: Time1  ! The value of the ocean model's time at the start of a call to step_MOM.
   integer :: index_bnds(4)  ! The computational domain index bounds in the ice-ocean boundary type.
   real :: weight            ! Flux accumulation weight of the current fluxes.
@@ -563,6 +567,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
     call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%restart_CSp)
   endif
 
+  Time_thermo_start = OS%Time
   Time_seg_start = OS%Time ; if (do_dyn) Time_seg_start = OS%Time_dyn
   Time1 = Time_seg_start
 
@@ -576,7 +581,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
                   reset_therm=Ocn_fluxes_used)
   elseif (OS%single_step_call) then
     call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves)
-  else
+  else  ! Step both the dynamics and thermodynamics with separate calls.
     n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001)
     dt_dyn = dt_coupling / real(n_max)
     thermo_does_span_coupling = (OS%thermo_spans_coupling .and. (OS%dt_therm > 1.5*dt_coupling))
@@ -636,7 +641,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
 
   if (do_dyn) OS%Time_dyn = Time_seg_start + Ocean_coupling_time_step
   if (do_dyn) OS%nstep = OS%nstep + 1
-  OS%Time = Time_seg_start  ! Reset the clock to compensate for shared pointers.
+  OS%Time = Time_thermo_start  ! Reset the clock to compensate for shared pointers.
   if (do_thermo) OS%Time = OS%Time + Ocean_coupling_time_step
   if (do_thermo) OS%nstep_thermo = OS%nstep_thermo + 1
 
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index 1b364d7ef0..21ec2e6dc6 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -734,8 +734,13 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, &
 
     !===========================================================================
     ! This is the second place where the diabatic processes and remapping could occur.
-    if (CS%t_dyn_rel_adv == 0.0 .and. do_thermo .and. .not.CS%diabatic_first) then
+    if ((CS%t_dyn_rel_adv==0.0) .and. do_thermo .and. (.not.CS%diabatic_first)) then
+
       dtdia = CS%t_dyn_rel_thermo
+      ! If the MOM6 dynamic and thermodynamic time stepping is being orchestrated
+      ! by the coupler, the value of diabatic_first does not matter.
+      if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) dtdia = dt
+
       if (CS%thermo_spans_coupling .and. (CS%dt_therm > 1.5*cycle_time) .and. &
           (abs(dt_therm - dtdia) > 1e-6*dt_therm)) then
         call MOM_error(FATAL, "step_MOM: Mismatch between dt_therm and dtdia "//&
@@ -750,7 +755,13 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, &
       call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, &
                            Time_local, .false., Waves=Waves)
       CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia
-      CS%t_dyn_rel_thermo = 0.0
+
+      if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) then
+        ! The diabatic processes are now ahead of the dynamics by dtdia.
+        CS%t_dyn_rel_thermo = -dtdia
+      else ! The diabatic processes and the dynamics are synchronized.
+        CS%t_dyn_rel_thermo = 0.0
+      endif
 
       if (dtdia > dt) & ! Reset CS%Time to its previous value.
         CS%Time = Time_start + real_to_time(rel_time - 0.5*dt)
diff --git a/src/framework/MOM_safe_alloc.F90 b/src/framework/MOM_safe_alloc.F90
index 75f5fda74e..47dd8376a3 100644
--- a/src/framework/MOM_safe_alloc.F90
+++ b/src/framework/MOM_safe_alloc.F90
@@ -10,13 +10,14 @@ module MOM_safe_alloc
 
 !> Allocate a pointer to a 1-d, 2-d or 3-d array
 interface safe_alloc_ptr
-  module procedure safe_alloc_ptr_3d_2arg, safe_alloc_ptr_2d_2arg
+  module procedure safe_alloc_ptr_3d_3arg,  safe_alloc_ptr_3d_6arg, safe_alloc_ptr_2d_2arg
   module procedure safe_alloc_ptr_3d, safe_alloc_ptr_2d, safe_alloc_ptr_1d
 end interface safe_alloc_ptr
 
 !> Allocate a 2-d or 3-d allocatable array
 interface safe_alloc_alloc
   module procedure safe_alloc_allocatable_3d, safe_alloc_allocatable_2d
+  module procedure safe_alloc_allocatable_3d_6arg
 end interface safe_alloc_alloc
 
 !   This combined interface might work with a later version of Fortran, but
@@ -57,7 +58,7 @@ subroutine safe_alloc_ptr_2d_2arg(ptr, ni, nj)
 end subroutine safe_alloc_ptr_2d_2arg
 
 !> Allocate a pointer to a 3-d array based on its dimension sizes
-subroutine safe_alloc_ptr_3d_2arg(ptr, ni, nj, nk)
+subroutine safe_alloc_ptr_3d_3arg(ptr, ni, nj, nk)
   real, dimension(:,:,:), pointer :: ptr !< A pointer to allocate
   integer, intent(in) :: ni !< The size of the 1st dimension of the array
   integer, intent(in) :: nj !< The size of the 2nd dimension of the array
@@ -66,7 +67,7 @@ subroutine safe_alloc_ptr_3d_2arg(ptr, ni, nj, nk)
     allocate(ptr(ni,nj,nk))
     ptr(:,:,:) = 0.0
   endif
-end subroutine safe_alloc_ptr_3d_2arg
+end subroutine safe_alloc_ptr_3d_3arg
 
 !> Allocate a pointer to a 2-d array based on its index starting and ending values
 subroutine safe_alloc_ptr_2d(ptr, is, ie, js, je)
@@ -95,6 +96,22 @@ subroutine safe_alloc_ptr_3d(ptr, is, ie, js, je, nk)
   endif
 end subroutine safe_alloc_ptr_3d
 
+!> Allocate a pointer to a 3-d array based on its index starting and ending values
+subroutine safe_alloc_ptr_3d_6arg(ptr, is, ie, js, je, ks, ke)
+  real, dimension(:,:,:), pointer :: ptr !< A pointer to allocate
+  integer, intent(in) :: is !< The start index to allocate for the 1st dimension
+  integer, intent(in) :: ie !< The end index to allocate for the 1st dimension
+  integer, intent(in) :: js !< The start index to allocate for the 2nd dimension
+  integer, intent(in) :: je !< The end index to allocate for the 2nd dimension
+  integer, intent(in) :: ks !< The start index to allocate for the 3rd dimension
+  integer, intent(in) :: ke !< The end index to allocate for the 3rd dimension
+  if (.not.associated(ptr)) then
+    allocate(ptr(is:ie,js:je,ks:ke))
+    ptr(:,:,:) = 0.0
+  endif
+end subroutine safe_alloc_ptr_3d_6arg
+
+
 !> Allocate a 2-d allocatable array based on its index starting and ending values
 subroutine safe_alloc_allocatable_2d(ptr, is, ie, js, je)
   real, dimension(:,:), allocatable :: ptr !< An allocatable array to allocate
@@ -109,6 +126,7 @@ subroutine safe_alloc_allocatable_2d(ptr, is, ie, js, je)
 end subroutine safe_alloc_allocatable_2d
 
 !> Allocate a 3-d allocatable array based on its index starting and ending values
+!! and k-index size
 subroutine safe_alloc_allocatable_3d(ptr, is, ie, js, je, nk)
   real, dimension(:,:,:), allocatable :: ptr !< An allocatable array to allocate
   integer, intent(in) :: is !< The start index to allocate for the 1st dimension
@@ -122,4 +140,19 @@ subroutine safe_alloc_allocatable_3d(ptr, is, ie, js, je, nk)
   endif
 end subroutine safe_alloc_allocatable_3d
 
+!> Allocate a 3-d allocatable array based on its 6 index starting and ending values
+subroutine safe_alloc_allocatable_3d_6arg(ptr, is, ie, js, je, ks, ke)
+  real, dimension(:,:,:), allocatable :: ptr !< An allocatable array to allocate
+  integer, intent(in) :: is !< The start index to allocate for the 1st dimension
+  integer, intent(in) :: ie !< The end index to allocate for the 1st dimension
+  integer, intent(in) :: js !< The start index to allocate for the 2nd dimension
+  integer, intent(in) :: je !< The end index to allocate for the 2nd dimension
+  integer, intent(in) :: ks !< The start index to allocate for the 3rd dimension
+  integer, intent(in) :: ke !< The end index to allocate for the 3rd dimension
+  if (.not.allocated(ptr)) then
+    allocate(ptr(is:ie,js:je,ks:ke))
+    ptr(:,:,:) = 0.0
+  endif
+end subroutine safe_alloc_allocatable_3d_6arg
+
 end module MOM_safe_alloc
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 200d3efdf7..24a529716d 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -422,6 +422,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
   eaml => eatr ; ebml => ebtr
 
   ! inverse time step
+  if (dt == 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
+        "diabatic was called with a zero length timestep.")
+  if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
+        "diabatic was called with a negative timestep.")
   Idt = 1.0 / dt
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
@@ -1301,6 +1305,10 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en
   eaml => eatr ; ebml => ebtr
 
   ! inverse time step
+  if (dt == 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
+        "legacy_diabatic was called with a zero length timestep.")
+  if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
+        "legacy_diabatic was called with a negative timestep.")
   Idt = 1.0 / dt
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
@@ -2506,7 +2514,7 @@ subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV,
   integer :: i, j, k, is, ie, js, je, nz
 
   is  = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
-  Idt = 1/dt
+  Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt
   work_3d(:,:,:) = 0.0
   work_2d(:,:)   = 0.0
 
@@ -2596,7 +2604,7 @@ subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old,
   integer :: i, j, k, is, ie, js, je, nz
 
   is  = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
-  Idt = 1/dt
+  Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt
   work_3d(:,:,:) = 0.0
   work_2d(:,:)   = 0.0
 
@@ -2683,7 +2691,7 @@ subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, CS)
   integer :: i, j, k, is, ie, js, je, nz
 
   is  = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
-  Idt = 1/dt
+  Idt = 0.0 ; if (dt > 0.0) Idt = 1. / dt
 
   ! temperature tendency
   if (CS%id_frazil_temp_tend > 0) then
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 5400f8c1df..e01374b5c6 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -669,7 +669,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
     allocate(CS%a1_shelf_v(G%isd:G%ied,G%JsdB:G%JedB)) ; CS%a1_shelf_v(:,:)=0.0
   endif
 
-  !$OMP parallel do default(private) shared(G,GV,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, &
+  !$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, &
   !$OMP                                     OBC,h_neglect,dt,I_valBL,Kv_u) &
   !$OMP                     firstprivate(i_hbbl)
   do j=G%Jsc,G%Jec
@@ -836,7 +836,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
 
 
   ! Now work on v-points.
-  !$OMP parallel do default(private) shared(G,GV,CS,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, &
+  !$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, &
   !$OMP                                  OBC,h_neglect,dt,I_valBL,Kv_v) &
   !$OMP                     firstprivate(i_hbbl)
   do J=Jsq,Jeq