Fixed small bug that would come up if using multiple subfaults with s…

…patially varying slip
cjvogl · Jan 9, 2018 · fd7095d · fd7095d
1 parent 61faacc
commit fd7095d
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 6 deletions.
diff --git a/2d/b4step2.f90 b/2d/b4step2.f90
@@ -3,15 +3,15 @@ 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
     real(kind=8), intent(in) :: xlower,ylower,dx,dy,t,dt
     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
@@ -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