diff --git a/src/SIS_dyn_trans.F90 b/src/SIS_dyn_trans.F90
index 65fe177f..e6605cb5 100644
--- a/src/SIS_dyn_trans.F90
+++ b/src/SIS_dyn_trans.F90
@@ -350,7 +350,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
   if ((CS%dt_ice_dyn > 0.0) .and. (CS%dt_ice_dyn < dt_slow)) &
     ndyn_steps = nadv_cycle * max(CEILING(dt_slow/(nadv_cycle*CS%dt_ice_dyn) - 0.000001), 1)
   dt_slow_dyn = dt_slow / ndyn_steps
-  if (CS%adv_substeps > 0) dt_adv = dt_slow_dyn / real(CS%adv_substeps)
+  dt_adv = dt_slow_dyn / real(CS%adv_substeps)
   nts_last = (ndyn_steps/nadv_cycle)*CS%adv_substeps
   if (CS%merged_cont .and. (CS%nts == 0) .and. (nts_last > CS%max_nts)) &
     call increase_max_tracer_step_memory(CS, G, nts_last)
@@ -675,13 +675,13 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G,
   real :: dt_slow_dyn  ! The slow dynamics timestep [s].
   real :: dt_adv       ! The advective timestep [s].
   real :: Idt_slow     ! The inverse of dt_slow [s-1].
+  type(time_type) :: Time_cycle_start ! The model's time at the start of an advective cycle.
   real, dimension(SZI_(G),SZJ_(G)) :: &
     rdg_rate  ! A ridging rate [s-1], this will be calculated from the strain rates in the dynamics.
   integer :: i, j, k, n, isc, iec, jsc, jec, ncat
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
   integer :: ndyn_steps, nds ! The number of dynamic steps in this call.
   integer :: nadv_cycle, nac ! The number of tracer advective cycles in this call.
-  integer :: nts_last ! The number of tracer advection steps before updating IST.
   character(len=256) :: mesg
 
   real, parameter :: T_0degC = 273.15 ! 0 degrees C in Kelvin
@@ -697,24 +697,21 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G,
   ndyn_steps = 1 ; nadv_cycle = 1
   if (CS%merged_cont .and. (CS%dt_advect > 0.0) .and. (CS%dt_advect < dt_slow)) &
     nadv_cycle = max(CEILING(dt_slow/CS%dt_advect - 1e-9), 1)
-  if ((CS%dt_ice_dyn > 0.0) .and. (CS%dt_ice_dyn < dt_slow)) &
+  if ((CS%dt_ice_dyn > 0.0) .and. (CS%dt_ice_dyn*nadv_cycle < dt_slow)) &
     ndyn_steps = max(CEILING(dt_slow/(nadv_cycle*CS%dt_ice_dyn) - 0.000001), 1)
   dt_slow_dyn = dt_slow / (nadv_cycle * ndyn_steps)
-  if (CS%adv_substeps > 0) dt_adv = dt_slow_dyn / real(CS%adv_substeps)
-  nts_last = ndyn_steps*CS%adv_substeps
-  if (CS%merged_cont .and. (CS%nts == 0) .and. (nts_last > CS%max_nts)) &
-    call increase_max_tracer_step_memory(CS, G, nts_last)
+  dt_adv = dt_slow_dyn / real(CS%adv_substeps)
 
   complete_ice_cover = 1.0 - 2.0*ncat*epsilon(complete_ice_cover)
 
-  do nac=0,nadv_cycle-1 ; do nds=1,ndyn_steps
-
-    call enable_SIS_averaging(dt_slow_dyn, CS%Time - real_to_time((ndyn_steps*(nadv_cycle-nac)-nds)*dt_slow_dyn), CS%diag)
-
+  do nac=1,nadv_cycle
     ! Convert the category-resolved ice state into the simplified 2-d ice state.
     ! This should be called after a thermodynamic step or if ice_transport was called.
     if (CS%nts == 0) then
       call mpp_clock_begin(iceClock4)
+      if (ndyn_steps*CS%adv_substeps > CS%max_nts) &
+        call increase_max_tracer_step_memory(CS, G, ndyn_steps*CS%adv_substeps)
+
       CS%mca_step(:,:,0) = 0.0 ; CS%mi_sum(:,:) = 0.0 ; CS%ice_cover(:,:) = 0.0
       !$OMP parallel do default(shared)
       do j=jsd,jed ; do k=1,ncat ; do i=isd,ied
@@ -745,188 +742,190 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G,
       call mpp_clock_end(iceClock4)
     endif
 
-    ! This is the end of the code that might be called as 2-d dynamics.
-    !   Variables updated here: CS%ice_cover, CS%[uv]_ice_[BC], CS%mca_step, CS%mi_sum,
-    !      IOF (stresses), CS%[uv]h_step, CS%nts, CS%SIS_[BC]_dyn_CSp
-    !   Variables used with intent in here: FIA, G, OSS, dt_slow_dyn, CS
-    !   Local variables: ice_free
-
-
-    call mpp_clock_begin(iceClock4)
-    do j=jsd,jed ; do i=isd,ied ; ice_free(i,j) = max(1.0 - CS%ice_cover(i,j), 0.0) ; enddo ; enddo
-
-    ! In the dynamics code, only the ice velocities are changed, and the ice-ocean
-    ! stresses are calculated.  The gravity wave dynamics (i.e. the continuity
-    ! equation) are not included in the dynamics (yet).  All of the thickness categories
-    ! are merged together.
-    if (CS%Cgrid_dyn) then
-
-      ! Correct the wind stresses for changes in the fractional ice-coverage and set
-      ! the wind stresses on the ice and the open ocean for a C-grid staggering.
-      ! This block of code must be executed if ice_cover and ice_free or the various wind
-      ! stresses were updated.
-      call set_wind_stresses_C(FIA, CS%ice_cover, ice_free, WindStr_x_Cu, WindStr_y_Cv, &
-                               WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, complete_ice_cover)
-
-      if (CS%debug) then
-        call uvchksum("Before SIS_C_dynamics [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
-        call hchksum(ice_free, "ice_free before SIS_C_dynamics", G%HI)
-        call hchksum(CS%mca_step(:,:,CS%nts), "misp_sum before SIS_C_dynamics", G%HI)
-        call hchksum(CS%mi_sum, "mi_sum before SIS_C_dynamics", G%HI)
-        call hchksum(OSS%sea_lev, "sea_lev before SIS_C_dynamics", G%HI, haloshift=1)
-        call hchksum(CS%ice_cover, "ice_cover before SIS_C_dynamics", G%HI, haloshift=1)
-        call uvchksum("[uv]_ocn before SIS_C_dynamics", OSS%u_ocn_C, OSS%v_ocn_C, G, halos=1)
-        call uvchksum("WindStr_[xy] before SIS_C_dynamics", WindStr_x_Cu, WindStr_y_Cv, G, halos=1)
-!        call hchksum_pair("WindStr_[xy]_A before SIS_C_dynamics", WindStr_x_A, WindStr_y_A, G, halos=1)
-      endif
+    Time_cycle_start = CS%Time - real_to_time((nadv_cycle-(nac-1))*ndyn_steps*dt_slow_dyn)
 
-      !### Ridging needs to be added with C-grid dynamics.
-      rdg_rate(:,:) = 0.0
-      call mpp_clock_begin(iceClocka)
-      call SIS_C_dynamics(CS%ice_cover, CS%mca_step(:,:,CS%nts), CS%mi_sum, CS%u_ice_C, CS%v_ice_C, &
-                          OSS%u_ocn_C, OSS%v_ocn_C, WindStr_x_Cu, WindStr_y_Cv, OSS%sea_lev, &
-                          str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, CS%SIS_C_dyn_CSp)
-      call mpp_clock_end(iceClocka)
+    do nds=1,ndyn_steps
+      ! This is the start of the code that might be called as 2-d dynamics.
+      !   Variables updated here: CS%ice_cover, CS%[uv]_ice_[BC], CS%mca_step, CS%mi_sum,
+      !      IOF (stresses), CS%[uv]h_step, CS%nts, CS%SIS_[BC]_dyn_CSp
+      !   Variables used with intent in here: FIA, G, OSS, dt_slow_dyn, CS
+      !   Local variables: ice_free
 
-      if (CS%debug) call uvchksum("After ice_dynamics [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
+      call mpp_clock_begin(iceClock4)
+      call enable_SIS_averaging(dt_slow_dyn, Time_cycle_start + real_to_time(nds*dt_slow_dyn), CS%diag)
+      do j=jsd,jed ; do i=isd,ied ; ice_free(i,j) = max(1.0 - CS%ice_cover(i,j), 0.0) ; enddo ; enddo
 
-      call mpp_clock_begin(iceClockb)
-      call pass_vector(CS%u_ice_C, CS%v_ice_C, G%Domain, stagger=CGRID_NE)
-      call pass_vector(str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, G%Domain, stagger=CGRID_NE)
-      call mpp_clock_end(iceClockb)
+      ! In the dynamics code, only the ice velocities are changed, and the ice-ocean
+      ! stresses are calculated.  The gravity wave dynamics (i.e. the continuity
+      ! equation) are not included in the dynamics (yet).  All of the thickness categories
+      ! are merged together.
+      if (CS%Cgrid_dyn) then
 
-      ! Dynamics diagnostics
-      call mpp_clock_begin(iceClockc)
-      if (CS%id_fax>0) call post_data(CS%id_fax, WindStr_x_Cu, CS%diag)
-      if (CS%id_fay>0) call post_data(CS%id_fay, WindStr_y_Cv, CS%diag)
-      if (CS%debug) call uvchksum("Before set_ocean_top_stress_Cgrid [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
+        ! Correct the wind stresses for changes in the fractional ice-coverage and set
+        ! the wind stresses on the ice and the open ocean for a C-grid staggering.
+        ! This block of code must be executed if ice_cover and ice_free or the various wind
+        ! stresses were updated.
+        call set_wind_stresses_C(FIA, CS%ice_cover, ice_free, WindStr_x_Cu, WindStr_y_Cv, &
+                                 WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, complete_ice_cover)
+
+        if (CS%debug) then
+          call uvchksum("Before SIS_C_dynamics [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
+          call hchksum(ice_free, "ice_free before SIS_C_dynamics", G%HI)
+          call hchksum(CS%mca_step(:,:,CS%nts), "misp_sum before SIS_C_dynamics", G%HI)
+          call hchksum(CS%mi_sum, "mi_sum before SIS_C_dynamics", G%HI)
+          call hchksum(OSS%sea_lev, "sea_lev before SIS_C_dynamics", G%HI, haloshift=1)
+          call hchksum(CS%ice_cover, "ice_cover before SIS_C_dynamics", G%HI, haloshift=1)
+          call uvchksum("[uv]_ocn before SIS_C_dynamics", OSS%u_ocn_C, OSS%v_ocn_C, G, halos=1)
+          call uvchksum("WindStr_[xy] before SIS_C_dynamics", WindStr_x_Cu, WindStr_y_Cv, G, halos=1)
+  !        call hchksum_pair("WindStr_[xy]_A before SIS_C_dynamics", WindStr_x_A, WindStr_y_A, G, halos=1)
+        endif
 
-      ! Store all mechanical ocean forcing.
-      call set_ocean_top_stress_C2(IOF, WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, &
-                                   str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, ice_free, CS%ice_cover, G)
-      call mpp_clock_end(iceClockc)
+        !### Ridging needs to be added with C-grid dynamics.
+        rdg_rate(:,:) = 0.0
+        call mpp_clock_begin(iceClocka)
+        call SIS_C_dynamics(CS%ice_cover, CS%mca_step(:,:,CS%nts), CS%mi_sum, CS%u_ice_C, CS%v_ice_C, &
+                            OSS%u_ocn_C, OSS%v_ocn_C, WindStr_x_Cu, WindStr_y_Cv, OSS%sea_lev, &
+                            str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, CS%SIS_C_dyn_CSp)
+        call mpp_clock_end(iceClocka)
 
-    else ! B-grid dynamics.
+        if (CS%debug) call uvchksum("After ice_dynamics [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
 
-      ! Correct the wind stresses for changes in the fractional ice-coverage and set
-      ! the wind stresses on the ice and the open ocean for a C-grid staggering.
-      ! This block of code must be executed if ice_cover and ice_free or the various wind
-      ! stresses were updated.
+        call mpp_clock_begin(iceClockb)
+        call pass_vector(CS%u_ice_C, CS%v_ice_C, G%Domain, stagger=CGRID_NE)
+        call pass_vector(str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, G%Domain, stagger=CGRID_NE)
+        call mpp_clock_end(iceClockb)
 
-      call set_wind_stresses_B(FIA, CS%ice_cover, ice_free, WindStr_x_B, WindStr_y_B, &
-                               WindStr_x_ocn_B, WindStr_y_ocn_B, G, complete_ice_cover)
+        ! Dynamics diagnostics
+        call mpp_clock_begin(iceClockc)
+        if (CS%id_fax>0) call post_data(CS%id_fax, WindStr_x_Cu, CS%diag)
+        if (CS%id_fay>0) call post_data(CS%id_fay, WindStr_y_Cv, CS%diag)
+        if (CS%debug) call uvchksum("Before set_ocean_top_stress_Cgrid [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
 
-      if (CS%debug) then
-        call Bchksum_pair("[uv]_ice_B before dynamics", CS%u_ice_B, CS%v_ice_B, G)
-        call hchksum(ice_free, "ice_free before ice_dynamics", G%HI)
-        call hchksum(CS%mca_step(:,:,CS%nts), "misp_sum before ice_dynamics", G%HI)
-        call hchksum(CS%mi_sum, "mi_sum before ice_dynamics", G%HI)
-        call hchksum(OSS%sea_lev, "sea_lev before ice_dynamics", G%HI, haloshift=1)
-        call Bchksum_pair("[uv]_ocn before ice_dynamics", OSS%u_ocn_B, OSS%v_ocn_B, G)
-        call Bchksum_pair("WindStr_[xy]_B before ice_dynamics", WindStr_x_B, WindStr_y_B, G, halos=1)
-      endif
+        ! Store all mechanical ocean forcing.
+        call set_ocean_top_stress_C2(IOF, WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, &
+                                     str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, ice_free, CS%ice_cover, G)
+        call mpp_clock_end(iceClockc)
+
+      else ! B-grid dynamics.
+
+        ! Correct the wind stresses for changes in the fractional ice-coverage and set
+        ! the wind stresses on the ice and the open ocean for a C-grid staggering.
+        ! This block of code must be executed if ice_cover and ice_free or the various wind
+        ! stresses were updated.
+
+        call set_wind_stresses_B(FIA, CS%ice_cover, ice_free, WindStr_x_B, WindStr_y_B, &
+                                 WindStr_x_ocn_B, WindStr_y_ocn_B, G, complete_ice_cover)
+
+        if (CS%debug) then
+          call Bchksum_pair("[uv]_ice_B before dynamics", CS%u_ice_B, CS%v_ice_B, G)
+          call hchksum(ice_free, "ice_free before ice_dynamics", G%HI)
+          call hchksum(CS%mca_step(:,:,CS%nts), "misp_sum before ice_dynamics", G%HI)
+          call hchksum(CS%mi_sum, "mi_sum before ice_dynamics", G%HI)
+          call hchksum(OSS%sea_lev, "sea_lev before ice_dynamics", G%HI, haloshift=1)
+          call Bchksum_pair("[uv]_ocn before ice_dynamics", OSS%u_ocn_B, OSS%v_ocn_B, G)
+          call Bchksum_pair("WindStr_[xy]_B before ice_dynamics", WindStr_x_B, WindStr_y_B, G, halos=1)
+        endif
 
-      rdg_rate(:,:) = 0.0
-      call mpp_clock_begin(iceClocka)
-      call SIS_B_dynamics(CS%ice_cover, CS%mca_step(:,:,CS%nts), CS%mi_sum, CS%u_ice_B, CS%v_ice_B, &
-                          OSS%u_ocn_B, OSS%v_ocn_B, WindStr_x_B, WindStr_y_B, OSS%sea_lev, &
-                          str_x_ice_ocn_B, str_y_ice_ocn_B, CS%do_ridging, &
-                          rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, CS%SIS_B_dyn_CSp)
-      call mpp_clock_end(iceClocka)
+        rdg_rate(:,:) = 0.0
+        call mpp_clock_begin(iceClocka)
+        call SIS_B_dynamics(CS%ice_cover, CS%mca_step(:,:,CS%nts), CS%mi_sum, CS%u_ice_B, CS%v_ice_B, &
+                            OSS%u_ocn_B, OSS%v_ocn_B, WindStr_x_B, WindStr_y_B, OSS%sea_lev, &
+                            str_x_ice_ocn_B, str_y_ice_ocn_B, CS%do_ridging, &
+                            rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, CS%SIS_B_dyn_CSp)
+        call mpp_clock_end(iceClocka)
+
+        if (CS%debug) call Bchksum_pair("After dynamics [uv]_ice_B", CS%u_ice_B, CS%v_ice_B, G)
+
+        call mpp_clock_begin(iceClockb)
+        call pass_vector(CS%u_ice_B, CS%v_ice_B, G%Domain, stagger=BGRID_NE)
+        call mpp_clock_end(iceClockb)
+
+        ! Dynamics diagnostics
+        call mpp_clock_begin(iceClockc)
+        if ((CS%id_fax>0) .or. (CS%id_fay>0)) then
+          !$OMP parallel do default(shared) private(ps_vel)
+          do J=jsc-1,jec ; do I=isc-1,iec
+            ps_vel = (1.0 - G%mask2dBu(I,J)) + 0.25*G%mask2dBu(I,J) * &
+                  ((ice_free(i+1,j+1) + ice_free(i,j)) + &
+                   (ice_free(i+1,j) + ice_free(i,j+1)) )
+            diagVarBx(I,J) = ps_vel * WindStr_x_ocn_B(I,J) + (1.0-ps_vel) * WindStr_x_B(I,J)
+            diagVarBy(I,J) = ps_vel * WindStr_y_ocn_B(I,J) + (1.0-ps_vel) * WindStr_y_B(I,J)
+          enddo ; enddo
 
-      if (CS%debug) call Bchksum_pair("After dynamics [uv]_ice_B", CS%u_ice_B, CS%v_ice_B, G)
+          if (CS%id_fax>0) call post_data(CS%id_fax, diagVarBx, CS%diag)
+          if (CS%id_fay>0) call post_data(CS%id_fay, diagVarBy, CS%diag)
+        endif
 
-      call mpp_clock_begin(iceClockb)
-      call pass_vector(CS%u_ice_B, CS%v_ice_B, G%Domain, stagger=BGRID_NE)
-      call mpp_clock_end(iceClockb)
+        if (CS%debug) call Bchksum_pair("Before set_ocean_top_stress_Bgrid [uv]_ice_B", CS%u_ice_B, CS%v_ice_B, G)
+        ! Store all mechanical ocean forcing.
+        call set_ocean_top_stress_B2(IOF, WindStr_x_ocn_B, WindStr_y_ocn_B, &
+                                     str_x_ice_ocn_B, str_y_ice_ocn_B, ice_free, CS%ice_cover, G)
+        call mpp_clock_end(iceClockc)
 
-      ! Dynamics diagnostics
-      call mpp_clock_begin(iceClockc)
-      if ((CS%id_fax>0) .or. (CS%id_fay>0)) then
-        !$OMP parallel do default(shared) private(ps_vel)
-        do J=jsc-1,jec ; do I=isc-1,iec
-          ps_vel = (1.0 - G%mask2dBu(I,J)) + 0.25*G%mask2dBu(I,J) * &
-                ((ice_free(i+1,j+1) + ice_free(i,j)) + &
-                 (ice_free(i+1,j) + ice_free(i,j+1)) )
-          diagVarBx(I,J) = ps_vel * WindStr_x_ocn_B(I,J) + (1.0-ps_vel) * WindStr_x_B(I,J)
-          diagVarBy(I,J) = ps_vel * WindStr_y_ocn_B(I,J) + (1.0-ps_vel) * WindStr_y_B(I,J)
+        ! Convert the velocities to C-grid points for use in transport.
+        do j=jsc,jec ; do I=isc-1,iec
+          CS%u_ice_C(I,j) = 0.5 * ( CS%u_ice_B(I,J-1) + CS%u_ice_B(I,J) )
         enddo ; enddo
+        do J=jsc-1,jec ; do i=isc,iec
+          CS%v_ice_C(i,J) = 0.5 * ( CS%v_ice_B(I-1,J) + CS%v_ice_B(I,J) )
+        enddo ; enddo
+      endif ! End of B-grid dynamics
 
-        if (CS%id_fax>0) call post_data(CS%id_fax, diagVarBx, CS%diag)
-        if (CS%id_fay>0) call post_data(CS%id_fay, diagVarBy, CS%diag)
-      endif
-
-      if (CS%debug) call Bchksum_pair("Before set_ocean_top_stress_Bgrid [uv]_ice_B", CS%u_ice_B, CS%v_ice_B, G)
-      ! Store all mechanical ocean forcing.
-      call set_ocean_top_stress_B2(IOF, WindStr_x_ocn_B, WindStr_y_ocn_B, &
-                                   str_x_ice_ocn_B, str_y_ice_ocn_B, ice_free, CS%ice_cover, G)
-      call mpp_clock_end(iceClockc)
-
-      ! Convert the velocities to C-grid points for use in transport.
-      do j=jsc,jec ; do I=isc-1,iec
-        CS%u_ice_C(I,j) = 0.5 * ( CS%u_ice_B(I,J-1) + CS%u_ice_B(I,J) )
-      enddo ; enddo
-      do J=jsc-1,jec ; do i=isc,iec
-        CS%v_ice_C(i,J) = 0.5 * ( CS%v_ice_B(I-1,J) + CS%v_ice_B(I,J) )
-      enddo ; enddo
-    endif ! End of B-grid dynamics
+      if (CS%debug) call uvchksum("Before ice_transport [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
+      call enable_SIS_averaging(dt_slow_dyn, Time_cycle_start + real_to_time(nds*dt_slow_dyn), CS%diag)
 
-    if (CS%debug) call uvchksum("Before ice_transport [uv]_ice_C", CS%u_ice_C, CS%v_ice_C, G)
-    call enable_SIS_averaging(dt_slow_dyn, CS%Time - real_to_time((ndyn_steps*(nadv_cycle-nac)-nds)*dt_slow_dyn), CS%diag)
-
-    ! Update the integrated ice mass and store the transports in each step.
-    if (CS%nts+CS%adv_substeps > CS%max_nts) call SIS_error(FATAL, &
-      "Attempting to store more advective substeps than allocated space allows.  Increase MAX_NTS.")
-    do n = CS%nts+1, CS%nts+CS%adv_substeps
-      if (n < nts_last) then
-        ! Some of the work is not needed for the last step before cat_ice_transport.
-        call summed_continuity(CS%u_ice_C, CS%v_ice_C, CS%mca_step(:,:,n-1), CS%mca_step(:,:,n), &
-                               CS%uh_step(:,:,n), CS%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp, &
-                               h_ice=CS%mi_sum)
-        call ice_cover_transport(CS%u_ice_C, CS%v_ice_C, CS%ice_cover, dt_adv, G, IG, CS%cover_trans_CSp)
-        call pass_var(CS%mi_sum, G%Domain, complete=.false.)
-        call pass_var(CS%ice_cover, G%Domain, complete=.false.)
-        call pass_var(CS%mca_step(:,:,n), G%Domain, complete=.true.)
-      else
-        call summed_continuity(CS%u_ice_C, CS%v_ice_C, CS%mca_step(:,:,n-1), CS%mca_step(:,:,n), &
-                               CS%uh_step(:,:,n), CS%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp)
-      endif
-    enddo
-    CS%nts = CS%nts + CS%adv_substeps
-    call mpp_clock_end(iceClock4)
+      ! Update the integrated ice mass and store the transports in each step.
+      if (CS%nts+CS%adv_substeps > CS%max_nts) call SIS_error(FATAL, &
+        "Attempting to store more advective substeps than allocated space allows.  Increase MAX_NTS.")
+      do n = CS%nts+1, CS%nts+CS%adv_substeps
+        if (n < ndyn_steps*CS%adv_substeps) then
+          ! Some of the work is not needed for the last step before cat_ice_transport.
+          call summed_continuity(CS%u_ice_C, CS%v_ice_C, CS%mca_step(:,:,n-1), CS%mca_step(:,:,n), &
+                                 CS%uh_step(:,:,n), CS%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp, &
+                                 h_ice=CS%mi_sum)
+          call ice_cover_transport(CS%u_ice_C, CS%v_ice_C, CS%ice_cover, dt_adv, G, IG, CS%cover_trans_CSp)
+          call pass_var(CS%mi_sum, G%Domain, complete=.false.)
+          call pass_var(CS%ice_cover, G%Domain, complete=.false.)
+          call pass_var(CS%mca_step(:,:,n), G%Domain, complete=.true.)
+        else
+          call summed_continuity(CS%u_ice_C, CS%v_ice_C, CS%mca_step(:,:,n-1), CS%mca_step(:,:,n), &
+                                 CS%uh_step(:,:,n), CS%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp)
+        endif
+      enddo
+      CS%nts = CS%nts + CS%adv_substeps
+      call mpp_clock_end(iceClock4)
 
+    enddo ! nds=1,ndyn_steps
     ! This is the end of the code that might be called as 2-d dynamics.
 
     call mpp_clock_begin(iceClock8)
-    if (nds == ndyn_steps) then
-      ! Do the transport of mass and tracers by category and vertical layer.
-      call ice_cat_transport(CS%CAS, IST%TrReg, dt_slow_dyn, CS%nts, G, IG, &
-                             CS%SIS_transport_CSp, mca_tot=CS%mca_step(:,:,0:CS%nts), &
-                             uh_tot=CS%uh_step(:,:,1:CS%nts), vh_tot=CS%vh_step(:,:,1:CS%nts))
-      ! Convert the cell-averaged state back to the ice-state type, adjusting the
-      ! category mass distributions, doing ridging, and updating the partition sizes.
-      call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp, &
-                                rdg_rate=rdg_rate)
-
-      CS%nts = 0 ! There is no outstanding transport to be done.
-
-      ! Copy the velocities back to the ice state type
-      if (CS%Cgrid_dyn) then
-        do j=jsd,jed ; do I=IsdB,IedB ; IST%u_ice_C(I,j) = CS%u_ice_C(I,j) ; enddo ; enddo
-        do J=JsdB,JedB ; do i=isd,ied ; IST%v_ice_C(i,J) = CS%v_ice_C(i,J) ; enddo ; enddo
-      else
-        do J=JsdB,JedB ; do I=IsdB,IedB
-          IST%u_ice_B(I,J) = CS%u_ice_B(I,J) ; IST%v_ice_B(I,J) = CS%v_ice_B(I,J)
-        enddo ; enddo
-      endif
-      if (CS%column_check) &
-        call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
-                                message="      Post_transport")! , check_column=.true.)
+    ! Do the transport of mass and tracers by category and vertical layer.
+    call ice_cat_transport(CS%CAS, IST%TrReg, dt_slow_dyn, CS%nts, G, IG, &
+                           CS%SIS_transport_CSp, mca_tot=CS%mca_step(:,:,0:CS%nts), &
+                           uh_tot=CS%uh_step(:,:,1:CS%nts), vh_tot=CS%vh_step(:,:,1:CS%nts))
+    ! Convert the cell-averaged state back to the ice-state type, adjusting the
+    ! category mass distributions, doing ridging, and updating the partition sizes.
+    call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp, &
+                              rdg_rate=rdg_rate)
+
+    CS%nts = 0 ! There is no outstanding transport to be done.
+
+    ! Copy the velocities back to the ice state type
+    if (CS%Cgrid_dyn) then
+      do j=jsd,jed ; do I=IsdB,IedB ; IST%u_ice_C(I,j) = CS%u_ice_C(I,j) ; enddo ; enddo
+      do J=JsdB,JedB ; do i=isd,ied ; IST%v_ice_C(i,J) = CS%v_ice_C(i,J) ; enddo ; enddo
+    else
+      do J=JsdB,JedB ; do I=IsdB,IedB
+        IST%u_ice_B(I,J) = CS%u_ice_B(I,J) ; IST%v_ice_B(I,J) = CS%v_ice_B(I,J)
+      enddo ; enddo
     endif
+    if (CS%column_check) &
+      call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
+                              message="      Post_transport")! , check_column=.true.)
 
     call mpp_clock_end(iceClock8)
 
-  enddo ; enddo ! nds=1,ndyn_steps*nadv_cycle
+  enddo ! nac=0,nadv_cycle-1
   call finish_ocean_top_stresses(IOF, G)
 
   ! Set appropriate surface quantities in categories with no ice.
@@ -2446,6 +2445,7 @@ subroutine SIS_dyn_trans_init(Time, G, IG, param_file, diag, CS, output_dir, Tim
     call get_param(param_file, mdl, "NSTEPS_ADV", CS%adv_substeps, &
                  "The number of advective iterations for each slow dynamics \n"//&
                  "time step.", default=1)
+    if (CS%adv_substeps < 1) CS%adv_substeps = 1
 
     call get_param(param_file, mdl, "ICEBERG_WINDSTRESS_BUG", CS%berg_windstress_bug, &
                  "If true, use older code that applied an old ice-ocean \n"//&
@@ -2511,9 +2511,9 @@ subroutine SIS_dyn_trans_init(Time, G, IG, param_file, diag, CS, output_dir, Tim
       CS%nts = 0 ; CS%max_nts = 0
       call safe_alloc_alloc(CS%mi_sum, G%isd, G%ied, G%jsd, G%jed)
       call safe_alloc_alloc(CS%ice_cover, G%isd, G%ied, G%jsd, G%jed)
-      max_nts = max(CS%adv_substeps, 1)
+      max_nts = CS%adv_substeps
       if ((CS%dt_ice_dyn > 0.0) .and. (CS%dt_advect > CS%dt_ice_dyn)) &
-        max_nts = max(CS%adv_substeps, 1) * max(CEILING(CS%dt_advect/CS%dt_ice_dyn - 1e-6), 1)
+        max_nts = CS%adv_substeps * max(CEILING(CS%dt_advect/CS%dt_ice_dyn - 1e-6), 1)
       call increase_max_tracer_step_memory(CS, G, max_nts)
 
       call safe_alloc_alloc(CS%u_ice_C, G%IsdB, G%IedB, G%jsd, G%jed)