diff --git a/src/util/util_add_to_layer.f90 b/src/util/util_add_to_layer.f90
index a99dd9dd..04ca287c 100644
--- a/src/util/util_add_to_layer.f90
+++ b/src/util/util_add_to_layer.f90
@@ -47,6 +47,7 @@ subroutine util_add_to_layer(user,surfi,crater)
    surfi%diam(layer) = crater%fcrat
    surfi%xl(layer) = crater%xl
    surfi%yl(layer) = crater%yl
+   surfi%timestamp(layer) = crater%timestamp
 else
    write(*,*)
    write(*,*) 'WARNING! No free layer to add crater pixel. Consider increasing NUMLAYERS'
diff --git a/src/util/util_remove_from_layer.f90 b/src/util/util_remove_from_layer.f90
index b9e773b7..22bb26d4 100644
--- a/src/util/util_remove_from_layer.f90
+++ b/src/util/util_remove_from_layer.f90
@@ -32,5 +32,6 @@ subroutine util_remove_from_layer(surfi,layer)
 surfi%diam(layer) = 0.0_DP
 surfi%xl(layer) = 0.0_SP
 surfi%yl(layer) = 0.0_SP
+surfi%timestamp(layer) = 0.0_SP
 
 end subroutine util_remove_from_layer
diff --git a/src/util/util_search.f90 b/src/util/util_search.f90
index caef13b8..61d3e415 100644
--- a/src/util/util_search.f90
+++ b/src/util/util_search.f90
@@ -92,6 +92,75 @@ subroutine util_search_double(arr,ind,n,val,klo)
 
 end subroutine util_search_double
 
+subroutine util_search_double_1(arr,ind,n,val,klo)
+use module_globals
+use module_util, EXCEPT_THIS_ONE => util_search_double_1
+implicit none
+
+   ! Arguments
+   integer(I4B),intent(in) :: ind,n
+   real(DP),dimension(:),intent(in) :: arr
+   real(DP),intent(in) :: val
+   integer(I4B),intent(inout) :: klo
+
+   ! Internals
+   integer(I4B) :: k,khi,inc
+   logical :: uporder ! true if arr is in ascending order
+
+   ! Executable code
+   uporder = arr(1) <= arr(n) 
+   ! First hunt for the starting value
+   if (.not.((klo <= 0).or.(klo > n))) then  ! Only proceed if guess is useful
+      inc = 1
+      if ((val >= arr(klo)) .eqv. uporder) then
+         do 
+            khi = klo + inc 
+            if (khi > n) then
+               khi = n + 1
+            else if ((val >= arr(khi)) .eqv. uporder) then
+               klo = khi
+               inc = 2 * inc
+               cycle
+            end if
+            exit
+         end do
+      else
+         khi = klo
+         do
+            klo = khi - inc
+            if (klo < 1 ) then
+               klo = 0
+            else if ((val < arr(klo)) .eqv. uporder) then
+               khi = klo
+               inc = 2 * inc
+               cycle
+            end if
+            exit
+         end do
+       end if
+   else
+      klo = 0 
+      khi = n + 1
+   end if
+
+   ! Start the bisection method search
+   do
+      if (khi - klo == 1) then
+         if (val == arr(n)) klo = n-1
+         if (val == arr(1)) klo = 1
+         return
+      end if
+      k = (khi + klo) / 2
+      if ((val >= arr(k)) .eqv. uporder) then
+         klo = k
+      else
+         khi = k
+      end if
+   end do
+
+end subroutine util_search_double_1
+
+
 subroutine util_search_int(arr,ind,n,val,klo)
 use module_globals
 use module_util, EXCEPT_THIS_ONE => util_search_int
diff --git a/src/util/util_sort_layer.f90 b/src/util/util_sort_layer.f90
index a230c4c5..c3d055b1 100644
--- a/src/util/util_sort_layer.f90
+++ b/src/util/util_sort_layer.f90
@@ -31,7 +31,7 @@ subroutine util_sort_layer(user,surf,crater)
    ! Internal variables
    integer(I4B),dimension(user%numlayers) :: ind
    real(DP),dimension(user%numlayers) :: tempdiam
-   real(SP),dimension(user%numlayers) :: tempxpos,tempypos
+   real(SP),dimension(user%numlayers) :: tempxpos,tempypos,temptime
    integer(I4B) :: i,j,k,inc,incsq,mx,my,iradsq
 
    inc = min(crater%rimdispx,(user%gridsize - 1)/2)
@@ -54,6 +54,7 @@ subroutine util_sort_layer(user,surf,crater)
             tempdiam = surf(mx,my)%diam(1:user%numlayers)
             tempxpos = surf(mx,my)%xl(1:user%numlayers)
             tempypos = surf(mx,my)%yl(1:user%numlayers)
+            temptime = surf(mx,my)%timestamp(1:user%numlayers)
 
             ! Sort the layers by crater diameter
             call util_mrgrnk(surf(mx,my)%diam(1:user%numlayers),ind)
@@ -62,6 +63,7 @@ subroutine util_sort_layer(user,surf,crater)
                surf(mx,my)%diam(k)=tempdiam(ind(k))
                surf(mx,my)%xl(k)=tempxpos(ind(k))
                surf(mx,my)%yl(k)=tempypos(ind(k))
+               surf(mx,my)%timestamp(k)=temptime(ind(k))
             end do
          end if
       end do
diff --git a/test/CTEM b/test/CTEM
deleted file mode 100755
index 611133f6..00000000
Binary files a/test/CTEM and /dev/null differ