[cig-commits] r19879 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src

[surendra at geodynamics.org]

Author: surendra
Date: 2012-03-26 16:57:04 -0700 (Mon, 26 Mar 2012)
New Revision: 19879

Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Added Fload distribution for TPV10x

Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-03-26 22:23:02 UTC (rev 19878)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-03-26 23:57:04 UTC (rev 19879)
@@ -548,6 +548,7 @@
   real(kind=CUSTOM_REAL) :: b(size(a))
   character(len=20) :: shape
   real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
+  real(kind=CUSTOM_REAL) :: r1(size(a))
   integer :: i
 
   NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
@@ -570,6 +571,13 @@
     select case(shape)
     case ('circle')
       b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
+    case ('circle-exp')
+      r1 = sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 )
+      where(r1<r)
+        b =exp (r1**2/(r1**2 - r**2) )
+      elsewhere
+        b =0._CUSTOM_REAL
+      endwhere
     case ('ellipse')
       b = heaviside( 1e0_CUSTOM_REAL - sqrt( (coord(1,:)-xc)**2/lx**2 + (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2 ) )
     case ('square')
@@ -713,24 +721,24 @@
     !WARNING : ad hoc for SCEC benchmark TPV10x
     TxExt = 0._CUSTOM_REAL
     TLoad = 1.0_CUSTOM_REAL
-    DTau0 = 25e6_CUSTOM_REAL
-    time = it*bc%dt
-    if (time==0._CUSTOM_REAL) then
-      GLoad = 0._CUSTOM_REAL
-    elseif (time <= TLoad) then
+    !DTau0 = 25e6_CUSTOM_REAL
+    time = it*bc%dt !time will never be zero. it starts from 1
+    if (time <= TLoad) then
       GLoad = exp( (time-TLoad)*(time-Tload) / (time*(time-2.0_CUSTOM_REAL*TLoad)) )
     else
       GLoad = 1.0_CUSTOM_REAL
     endif
-    TxExt = DTau0 * bc%asp%Fload * GLoad
-    T(1,:) = T(1,:) + TxExt
+    TxExt = (bc%asp%AspTx - bc%T0(1,:)) * bc%asp%Fload * GLoad + bc%T0(1,:)
+    !TxExt = DTau0 * bc%asp%Fload * GLoad
+    !T(1,:) = T(1,:) + TxExt
 
     Vf = Vnorm_old
     theta_old = bc%rsf%theta
     call rsf_update_state(Vf,bc%dt,bc%rsf)
 
     ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
-    tStick = sqrt(T(1,:)**2 + T(2,:)**2)
+    tStick = sqrt((T(1,:) - bc%T0(1,:) +TxExt )**2 + T(2,:)**2)
+    !tStick = sqrt(T(1,:)**2 + T(2,:)**2)
     do iNode=1,bc%nglob
       call NRsearch(Vf1(iNode),tau1(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),-T(3,iNode),log(bc%rsf%theta(iNode)),bc%Z(iNode),ta(iNode),tStick(iNode),ierr)