diff --git a/2d/b4step2.f90 b/2d/b4step2.f90
index e0e07b4..0b33930 100644
--- a/2d/b4step2.f90
+++ b/2d/b4step2.f90
@@ -3,7 +3,7 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux)
 
 !   set slip before call to step2
 
-    use fault_module, only: center, xcb, nsubfaults, subfaults, final_rupture_rise_time, LAT2METER
+    use fault_module, only: center, xcb, nsubfaults, subfaults, final_rupture_rise_time
 
     implicit none
     integer, intent(in) :: mbc,mx,my,meqn,maux
@@ -11,7 +11,7 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux)
     real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
     real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)
 
-    integer :: i, j, k 
+    integer :: i, j, k
     real(kind=8) :: xcell, ycell
 
     aux(13,:,:) = 0.d0
@@ -26,8 +26,8 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux)
             if (xcb(1)-1.d-10 <= xcell - 0.5d0*dx .and. xcell + 0.5d0*dx <= xcb(2)+1.d-10) then
 
               do k=1,nsubfaults
-                if (subfaults(k)%longitude*LAT2METER <= xcell .and. &
-                    xcell <= subfaults(k)%longitude*LAT2METER + subfaults(k)%width .and. &
+                if (subfaults(k)%xcb(1) <= xcell .and. &
+                    xcell <= subfaults(k)%xcb(2) .and. &
                     subfaults(k)%rupture_time <= t .and. &
                     t <= subfaults(k)%rupture_time + subfaults(k)%rise_time) then
                   aux(13,i,j) = subfaults(k)%slip/subfaults(k)%rise_time
@@ -37,7 +37,7 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux)
 
             end if
           end do
- 
+
         end if
       end do
 
diff --git a/2d/fault_module.f90 b/2d/fault_module.f90
index 038167e..9fc68dc 100644
--- a/2d/fault_module.f90
+++ b/2d/fault_module.f90
@@ -5,6 +5,7 @@ module fault_module
     integer :: nsubfaults
     type subfault
       real(kind=8) :: width, depth, slip, longitude, rupture_time, rise_time
+      real(kind=8) :: xcb(2)
     end type subfault
     type(subfault), allocatable :: subfaults(:)
     real(kind=8) :: center(2), theta, xcb(2), final_rupture_rise_time
@@ -51,11 +52,22 @@ subroutine load_fault(fname)
 
         center(1) = 0.5d0*(xp1 + xp2)
         center(2) = 0.5d0*(yp1 + yp2)
-        total_width = dsqrt((xp2-xp1)**2 + (yp2-yp1)**2)      
+        total_width = dsqrt((xp2-xp1)**2 + (yp2-yp1)**2)
 
         xcb(1) = center(1) - 0.5d0*total_width
         xcb(2) = center(1) + 0.5d0*total_width
 
+        ! Determine subfault location in computational domain
+        print *, "*******WARNING*******"
+        print *, "Subfaults are assumed to be specified in top-center coords"
+        print *, "*********************"
+        do i=1,nsubfaults
+          subfaults(i)%xcb(1) = xcb(1) + dsqrt( &
+                (subfaults(i)%longitude*LAT2METER - xp1)**2 + &
+                (-subfaults(i)%depth - yp1)**2 )
+          subfaults(i)%xcb(2) = subfaults(i)%xcb(1) + subfaults(i)%width
+        end do
+
         close(7)
 
     end subroutine load_fault