diff --git a/vlaplex.f90 b/vlaplex.f90
index 44550cf..5031648 100644
--- a/vlaplex.f90
+++ b/vlaplex.f90
@@ -279,6 +279,9 @@ program VlaPlEx
     sum_ni = 0.0_dp
     ! Advect in the r direction
     do iz = 1, nz
+      if (all(n_i(iz,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
+        cycle
+      end if
       !$omp parallel do
       do i = 1, nr
         ! Advect negative velocity
@@ -383,6 +386,9 @@ program VlaPlEx
     f_i_old = f_i
 
     do iz = 1, nz
+      if (all(n_i(iz,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
+        cycle
+      end if
       ! Advect in the v direction
       ! i = 1, v<0
       i = 1
@@ -417,6 +423,9 @@ program VlaPlEx
     ! Reset values for next iteration
     f_i_old = f_i
     do iz = 1, nz
+      if (all(n_i(iz,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
+          cycle
+      end if
       fCum_i(iz,:)  = fCum_i(iz,:) + f_i_old(iz,rCum_index,:)
     end do
     ! Write output