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

surendra at geodynamics.org surendra at geodynamics.org
Mon Mar 26 22:48:16 PDT 2012 

Author: surendra
Date: 2012-03-26 22:48:15 -0700 (Mon, 26 Mar 2012)
New Revision: 19880

Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Removed AspTx

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 23:57:04 UTC (rev 19879)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-03-27 05:48:15 UTC (rev 19880)
@@ -76,7 +76,7 @@
 
   type asperity_type
     private
-    real(kind=CUSTOM_REAL), dimension(:), pointer   :: AspTx=>null(), Fload=>null()
+    real(kind=CUSTOM_REAL), dimension(:), pointer   :: Fload=>null()
   end type asperity_type
 
   type bc_dynflt_type
@@ -230,8 +230,10 @@
   real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
   real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init
   integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init
-  real(kind=CUSTOM_REAL) :: AspTx,Fload
-  integer :: nAsp,nFload
+  real(kind=CUSTOM_REAL) :: Fload
+  integer :: nFload
+  real(kind=CUSTOM_REAL) :: C,T
+  integer :: nC,nForcedRup
 
 
   integer, parameter :: IIN_NUC =270
@@ -253,9 +255,9 @@
 
 
   NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
-  NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
-  NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init
-  NAMELIST / ASP / AspTx,Fload,nAsp,nFload
+  NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc,C,T,nC,nForcedRup
+  NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init,C,T,nC,nForcedRup
+  NAMELIST / ASP / Fload,nFload
 
   read(IIN_BIN) bc%nspec,bc%nglob
   if (bc%nspec==0) return
@@ -407,6 +409,8 @@
     allocate( bc%swf%mud(bc%nglob) )
     allocate( bc%swf%Dc(bc%nglob) )
     allocate( bc%swf%theta(bc%nglob) )
+    allocate( bc%swf%C(bc%nglob) )
+    allocate( bc%swf%T(bc%nglob) )
     ! WARNING: if V_HEALING is negative we turn off healing
     bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
 
@@ -417,13 +421,22 @@
     nmud = 0
     ndc  = 0
 
+    C = 0._CUSTOM_REAL
+    T = HUGEVAL
+    nC = 0
+    nForcedRup = 0
+
     read(IIN_PAR, nml=SWF)
     bc%swf%mus = mus
     bc%swf%mud = mud
     bc%swf%Dc  = dc
+    bc%swf%C   = C
+    bc%swf%T   = T
     call init_2d_distribution(bc%swf%mus,bc%coord,IIN_PAR,nmus)
     call init_2d_distribution(bc%swf%mud,bc%coord,IIN_PAR,nmud) 
     call init_2d_distribution(bc%swf%Dc ,bc%coord,IIN_PAR,ndc)
+    call init_2d_distribution(bc%swf%C  ,bc%coord,IIN_PAR,nC)
+    call init_2d_distribution(bc%swf%T  ,bc%coord,IIN_PAR,nForcedRup)
 
     bc%swf%theta = 0e0_CUSTOM_REAL
     allocate(bc%MU(bc%nglob))
@@ -439,11 +452,13 @@
     allocate( bc%rsf%L(bc%nglob) )
     allocate( bc%rsf%V_init(bc%nglob) )
     allocate( bc%rsf%theta(bc%nglob) )
+    allocate( bc%rsf%C(bc%nglob) )
+    allocate( bc%rsf%T(bc%nglob) )
 
     V0 =1.e-6_CUSTOM_REAL
     f0 =0.6_CUSTOM_REAL
-    a =0.0220_CUSTOM_REAL  !0.0080_CUSTOM_REAL
-    b =0.0050_CUSTOM_REAL  !0.0120_CUSTOM_REAL
+    a =0.0080_CUSTOM_REAL  !0.0080_CUSTOM_REAL
+    b =0.0040_CUSTOM_REAL  !0.0120_CUSTOM_REAL
     L =0.0135_CUSTOM_REAL
     V_init =1.e-12_CUSTOM_REAL
     theta_init =1.084207680000000e+09_CUSTOM_REAL
@@ -456,6 +471,11 @@
     nV_init =0
     ntheta_init =0
 
+    C = 0._CUSTOM_REAL
+    T = HUGEVAL
+    nC = 0
+    nForcedRup = 0
+
     read(IIN_PAR, nml=RSF)
     bc%rsf%V0 = V0
     bc%rsf%f0 = f0
@@ -464,6 +484,8 @@
     bc%rsf%L = L
     bc%rsf%V_init = V_init
     bc%rsf%theta = theta_init
+    bc%rsf%C   = C
+    bc%rsf%T   = T
     call init_2d_distribution(bc%rsf%V0,bc%coord,IIN_PAR,nV0)
     call init_2d_distribution(bc%rsf%f0,bc%coord,IIN_PAR,nf0)
     call init_2d_distribution(bc%rsf%a,bc%coord,IIN_PAR,na)
@@ -471,26 +493,22 @@
     call init_2d_distribution(bc%rsf%L,bc%coord,IIN_PAR,nL)
     call init_2d_distribution(bc%rsf%V_init,bc%coord,IIN_PAR,nV_init)
     call init_2d_distribution(bc%rsf%theta,bc%coord,IIN_PAR,ntheta_init)
+    call init_2d_distribution(bc%rsf%C  ,bc%coord,IIN_PAR,nC)
+    call init_2d_distribution(bc%rsf%T  ,bc%coord,IIN_PAR,nForcedRup)
 
-    allocate(bc%MU(bc%nglob)) !Not necessary?
-    bc%MU = rsf_mu(bc%rsf,bc%rsf%V_init)
-!    bc%T0(3,:) = bc%T0(3,:) + 1.e6_CUSTOM_REAL
+    allocate(bc%MU(bc%nglob)) 
+    bc%MU = rsf_mu(bc%rsf,bc%rsf%V_init) ! Should we switch to regularized form in rsf_mu?
     bc%T0(1,:) = bc%MU * bc%T0(3,:) 
     bc%T0(2,:) = 0._CUSTOM_REAL
 
     allocate( bc%asp )
-    allocate( bc%asp%AspTx(bc%nglob) )
     allocate( bc%asp%Fload(bc%nglob) )
 
-    AspTx =0.e0_CUSTOM_REAL
     Fload =0.e0_CUSTOM_REAL
-    nAsp =0
     nFload =0
 
     read(IIN_PAR, nml=ASP)
-    bc%asp%AspTx =AspTx
     bc%asp%Fload =Fload
-    call init_2d_distribution(bc%asp%AspTx,bc%coord,IIN_PAR,nAsp)
     call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
 
   endif
@@ -657,6 +675,7 @@
   real(kind=CUSTOM_REAL), dimension(bc%nglob) :: TxExt
   real(kind=CUSTOM_REAL) :: TLoad,DTau0,GLoad
   real(kind=CUSTOM_REAL) :: time
+  real(kind=CUSTOM_REAL) :: psi
 
   half_dt = 0.5e0_CUSTOM_REAL*bc%dt
   Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
@@ -721,26 +740,25 @@
     !WARNING : ad hoc for SCEC benchmark TPV10x
     TxExt = 0._CUSTOM_REAL
     TLoad = 1.0_CUSTOM_REAL
-    !DTau0 = 25e6_CUSTOM_REAL
+    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 = (bc%asp%AspTx - bc%T0(1,:)) * bc%asp%Fload * GLoad + bc%T0(1,:)
-    !TxExt = DTau0 * bc%asp%Fload * GLoad
-    !T(1,:) = T(1,:) + TxExt
+    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,:) - bc%T0(1,:) +TxExt )**2 + T(2,:)**2)
-    !tStick = sqrt(T(1,:)**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)
+      psi = log( bc%rsf%theta(iNode) * bc%rsf%V_init(iNode) / bc%rsf%L(iNode) )
+      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),psi,bc%Z(iNode),ta(iNode),tStick(iNode),ierr)
 
       ! If Vf1 or tau1 is negative, set Vf1 = 0 and tau1 = TauStick
       if (Vf1(iNode) < 0._CUSTOM_REAL .or. tau1(iNode) < 0._CUSTOM_REAL .or. ierr == 1) then
@@ -763,7 +781,8 @@
 
     ! NR search 2nd loop
     do iNode=1,bc%nglob
-      call NRsearch(Vf2(iNode),tau2(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)
+      psi = log( bc%rsf%theta(iNode) * bc%rsf%V_init(iNode) / bc%rsf%L(iNode) )
+      call NRsearch(Vf2(iNode),tau2(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),-T(3,iNode),psi,bc%Z(iNode),ta(iNode),tStick(iNode),ierr)
 
       if (Vf2(iNode) < 0._CUSTOM_REAL .or. tau2(iNode) < 0._CUSTOM_REAL .or. ierr == 1) then
         tau2(iNode) = tStick(iNode)

More information about the CIG-COMMITS mailing list