[cig-commits] r19940 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
ampuero at geodynamics.org
ampuero at geodynamics.org
Thu Apr 12 14:26:05 PDT 2012
Author: ampuero
Date: 2012-04-12 14:26:05 -0700 (Thu, 12 Apr 2012)
New Revision: 19940
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
Log:
removed reallocation of T0 in TPV16
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-04-12 18:49:24 UTC (rev 19939)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-04-12 21:26:05 UTC (rev 19940)
@@ -1,1410 +1,1409 @@
-!=====================================================================
-!
-! s p e c f e m 3 d v e r s i o n 1 . 4
-! ---------------------------------------
-!
-! dimitri komatitsch and jeroen tromp
-! seismological laboratory - california institute of technology
-! (c) california institute of technology september 2006
-!
-! this program is free software; you can redistribute it and/or modify
-! it under the terms of the gnu general public license as published by
-! the free software foundation; either version 2 of the license, or
-! (at your option) any later version.
-!
-! this program is distributed in the hope that it will be useful,
-! but without any warranty; without even the implied warranty of
-! merchantability or fitness for a particular purpose. see the
-! gnu general public license for more details.
-!
-! you should have received a copy of the gnu general public license along
-! with this program; if not, write to the free software foundation, inc.,
-! 51 franklin street, fifth floor, boston, ma 02110-1301 usa.
-!
-!===============================================================================================================
-
-! This module was written by:
-! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
-! Surendra Nadh Somala : Added Heterogenous initial stress capabilities (based on TPV16)
-! Surendra Nadh Somala : Added Rate and State Friction
-
-module fault_solver
-
- implicit none
-
- include 'constants.h'
-
- private
-
- ! outputs on selected fault nodes at every time step:
- ! slip, slip velocity, fault stresses
- type dataT_type
- integer :: npoin
- integer, dimension(:), pointer :: iglob ! on-fault global index of output nodes
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1,v1,t1,d2,v2,t2,t3,theta
- character(len=70), dimension(:), pointer :: name
- end type dataT_type
-
-
- ! outputs at selected times for all fault nodes:
- ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
- ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
- ! process zone time = first time when slip = Dc
- type dataXZ_type
- real(kind=CUSTOM_REAL), dimension(:), pointer :: stg, sta, d1, d2, v1, v2, &
- t1, t2, t3, tRUP,tPZ
- real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord,ycoord,zcoord
- integer :: npoin
- end type dataXZ_type
-
- type swf_type
- private
- integer :: kind
- logical :: healing = .false.
- real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), &
- theta=>null(), T=>null(), C=>null()
- end type swf_type
-
- type rsf_type
- private
- integer :: StateLaw = 1 ! 1=ageing law, 2=slip law
- real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), &
- V_init=>null(), &
- a=>null(), b=>null(), theta=>null(), &
- T=>null(), C=>null()
- end type rsf_type
-
- type asperity_type
- private
- real(kind=CUSTOM_REAL), dimension(:), pointer :: Fload=>null()
- end type asperity_type
-
- type bc_dynflt_type
- private
- integer :: nspec,nglob
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0,T,V,D
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: coord
- real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R
- real(kind=CUSTOM_REAL), dimension(:), pointer :: MU,B,invM1,invM2,Z
- real(kind=CUSTOM_REAL) :: dt
- integer, dimension(:), pointer :: ibulk1, ibulk2
- type(swf_type), pointer :: swf => null()
- type(rsf_type), pointer :: rsf => null()
- type(asperity_type), pointer :: asp => null()
- logical :: allow_opening = .false. ! default : do not allow opening
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ
- end type bc_dynflt_type
-
- type(bc_dynflt_type), allocatable, save :: faults(:)
-
- !slip velocity threshold for healing
- !WARNING: not very robust
- real(kind=CUSTOM_REAL), save :: V_HEALING
-
- !slip velocity threshold for definition of rupture front
- real(kind=CUSTOM_REAL), save :: V_RUPT
-
- !Number of time steps defined by the user : NTOUT
- integer, save :: NTOUT,NSNAP
-
- logical, save :: SIMULATION_TYPE_DYN = .false.
-
- logical, save :: TPV16 = .false.
-
- logical, save :: Rate_AND_State = .true.
-
- real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
-
- public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
- SIMULATION_TYPE_DYN
-
-
-contains
-
-!=====================================================================
-! BC_DYNFLT_init initializes dynamic faults
-!
-! prname fault database is read from file prname_fault_db.bin
-! Minv inverse mass matrix
-! dt global time step
-!
-subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt)
-
- character(len=256), intent(in) :: prname ! 'proc***'
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
- double precision, intent(in) :: DTglobal
- integer, intent(in) :: nt
-
- real(kind=CUSTOM_REAL) :: dt
- integer :: iflt,ier,dummy_idfault
- integer :: nbfaults
- integer :: size_Kelvin_Voigt
- integer :: SIMULATION_TYPE
- character(len=256) :: filename
- integer, parameter :: IIN_PAR =151
- integer, parameter :: IIN_BIN =170
-
- NAMELIST / BEGIN_FAULT / dummy_idfault
-
- dummy_idfault = 0
-
- open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults',status='old',iostat=ier)
- if( ier /= 0 ) then
- write(6,*) 'File Par_file_faults not found: assume no faults'
- close(IIN_PAR)
- return
- endif
-
- dt = real(DTglobal)
- filename = prname(1:len_trim(prname))//'fault_db.bin'
- open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- write(6,*) 'File ',trim(filename),' not found. Abort'
- stop
- endif
- ! WARNING TO DO: should be an MPI abort
-
- read(IIN_PAR,*) nbfaults
- if (nbfaults==0) return
- ! Reading etas of each fault
- do iflt = 1,nbfaults
- read(IIN_PAR,*) ! etas
- enddo
- read(IIN_PAR,*) SIMULATION_TYPE
- if ( SIMULATION_TYPE /= 1 ) then
- close(IIN_BIN)
- close(IIN_PAR)
- return
- endif
- SIMULATION_TYPE_DYN = .true.
- read(IIN_PAR,*) NTOUT
- read(IIN_PAR,*) NSNAP
- read(IIN_PAR,*) V_HEALING
- read(IIN_PAR,*) V_RUPT
-
- read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
- allocate( faults(nbfaults) )
- do iflt=1,nbfaults
- read(IIN_PAR,nml=BEGIN_FAULT,end=100)
- call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
- enddo
- close(IIN_BIN)
- close(IIN_PAR)
-
- filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
- open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- write(6,*) 'File ',trim(filename),' not found. Abort'
- stop
- endif
- read(IIN_BIN) size_Kelvin_Voigt
- if (size_Kelvin_Voigt > 0) then
- allocate(Kelvin_Voigt_eta(size_Kelvin_Voigt))
- read(IIN_BIN) Kelvin_Voigt_eta
- endif
- close(IIN_BIN)
- return
-100 stop 'Did not find BEGIN_FAULT block #'
- ! WARNING TO DO: should be an MPI abort
-
-end subroutine BC_DYNFLT_init
-
-!---------------------------------------------------------------------
-
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
-
- type(bc_dynflt_type), intent(inout) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
- integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: dt
-
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
- integer, dimension(:,:), allocatable :: ibool1
- real(kind=CUSTOM_REAL) :: norm
- real(kind=CUSTOM_REAL) :: S1,S2,S3
- integer :: n1,n2,n3
- real(kind=CUSTOM_REAL) :: mus,mud,dc
- integer :: nmus,nmud,ndc,ij,k,e
- 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) :: C,T
- integer :: nC,nForcedRup
-
-
- integer, parameter :: IIN_NUC =270
- integer :: ipar
- integer :: ier
-
- integer :: relz_num,sub_relz_num
- integer :: num_cell_str,num_cell_dip
- real(kind=CUSTOM_REAL) :: siz_str,siz_dip
- integer :: hypo_cell_str,hypo_cell_dip
- real(kind=CUSTOM_REAL) :: hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
-
- integer, dimension(:), allocatable :: inp_nx,inp_nz
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
- dyn_fc,swcd,cohes,tim_forcedRup
-
- integer :: iLoc
- real(kind=CUSTOM_REAL) :: minX
-
-
-
- real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
- real(kind=CUSTOM_REAL) :: x,z
- logical :: c1,c2,c3,c4
- real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
- integer iglob
-
- double precision, dimension(:), allocatable :: mu1,mu2,mu3
-
-
- NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- 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
- allocate( bc%ibulk1(bc%nglob) )
- allocate( bc%ibulk2(bc%nglob) )
- allocate( ibool1(NGLLSQUARE,bc%nspec) )
- allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
- allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
-
- allocate(bc%coord(3,(bc%nglob)))
- read(IIN_BIN) ibool1
- read(IIN_BIN) jacobian2Dw
- read(IIN_BIN) normal
- read(IIN_BIN) bc%ibulk1
- read(IIN_BIN) bc%ibulk2
- read(IIN_BIN) bc%coord(1,:)
- read(IIN_BIN) bc%coord(2,:)
- read(IIN_BIN) bc%coord(3,:)
- bc%dt = dt
-
- allocate( bc%B(bc%nglob) )
- bc%B = 0e0_CUSTOM_REAL
- allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
- nx = 0e0_CUSTOM_REAL
- ny = 0e0_CUSTOM_REAL
- nz = 0e0_CUSTOM_REAL
- do e=1,bc%nspec
- do ij = 1,NGLLSQUARE
- k = ibool1(ij,e)
- nx(k) = nx(k) + normal(1,ij,e)
- ny(k) = ny(k) + normal(2,ij,e)
- nz(k) = nz(k) + normal(3,ij,e)
- bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
- enddo
- enddo
- do k=1,bc%nglob
- norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
- nx(k) = nx(k) / norm
- ny(k) = ny(k) / norm
- nz(k) = nz(k) / norm
- enddo
-
- allocate( bc%R(3,3,bc%nglob) )
- call compute_R(bc%R,bc%nglob,nx,ny,nz)
-
- ! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
- allocate(bc%invM1(bc%nglob))
- allocate(bc%invM2(bc%nglob))
- bc%invM1 = Minv(bc%ibulk1)
- bc%invM2 = Minv(bc%ibulk2)
-
- ! Fault impedance, Z in : Trac=T_Stick-Z*dV
- ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
- ! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
- ! NOTE: same Bi on both sides, see note above
- allocate(bc%Z(bc%nglob))
- bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
-
- allocate(bc%T(3,bc%nglob))
- allocate(bc%D(3,bc%nglob))
- allocate(bc%V(3,bc%nglob))
- bc%T = 0e0_CUSTOM_REAL
- bc%D = 0e0_CUSTOM_REAL
- bc%V = 0e0_CUSTOM_REAL
-
- ! Set initial fault stresses
- allocate(bc%T0(3,bc%nglob))
- S1 = 0e0_CUSTOM_REAL
- S2 = 0e0_CUSTOM_REAL
- S3 = 0e0_CUSTOM_REAL
- n1=0
- n2=0
- n3=0
- read(IIN_PAR, nml=INIT_STRESS)
- bc%T0(1,:) = S1
- bc%T0(2,:) = S2
- bc%T0(3,:) = S3
-
- call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
- call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
- call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
-
- !WARNING : Quick and dirty free surface condition at z=0
- ! do k=1,bc%nglob
- ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
- ! end do
-
-
- if (TPV16) then
-
- allocate(inp_nx(bc%nglob))
- allocate(inp_nz(bc%nglob))
- allocate(loc_str(bc%nglob))
- allocate(loc_dip(bc%nglob))
- allocate(sigma0(bc%nglob))
- allocate(tau0_str(bc%nglob))
- allocate(tau0_dip(bc%nglob))
- allocate(Rstress_str(bc%nglob))
- allocate(Rstress_dip(bc%nglob))
- allocate(static_fc(bc%nglob))
- allocate(dyn_fc(bc%nglob))
- allocate(swcd(bc%nglob))
- allocate(cohes(bc%nglob))
- allocate(tim_forcedRup(bc%nglob))
-
- open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
- read(IIN_NUC,*) relz_num,sub_relz_num
- read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
- read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
- do ipar=1,bc%nglob
- read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
- Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
- enddo
- close(IIN_NUC)
-
- allocate(bc%T0(3,bc%nglob))
- allocate( bc%swf )
- allocate( bc%swf%mus(bc%nglob) )
- allocate( bc%swf%mud(bc%nglob) )
- allocate( bc%swf%Dc(bc%nglob) )
- allocate( bc%swf%theta(bc%nglob) )
- allocate( bc%swf%T(bc%nglob) )
- allocate( bc%swf%C(bc%nglob) )
-
- minX = minval(bc%coord(1,:))
-
- do iLoc=1,bc%nglob
-
- ipar = minloc( (minX+loc_str(:)-bc%coord(1,iLoc))**2 + (-loc_dip(:)-bc%coord(3,iLoc))**2 , 1) !iloc_dip is negative of Z-coord
- bc%T0(3,iLoc) = -sigma0(ipar)
- bc%T0(1,iLoc) = tau0_str(ipar)
- bc%T0(2,iLoc) = tau0_dip(ipar)
-
- bc%swf%mus(iLoc) = static_fc(ipar)
- bc%swf%mud(iLoc) = dyn_fc(ipar)
- bc%swf%Dc(iLoc) = swcd(ipar)
- bc%swf%C(iLoc) = cohes(ipar)
- bc%swf%T(iLoc) = tim_forcedRup(ipar)
- enddo
-
- endif
-
-
- ! Set friction parameters and initialize friction variables
- ! Slip weakening friction
- if(.not. Rate_AND_State) then
- allocate( bc%swf )
- allocate( bc%swf%mus(bc%nglob) )
- 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)
-
- mus = 0.6e0_CUSTOM_REAL
- mud = 0.1e0_CUSTOM_REAL
- dc = 1e0_CUSTOM_REAL
- nmus = 0
- 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))
- bc%MU = swf_mu(bc%swf)
-
- ! Rate and state friction
- else
- allocate( bc%rsf )
- allocate( bc%rsf%V0(bc%nglob) )
- allocate( bc%rsf%f0(bc%nglob) )
- allocate( bc%rsf%a(bc%nglob) )
- allocate( bc%rsf%b(bc%nglob) )
- 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.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
-
- nV0 =0
- nf0 =0
- na =0
- nb =0
- nL =0
- 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
- bc%rsf%a = a
- bc%rsf%b = b
- 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)
- call init_2d_distribution(bc%rsf%b,bc%coord,IIN_PAR,nb)
- 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)
-
-
- W1=15000._CUSTOM_REAL
- W2=7500._CUSTOM_REAL
- w=3000._CUSTOM_REAL
- hypo_z = -7500._CUSTOM_REAL
- do iglob=1,bc%nglob
- x=bc%coord(1,iglob)
- z=bc%coord(3,iglob)
- c1=abs(x)<W1+w
- c2=abs(x)>W1
- c3=abs(z-hypo_z)<W2+w
- c4=abs(z-hypo_z)>W2
- if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
-
- if (c1 .and. c2) then
- b11 = w/(abs(x)-W1-w)
- b12 = w/(abs(x)-W1)
- B1 = 0.5 * (ONE + tanh(b11 + b12))
- elseif(abs(x)<=W1) then
- B1 = 1._CUSTOM_REAL
- else
- B1 = 0._CUSTOM_REAL
- endif
-
- if (c3 .and. c4) then
- b21 = w/(abs(z-hypo_z)-W2-w)
- b22 = w/(abs(z-hypo_z)-W2)
- B2 = 0.5 * (ONE + tanh(b21 + b22))
- elseif(abs(z-hypo_z)<=W2) then
- B2 = 1._CUSTOM_REAL
- else
- B2 = 0._CUSTOM_REAL
- endif
-
-
- bc%rsf%a(iglob) = 0.008 + 0.008 * (ONE - B1*B2)
- elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
- bc%rsf%a(iglob) = 0.008
- else
- bc%rsf%a(iglob) = 0.016
- endif
-
- bc%rsf%theta(iglob) = L/V0 * exp( (bc%rsf%a(iglob) * log(2.0*sinh(-S1/S3/bc%rsf%a(iglob))) - f0 - bc%rsf%a(iglob)*log(V_init/V0) )/b )
-
-
- enddo
-
-
-
- allocate(bc%MU(bc%nglob))
-
- bc%V(1,:) = bc%rsf%V_init
- allocate( bc%asp )
- allocate( bc%asp%Fload(bc%nglob) )
-
- Fload =0.e0_CUSTOM_REAL
- nFload =0
-
- read(IIN_PAR, nml=ASP)
- bc%asp%Fload =Fload
- call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
-
- endif
-
-
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc,bc%nglob)
-
-end subroutine init_one_fault
-
-!---------------------------------------------------------------------
-subroutine compute_R(R,nglob,nx,ny,nz)
-
- integer :: nglob
- real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
- real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
-
- real(kind=CUSTOM_REAL), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
-
- ! Percy , defining fault directions (in concordance with SCEC conventions) .
- ! fault coordinates (s,d,n) = (1,2,3)
- ! s = strike , d = dip , n = n.
- ! 1 = strike , 2 = dip , 3 = n.
- norm = sqrt(nx*nx+ny*ny)
- sx = ny/norm
- sy = -nx/norm
- sz = 0.e0_CUSTOM_REAL
-
- norm = sqrt(sy*sy*nz*nz+sx*sx*nz*nz+(sy*nx-ny*sx)*(nx*sy-ny*sx))
- dx = -sy*nz/norm
- dy = sx*nz/norm
- dz = (sy*nx-ny*sx)/norm
- !Percy, dz is always dipwards = -1/norm , because (nx*sy-ny*sx)= - 1
-
- R(1,1,:)=sx
- R(1,2,:)=sy
- R(1,3,:)=sz
- R(2,1,:)=dx
- R(2,2,:)=dy
- R(2,3,:)=dz
- R(3,1,:)=nx
- R(3,2,:)=ny
- R(3,3,:)=nz
-
-end subroutine compute_R
-
-!---------------------------------------------------------------------
-! adds a value to a fault parameter inside an area with prescribed shape
-subroutine init_2d_distribution(a,coord,iin,n)
-
- real(kind=CUSTOM_REAL), intent(inout) :: a(:)
- real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
- integer, intent(in) :: iin,n
-
- 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
-
- if (n==0) return
-
- do i=1,n
- shape = ''
- xc = 0e0_CUSTOM_REAL
- yc = 0e0_CUSTOM_REAL
- zc = 0e0_CUSTOM_REAL
- r = 0e0_CUSTOM_REAL
- l = 0e0_CUSTOM_REAL
- lx = 0e0_CUSTOM_REAL
- ly = 0e0_CUSTOM_REAL
- lz = 0e0_CUSTOM_REAL
- valh = 0e0_CUSTOM_REAL
-
- read(iin,DIST2D)
- 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')
- b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('rectangle')
- b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('rectangle-taper')
- b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
- case default
- stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
- end select
-
- where (b /= 0e0_CUSTOM_REAL) a = b
- enddo
-
-end subroutine init_2d_distribution
-
-!---------------------------------------------------------------------
-elemental function heaviside(x)
-
- real(kind=CUSTOM_REAL), intent(in) :: x
- real(kind=CUSTOM_REAL) :: heaviside
-
- if (x>=0e0_CUSTOM_REAL) then
- heaviside = 1e0_CUSTOM_REAL
- else
- heaviside = 0e0_CUSTOM_REAL
- endif
-
-end function heaviside
-
-!=====================================================================
-! adds boundary term Bt into Force array for each fault.
-!
-subroutine bc_dynflt_set3d_all(F,Vel,Dis)
-
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
- real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
-
- integer :: iflt
-
- if (.not. allocated(faults)) return
- do iflt=1,size(faults)
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
- enddo
-
-end subroutine bc_dynflt_set3d_all
-
-!---------------------------------------------------------------------
-subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
-
- use specfem_par, only:it,NSTEP
-
- real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
- type(bc_dynflt_type), intent(inout) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
- integer,intent(in) :: iflt
-
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,tnew
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, theta_new, Vnorm, Vnorm_old, dc
- real(kind=CUSTOM_REAL) :: half_dt
-
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: ta, Vf,Vf1,tau1,Vf2,tau2,Vfavg
- integer :: ierr,iNode
- 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,:))
-
- ! get predicted values
- dD = get_jump(bc,D) ! dD_predictor
- dV = get_jump(bc,V) ! dV_predictor
- dA = get_weighted_jump(bc,MxA) ! dA_free
-
- ! rotate to fault frame (tangent,normal)
- ! component 3 is normal to the fault
- dD = rotate(bc,dD,1)
- dV = rotate(bc,dV,1)
- dA = rotate(bc,dA,1)
-
- ! T_stick
- T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
- T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
- T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
-
- !Warning : dirty particular free surface condition z = 0.
- ! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
- ! do k=1,bc%nglob
- ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
- ! end do
-
- ! add initial stress
- T = T + bc%T0
-
- ! Solve for normal stress (negative is compressive)
- ! Opening implies free stress
- if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
-
- if(.not. Rate_AND_State) then ! Update slip weakening friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
- theta_old = bc%swf%theta
- call swf_update_state(bc%D,dD,bc%V,bc%swf)
-
- ! Update friction coeficient
- bc%MU = swf_mu(bc%swf)
-
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
- if (TPV16) then
- where (bc%swf%T <= it*bc%dt) bc%MU = bc%swf%mud
- endif
-
- ! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
-
- ! Solve for shear stress
- tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
- tnew = min(tStick,strength)
- tStick = max(tStick,1e0_CUSTOM_REAL)
- T(1,:) = tnew * T(1,:)/tStick
- T(2,:) = tnew * T(2,:)/tStick
-
- else ! Update rate and state friction:
-
- ! smooth loading within nucleation patch
- !WARNING : ad hoc for SCEC benchmark TPV10x
- TxExt = 0._CUSTOM_REAL
- TLoad = 1.0_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 = 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)
-
- tStick = sqrt(T(1,:)**2 + T(2,:)**2)
- do iNode=1,bc%nglob
- Vf1(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
- tau1(iNode) = tStick(iNode) - bc%Z(iNode)*Vf1(iNode)
- enddo
-
- ! Updating state variable: 2nd loop
- Vfavg = 0.5e0_CUSTOM_REAL*(Vf + Vf1)
- bc%rsf%theta = theta_old
- call rsf_update_state(Vfavg,bc%dt,bc%rsf)
-
- ! NR search 2nd loop
- do iNode=1,bc%nglob
- Vf2(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
- tau2(iNode) = tStick(iNode) - bc%Z(iNode)*Vf2(iNode)
- enddo
-
- tStick = max(tStick,1e0_CUSTOM_REAL)
- T(1,:) = tau2 * T(1,:)/tStick
- T(2,:) = tau2 * T(2,:)/tStick
-
- endif
-
- ! Save total tractions
- bc%T = T
-
- ! Subtract initial stress
- T = T - bc%T0
-
- ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
- dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
- dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
- dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
-
- ! Update slip and slip rate, in fault frame
- bc%D = dD
- bc%V = dV + half_dt*dA
-
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
-
- ! Add boundary term B*T to M*a
- MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
- MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
- MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
-
- MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
- MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
- MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
-
- !-- intermediate storage of outputs --
- Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- if(.not. Rate_AND_State) then
- theta_new = bc%swf%theta
- dc = bc%swf%dc
- else
- theta_new = bc%rsf%theta
- dc = bc%rsf%L
- endif
-
- call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
- Vnorm_old, Vnorm, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
-
- !-- outputs --
- ! write dataT every NTOUT time step or at the end of simulation
- if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
- ! write dataXZ every NSNAP time step
- if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
- if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
-
-end subroutine BC_DYNFLT_set3d
-
-!===============================================================
-function get_jump (bc,v) result(dv)
-
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
- real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
-
- ! diference between side 2 and side 1 of fault nodes. dv
- dv(1,:) = v(1,bc%ibulk2)-v(1,bc%ibulk1)
- dv(2,:) = v(2,bc%ibulk2)-v(2,bc%ibulk1)
- dv(3,:) = v(3,bc%ibulk2)-v(3,bc%ibulk1)
-
-end function get_jump
-
-!---------------------------------------------------------------------
-function get_weighted_jump (bc,f) result(da)
-
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
-
- real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
-
- ! diference between side 2 and side 1 of fault nodes. M-1 * F
- da(1,:) = bc%invM2*f(1,bc%ibulk2)-bc%invM1*f(1,bc%ibulk1)
- da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
- da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
-
-end function get_weighted_jump
-
-!----------------------------------------------------------------------
-function rotate(bc,v,fb) result(vr)
-
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: v(3,bc%nglob)
- integer, intent(in) :: fb
- real(kind=CUSTOM_REAL) :: vr(3,bc%nglob)
-
- ! Percy, tangential direction Vt, equation 7 of Pablo's notes in agreement with SPECFEM3D
-
- ! forward rotation
- if (fb==1) then
- vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(1,2,:)+v(3,:)*bc%R(1,3,:) ! vs
- vr(2,:) = v(1,:)*bc%R(2,1,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(2,3,:) ! vd
- vr(3,:) = v(1,:)*bc%R(3,1,:)+v(2,:)*bc%R(3,2,:)+v(3,:)*bc%R(3,3,:) ! vn
-
- ! backward rotation
- else
- vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(2,1,:)+v(3,:)*bc%R(3,1,:) !vx
- vr(2,:) = v(1,:)*bc%R(1,2,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(3,2,:) !vy
- vr(3,:) = v(1,:)*bc%R(1,3,:)+v(2,:)*bc%R(2,3,:)+v(3,:)*bc%R(3,3,:) !vz
-
- endif
-
-end function rotate
-
-
-!=====================================================================
-subroutine swf_update_state(dold,dnew,vold,f)
-
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
- type(swf_type), intent(inout) :: f
-
- real(kind=CUSTOM_REAL) :: vnorm
- integer :: k,npoin
-
- f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
-
- if (f%healing) then
- npoin = size(vold,2)
- do k=1,npoin
- vnorm = sqrt(vold(1,k)**2 + vold(2,k)**2)
- if (vnorm<V_HEALING) f%theta(k) = 0e0_CUSTOM_REAL
- enddo
- endif
-end subroutine swf_update_state
-
-!---------------------------------------------------------------------
-function swf_mu(f) result(mu)
-
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
-
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- mu = max( mu, f%mud)
-
-end function swf_mu
-
-!=====================================================================
-! Rate and state friction coefficient
-function rsf_mu(f,V) result(mu)
-
- type(rsf_type), intent(in) :: f
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
- real(kind=CUSTOM_REAL) :: mu(size(V))
- mu = f%a * asinh( V/2.0/f%V0 * exp((f%f0 + f%b*log(f%theta*f%V0/f%L))/f%a ) ) ! Regularized
-
-end function rsf_mu
-
-!---------------------------------------------------------------------
-subroutine funcd(x,fn,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
-
- real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
- double precision :: arg,fn,df,x
-
- arg = exp((f0+dble(b)*log(V0*theta/L))/a)/2._CUSTOM_REAL/V0
- fn = tStick - Z*x - a*Seff*asinh(x*arg)
- df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
-end subroutine funcd
-
-!---------------------------------------------------------------------
-function rtsafe(funcd,x1,x2,xacc,tStick,Seff,Z,f0,V0,a,b,L,theta)
-
- integer, parameter :: MAXIT=200
- real(kind=CUSTOM_REAL) :: x1,x2,xacc
- EXTERNAL funcd
- integer :: j
- !real(kind=CUSTOM_REAL) :: df,dx,dxold,f,fh,fl,temp,xh,xl
- double precision :: df,dx,dxold,f,fh,fl,temp,xh,xl,rtsafe
- real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
-
- call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- if( (fl>0 .and. fh>0) .or. (fl<0 .and. fh<0) ) stop 'root must be bracketed in rtsafe'
- if(fl==0.) then
- rtsafe=x2
- return
- elseif(fh==0.) then
- rtsafe=x2
- return
- elseif(fl<0) then
- xl=x1
- xh=x2
- else
- xh=x1
- xl=x2
- endif
-
- rtsafe=0.5d0*(x1+x2)
- dxold=abs(x2-x1)
- dx=dxold
- call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- do j=1,MAXIT
- if( ((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f)>0 .or. abs(2.*f)>abs(dxold*df) ) then
- dxold=dx
- dx=0.5d0*(xh-xl)
- rtsafe=xl+dx
- if(xl==rtsafe) return
- else
- dxold=dx
- dx=f/df
- temp=rtsafe
- rtsafe=rtsafe-dx
- if(temp==rtsafe) return
- endif
- if(abs(dx)<xacc) return
- call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- if(f<0.) then
- xl=rtsafe
- else
- xh=rtsafe
- endif
- enddo
- stop 'rtsafe exceeding maximum iterations'
- return
-
-end function rtsafe
-
-!---------------------------------------------------------------------
-subroutine rsf_update_state(V,dt,f)
-
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
- type(rsf_type), intent(inout) :: f
- real(kind=CUSTOM_REAL), intent(in) :: dt
-
- real(kind=CUSTOM_REAL) :: vDtL(size(V))
-
- vDtL = V * dt / f%L
-
- ! ageing law
- if (f%StateLaw == 1) then
- where(vDtL > 1.e-5_CUSTOM_REAL)
- f%theta = f%theta * exp(-vDtL) + f%L/V * (ONE - exp(-vDtL))
- elsewhere
- f%theta = f%theta * exp(-vDtL) + dt*( ONE - HALF*vDtL )
- endwhere
-
- ! slip law
- else
- f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
- endif
-
-end subroutine rsf_update_state
-
-!===============================================================
-! OUTPUTS
-subroutine init_dataT(DataT,coord,nglob,NT,iflt)
- ! NT = total number of time steps
-
- integer, intent(in) :: nglob,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
- type (dataT_type), intent(out) :: DataT
-
- real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
- integer :: i, iglob , IIN, ier, jflt, np, k
- character(len=70) :: tmpname
-
- ! 1. read fault output coordinates from user file,
- ! 2. define iglob: the fault global index of the node nearest to user
- ! requested coordinate
-
- IIN = 251
- open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
- read(IIN,*) np
- DataT%npoin =0
- do i=1,np
- read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
- if (jflt==iflt) DataT%npoin = DataT%npoin +1
- enddo
- close(IIN)
-
- if (DataT%npoin == 0) return
-
- allocate(DataT%iglob(DataT%npoin))
- allocate(DataT%name(DataT%npoin))
-
- open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
- if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
- read(IIN,*) np
- k = 0
- do i=1,np
- read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
- if (jflt/=iflt) cycle
- k = k+1
- DataT%name(k) = tmpname
- !search nearest node
- distkeep = huge(distkeep)
-
- do iglob=1,nglob
- dist = sqrt((coord(1,iglob)-xtarget)**2 &
- + (coord(2,iglob)-ytarget)**2 &
- + (coord(3,iglob)-ztarget)**2)
- if (dist < distkeep) then
- distkeep = dist
- DataT%iglob(k) = iglob
- endif
- enddo
- enddo
-
- ! 3. allocate arrays and set to zero
- allocate(DataT%d1(NT,DataT%npoin))
- allocate(DataT%v1(NT,DataT%npoin))
- allocate(DataT%t1(NT,DataT%npoin))
- allocate(DataT%d2(NT,DataT%npoin))
- allocate(DataT%v2(NT,DataT%npoin))
- allocate(DataT%t2(NT,DataT%npoin))
- allocate(DataT%t3(NT,DataT%npoin))
- allocate(DataT%theta(NT,DataT%npoin))
- DataT%d1 = 0e0_CUSTOM_REAL
- DataT%v1 = 0e0_CUSTOM_REAL
- DataT%t1 = 0e0_CUSTOM_REAL
- DataT%d2 = 0e0_CUSTOM_REAL
- DataT%v2 = 0e0_CUSTOM_REAL
- DataT%t2 = 0e0_CUSTOM_REAL
- DataT%t3 = 0e0_CUSTOM_REAL
- DataT%theta = 0e0_CUSTOM_REAL
-
- close(IIN)
-
-end subroutine init_dataT
-
-!---------------------------------------------------------------
-subroutine store_dataT(dataT,d,v,t,theta,itime)
-
- type(dataT_type), intent(inout) :: dataT
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
- integer, intent(in) :: itime
-
- integer :: i,k
-
- do i=1,dataT%npoin
- k = dataT%iglob(i)
- dataT%d1(itime,i) = d(1,k)
- dataT%d2(itime,i) = d(2,k)
- dataT%v1(itime,i) = v(1,k)
- dataT%v2(itime,i) = v(2,k)
- dataT%t1(itime,i) = t(1,k)
- dataT%t2(itime,i) = t(2,k)
- dataT%t3(itime,i) = t(3,k)
- dataT%theta(itime,i) = theta(k)
- enddo
-
-end subroutine store_dataT
-
-!-----------------------------------------------------------------
-subroutine write_dataT_all(nt)
-
- integer, intent(in) :: nt
-
- integer :: i
-
- if (.not.allocated(faults)) return
- do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
- enddo
-
-end subroutine write_dataT_all
-
-!------------------------------------------------------------------------
-subroutine SCEC_write_dataT(dataT,DT,NT)
-
- type(dataT_type), intent(in) :: dataT
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT
-
- integer :: i,k,IOUT
- character :: NTchar*5
- integer :: today(3), now(3)
-
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
-
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
-
- write(NTchar,1) NT
- NTchar = adjustl(NTchar)
-
-1 format(I5)
- do i=1,dataT%npoin
-
- open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "# problem=TPV102"
- write(IOUT,*) "# author=Surendra Nadh Somala"
- write(IOUT,1000) today(2), today(1), today(3), now
- write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
- write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
- write(IOUT,*) "# time_step=",DT
- write(IOUT,*) "# location=",trim(dataT%name(i))
- write(IOUT,*) "# Column #1 = Time (s)"
- write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
- write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
- write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
- write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
- write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
- write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
- write(IOUT,*) "# Column #8 = normal stress (MPa)"
- if(Rate_AND_State) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
- write(IOUT,*) "#"
- write(IOUT,*) "# The line below lists the names of the data fields:"
- if(.not. Rate_AND_State) then
- write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
- write(IOUT,*) "#"
- do k=1,NT
- write(IOUT,'(8(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
- -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
- dataT%t3(k,i)/1.0e6_CUSTOM_REAL
- enddo
- else
- write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
- write(IOUT,*) "#"
- do k=1,NT
- write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
- -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
- dataT%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
- enddo
- endif
- close(IOUT)
-
-
-1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
-
-
- enddo
-
-end subroutine SCEC_write_dataT
-
-!-------------------------------------------------------------------------------------------------
-subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
-
- type(dataXZ_type), intent(in) :: dataXZ
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT,iflt
-
- integer :: i,IOUT
- character(len=70) :: filename
- integer*4 today(3), now(3)
-
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
-
-
- write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
-
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
-
- open(IOUT,file=trim(filename),status='replace')
- write(IOUT,*) "# problem=TPV102"
- write(IOUT,*) "# author=Surendra Nadh Somala"
- write(IOUT,1000) today(2), today(1), today(3), now
- write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
- write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
- write(IOUT,*) "# Column #1 = horizontal coordinate, distance along strike (m)"
- write(IOUT,*) "# Column #2 = vertical coordinate, distance down-dip (m)"
- write(IOUT,*) "# Column #3 = rupture time (s)"
- write(IOUT,*) "# "
- write(IOUT,*) "j k t"
- do i = 1,size(dataXZ%tRUP)
- write(IOUT,'(3(E15.7))') dataXZ%xcoord(i), -dataXZ%zcoord(i), dataXZ%tRUP(i)
- end do
-
- close(IOUT)
-
-1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
-
-end subroutine SCEC_Write_RuptureTime
-
-!-------------------------------------------------------------------------------------------------
-subroutine init_dataXZ(DataXZ,bc,nglob)
-
- type(dataXZ_type), intent(inout) :: DataXZ
- type(bc_dynflt_type) :: bc
- integer, intent(in) :: nglob
-
- allocate(DataXZ%stg(nglob))
- if(.not. Rate_AND_State) then
- DataXZ%sta => bc%swf%theta
- else
- DataXZ%sta => bc%rsf%theta
- endif
- DataXZ%d1 => bc%d(1,:)
- DataXZ%d2 => bc%d(2,:)
- DataXZ%v1 => bc%v(1,:)
- DataXZ%v2 => bc%v(2,:)
- DataXZ%t1 => bc%t(1,:)
- DataXZ%t2 => bc%t(2,:)
- DataXZ%t3 => bc%t(3,:)
- DataXZ%xcoord => bc%coord(1,:)
- DataXZ%ycoord => bc%coord(2,:)
- DataXZ%zcoord => bc%coord(3,:)
- allocate(DataXZ%tRUP(nglob))
- allocate(DataXZ%tPZ(nglob))
-
- !Percy , setting up initial rupture time null for all faults.
- DataXZ%tRUP = 0e0_CUSTOM_REAL
- DataXZ%tPZ = 0e0_CUSTOM_REAL
-
-end subroutine init_dataXZ
-
-!---------------------------------------------------------------
-subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
-
- type(dataXZ_type), intent(inout) :: dataXZ
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: stg,dold,dnew,dc,vold,vnew
- real(kind=CUSTOM_REAL), intent(in) :: time,dt
-
- integer :: i
-
- ! "stg" : strength .
-
- dataXZ%stg = stg
-
- do i = 1,size(stg)
- ! process zone time = first time when slip = dc (break down process).
- ! with linear time interpolation
- if (dataXZ%tPZ(i)==0e0_CUSTOM_REAL) then
- if (dold(i)<=dc(i) .and. dnew(i) >= dc(i)) then
- dataXZ%tPZ(i) = time-dt*(dnew(i)-dc(i))/(dnew(i)-dold(i))
- endif
- endif
- ! rupture time = first time when slip velocity = vc
- ! with linear time interpolation
- ! vc should be pre-defined as input data .
-
- if (dataXZ%tRUP(i)==0e0_CUSTOM_REAL) then
- if (vold(i)<=V_RUPT .and. vnew(i)>=V_RUPT) dataXZ%tRUP(i)= time-dt*(vnew(i)-V_RUPT)/(vnew(i)-vold(i))
- endif
- enddo
-
- ! To do : add stress criteria (firs time strength is reached).
-
- ! note: the other arrays in dataXZ are pointers to arrays in bc
- ! they do not need to be updated here
-
-end subroutine store_dataXZ
-
-!---------------------------------------------------------------
-subroutine write_dataXZ(dataXZ,itime,iflt)
-
- type(dataXZ_type), intent(in) :: dataXZ
- integer, intent(in) :: itime,iflt
-
- character(len=70) :: filename
-
- write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
- ! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
- ! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
- ! compilers.
-
- ! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
-
- open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
-
- write(IOUT) dataXZ%xcoord
- write(IOUT) dataXZ%ycoord
- write(IOUT) dataXZ%zcoord
- write(IOUT) dataXZ%d1
- write(IOUT) dataXZ%d2
- write(IOUT) dataXZ%v1
- write(IOUT) dataXZ%v2
- write(IOUT) dataXZ%t1
- write(IOUT) dataXZ%t2
- write(IOUT) dataXZ%t3
- write(IOUT) dataXZ%sta
- write(IOUT) dataXZ%stg
- write(IOUT) dataXZ%tRUP
- write(IOUT) dataXZ%tPZ
- close(IOUT)
-
-end subroutine write_dataXZ
-
-
-end module fault_solver
+!=====================================================================
+!
+! s p e c f e m 3 d v e r s i o n 1 . 4
+! ---------------------------------------
+!
+! dimitri komatitsch and jeroen tromp
+! seismological laboratory - california institute of technology
+! (c) california institute of technology september 2006
+!
+! this program is free software; you can redistribute it and/or modify
+! it under the terms of the gnu general public license as published by
+! the free software foundation; either version 2 of the license, or
+! (at your option) any later version.
+!
+! this program is distributed in the hope that it will be useful,
+! but without any warranty; without even the implied warranty of
+! merchantability or fitness for a particular purpose. see the
+! gnu general public license for more details.
+!
+! you should have received a copy of the gnu general public license along
+! with this program; if not, write to the free software foundation, inc.,
+! 51 franklin street, fifth floor, boston, ma 02110-1301 usa.
+!
+!===============================================================================================================
+
+! This module was written by:
+! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
+! Surendra Nadh Somala : Added Heterogenous initial stress capabilities (based on TPV16)
+! Surendra Nadh Somala : Added Rate and State Friction
+
+module fault_solver
+
+ implicit none
+
+ include 'constants.h'
+
+ private
+
+ ! outputs on selected fault nodes at every time step:
+ ! slip, slip velocity, fault stresses
+ type dataT_type
+ integer :: npoin
+ integer, dimension(:), pointer :: iglob ! on-fault global index of output nodes
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1,v1,t1,d2,v2,t2,t3,theta
+ character(len=70), dimension(:), pointer :: name
+ end type dataT_type
+
+
+ ! outputs at selected times for all fault nodes:
+ ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
+ ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
+ ! process zone time = first time when slip = Dc
+ type dataXZ_type
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: stg, sta, d1, d2, v1, v2, &
+ t1, t2, t3, tRUP,tPZ
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord,ycoord,zcoord
+ integer :: npoin
+ end type dataXZ_type
+
+ type swf_type
+ private
+ integer :: kind
+ logical :: healing = .false.
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), &
+ theta=>null(), T=>null(), C=>null()
+ end type swf_type
+
+ type rsf_type
+ private
+ integer :: StateLaw = 1 ! 1=ageing law, 2=slip law
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), &
+ V_init=>null(), &
+ a=>null(), b=>null(), theta=>null(), &
+ T=>null(), C=>null()
+ end type rsf_type
+
+ type asperity_type
+ private
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: Fload=>null()
+ end type asperity_type
+
+ type bc_dynflt_type
+ private
+ integer :: nspec,nglob
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0,T,V,D
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: coord
+ real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: MU,B,invM1,invM2,Z
+ real(kind=CUSTOM_REAL) :: dt
+ integer, dimension(:), pointer :: ibulk1, ibulk2
+ type(swf_type), pointer :: swf => null()
+ type(rsf_type), pointer :: rsf => null()
+ type(asperity_type), pointer :: asp => null()
+ logical :: allow_opening = .false. ! default : do not allow opening
+ type(dataT_type) :: dataT
+ type(dataXZ_type) :: dataXZ
+ end type bc_dynflt_type
+
+ type(bc_dynflt_type), allocatable, save :: faults(:)
+
+ !slip velocity threshold for healing
+ !WARNING: not very robust
+ real(kind=CUSTOM_REAL), save :: V_HEALING
+
+ !slip velocity threshold for definition of rupture front
+ real(kind=CUSTOM_REAL), save :: V_RUPT
+
+ !Number of time steps defined by the user : NTOUT
+ integer, save :: NTOUT,NSNAP
+
+ logical, save :: SIMULATION_TYPE_DYN = .false.
+
+ logical, save :: TPV16 = .false.
+
+ logical, save :: Rate_AND_State = .true.
+
+ real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
+
+ public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
+ SIMULATION_TYPE_DYN
+
+
+contains
+
+!=====================================================================
+! BC_DYNFLT_init initializes dynamic faults
+!
+! prname fault database is read from file prname_fault_db.bin
+! Minv inverse mass matrix
+! dt global time step
+!
+subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt)
+
+ character(len=256), intent(in) :: prname ! 'proc***'
+ real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
+ double precision, intent(in) :: DTglobal
+ integer, intent(in) :: nt
+
+ real(kind=CUSTOM_REAL) :: dt
+ integer :: iflt,ier,dummy_idfault
+ integer :: nbfaults
+ integer :: size_Kelvin_Voigt
+ integer :: SIMULATION_TYPE
+ character(len=256) :: filename
+ integer, parameter :: IIN_PAR =151
+ integer, parameter :: IIN_BIN =170
+
+ NAMELIST / BEGIN_FAULT / dummy_idfault
+
+ dummy_idfault = 0
+
+ open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults',status='old',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File Par_file_faults not found: assume no faults'
+ close(IIN_PAR)
+ return
+ endif
+
+ dt = real(DTglobal)
+ filename = prname(1:len_trim(prname))//'fault_db.bin'
+ open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File ',trim(filename),' not found. Abort'
+ stop
+ endif
+ ! WARNING TO DO: should be an MPI abort
+
+ read(IIN_PAR,*) nbfaults
+ if (nbfaults==0) return
+ ! Reading etas of each fault
+ do iflt = 1,nbfaults
+ read(IIN_PAR,*) ! etas
+ enddo
+ read(IIN_PAR,*) SIMULATION_TYPE
+ if ( SIMULATION_TYPE /= 1 ) then
+ close(IIN_BIN)
+ close(IIN_PAR)
+ return
+ endif
+ SIMULATION_TYPE_DYN = .true.
+ read(IIN_PAR,*) NTOUT
+ read(IIN_PAR,*) NSNAP
+ read(IIN_PAR,*) V_HEALING
+ read(IIN_PAR,*) V_RUPT
+
+ read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
+ allocate( faults(nbfaults) )
+ do iflt=1,nbfaults
+ read(IIN_PAR,nml=BEGIN_FAULT,end=100)
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
+ enddo
+ close(IIN_BIN)
+ close(IIN_PAR)
+
+ filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
+ open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File ',trim(filename),' not found. Abort'
+ stop
+ endif
+ read(IIN_BIN) size_Kelvin_Voigt
+ if (size_Kelvin_Voigt > 0) then
+ allocate(Kelvin_Voigt_eta(size_Kelvin_Voigt))
+ read(IIN_BIN) Kelvin_Voigt_eta
+ endif
+ close(IIN_BIN)
+ return
+100 stop 'Did not find BEGIN_FAULT block #'
+ ! WARNING TO DO: should be an MPI abort
+
+end subroutine BC_DYNFLT_init
+
+!---------------------------------------------------------------------
+
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+
+ type(bc_dynflt_type), intent(inout) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
+ integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: dt
+
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
+ integer, dimension(:,:), allocatable :: ibool1
+ real(kind=CUSTOM_REAL) :: norm
+ real(kind=CUSTOM_REAL) :: S1,S2,S3
+ integer :: n1,n2,n3
+ real(kind=CUSTOM_REAL) :: mus,mud,dc
+ integer :: nmus,nmud,ndc,ij,k,e
+ 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) :: C,T
+ integer :: nC,nForcedRup
+
+
+ integer, parameter :: IIN_NUC =270
+ integer :: ipar
+ integer :: ier
+
+ integer :: relz_num,sub_relz_num
+ integer :: num_cell_str,num_cell_dip
+ real(kind=CUSTOM_REAL) :: siz_str,siz_dip
+ integer :: hypo_cell_str,hypo_cell_dip
+ real(kind=CUSTOM_REAL) :: hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+
+ integer, dimension(:), allocatable :: inp_nx,inp_nz
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
+ dyn_fc,swcd,cohes,tim_forcedRup
+
+ integer :: iLoc
+ real(kind=CUSTOM_REAL) :: minX
+
+
+
+ real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
+ real(kind=CUSTOM_REAL) :: x,z
+ logical :: c1,c2,c3,c4
+ real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
+ integer iglob
+
+ double precision, dimension(:), allocatable :: mu1,mu2,mu3
+
+
+ NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
+ 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
+ allocate( bc%ibulk1(bc%nglob) )
+ allocate( bc%ibulk2(bc%nglob) )
+ allocate( ibool1(NGLLSQUARE,bc%nspec) )
+ allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
+ allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
+
+ allocate(bc%coord(3,(bc%nglob)))
+ read(IIN_BIN) ibool1
+ read(IIN_BIN) jacobian2Dw
+ read(IIN_BIN) normal
+ read(IIN_BIN) bc%ibulk1
+ read(IIN_BIN) bc%ibulk2
+ read(IIN_BIN) bc%coord(1,:)
+ read(IIN_BIN) bc%coord(2,:)
+ read(IIN_BIN) bc%coord(3,:)
+ bc%dt = dt
+
+ allocate( bc%B(bc%nglob) )
+ bc%B = 0e0_CUSTOM_REAL
+ allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
+ nx = 0e0_CUSTOM_REAL
+ ny = 0e0_CUSTOM_REAL
+ nz = 0e0_CUSTOM_REAL
+ do e=1,bc%nspec
+ do ij = 1,NGLLSQUARE
+ k = ibool1(ij,e)
+ nx(k) = nx(k) + normal(1,ij,e)
+ ny(k) = ny(k) + normal(2,ij,e)
+ nz(k) = nz(k) + normal(3,ij,e)
+ bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
+ enddo
+ enddo
+ do k=1,bc%nglob
+ norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
+ nx(k) = nx(k) / norm
+ ny(k) = ny(k) / norm
+ nz(k) = nz(k) / norm
+ enddo
+
+ allocate( bc%R(3,3,bc%nglob) )
+ call compute_R(bc%R,bc%nglob,nx,ny,nz)
+
+ ! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
+ allocate(bc%invM1(bc%nglob))
+ allocate(bc%invM2(bc%nglob))
+ bc%invM1 = Minv(bc%ibulk1)
+ bc%invM2 = Minv(bc%ibulk2)
+
+ ! Fault impedance, Z in : Trac=T_Stick-Z*dV
+ ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
+ ! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
+ ! NOTE: same Bi on both sides, see note above
+ allocate(bc%Z(bc%nglob))
+ bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
+
+ allocate(bc%T(3,bc%nglob))
+ allocate(bc%D(3,bc%nglob))
+ allocate(bc%V(3,bc%nglob))
+ bc%T = 0e0_CUSTOM_REAL
+ bc%D = 0e0_CUSTOM_REAL
+ bc%V = 0e0_CUSTOM_REAL
+
+ ! Set initial fault stresses
+ allocate(bc%T0(3,bc%nglob))
+ S1 = 0e0_CUSTOM_REAL
+ S2 = 0e0_CUSTOM_REAL
+ S3 = 0e0_CUSTOM_REAL
+ n1=0
+ n2=0
+ n3=0
+ read(IIN_PAR, nml=INIT_STRESS)
+ bc%T0(1,:) = S1
+ bc%T0(2,:) = S2
+ bc%T0(3,:) = S3
+
+ call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
+ call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
+ call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
+
+ !WARNING : Quick and dirty free surface condition at z=0
+ ! do k=1,bc%nglob
+ ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
+ ! end do
+
+
+ if (TPV16) then
+
+ allocate(inp_nx(bc%nglob))
+ allocate(inp_nz(bc%nglob))
+ allocate(loc_str(bc%nglob))
+ allocate(loc_dip(bc%nglob))
+ allocate(sigma0(bc%nglob))
+ allocate(tau0_str(bc%nglob))
+ allocate(tau0_dip(bc%nglob))
+ allocate(Rstress_str(bc%nglob))
+ allocate(Rstress_dip(bc%nglob))
+ allocate(static_fc(bc%nglob))
+ allocate(dyn_fc(bc%nglob))
+ allocate(swcd(bc%nglob))
+ allocate(cohes(bc%nglob))
+ allocate(tim_forcedRup(bc%nglob))
+
+ open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
+ read(IIN_NUC,*) relz_num,sub_relz_num
+ read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
+ read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ do ipar=1,bc%nglob
+ read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
+ Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
+ enddo
+ close(IIN_NUC)
+
+ allocate( bc%swf )
+ allocate( bc%swf%mus(bc%nglob) )
+ allocate( bc%swf%mud(bc%nglob) )
+ allocate( bc%swf%Dc(bc%nglob) )
+ allocate( bc%swf%theta(bc%nglob) )
+ allocate( bc%swf%T(bc%nglob) )
+ allocate( bc%swf%C(bc%nglob) )
+
+ minX = minval(bc%coord(1,:))
+
+ do iLoc=1,bc%nglob
+
+ ipar = minloc( (minX+loc_str(:)-bc%coord(1,iLoc))**2 + (-loc_dip(:)-bc%coord(3,iLoc))**2 , 1) !iloc_dip is negative of Z-coord
+ bc%T0(3,iLoc) = -sigma0(ipar)
+ bc%T0(1,iLoc) = tau0_str(ipar)
+ bc%T0(2,iLoc) = tau0_dip(ipar)
+
+ bc%swf%mus(iLoc) = static_fc(ipar)
+ bc%swf%mud(iLoc) = dyn_fc(ipar)
+ bc%swf%Dc(iLoc) = swcd(ipar)
+ bc%swf%C(iLoc) = cohes(ipar)
+ bc%swf%T(iLoc) = tim_forcedRup(ipar)
+ enddo
+
+ endif
+
+
+ ! Set friction parameters and initialize friction variables
+ ! Slip weakening friction
+ if(.not. Rate_AND_State) then
+ allocate( bc%swf )
+ allocate( bc%swf%mus(bc%nglob) )
+ 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)
+
+ mus = 0.6e0_CUSTOM_REAL
+ mud = 0.1e0_CUSTOM_REAL
+ dc = 1e0_CUSTOM_REAL
+ nmus = 0
+ 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))
+ bc%MU = swf_mu(bc%swf)
+
+ ! Rate and state friction
+ else
+ allocate( bc%rsf )
+ allocate( bc%rsf%V0(bc%nglob) )
+ allocate( bc%rsf%f0(bc%nglob) )
+ allocate( bc%rsf%a(bc%nglob) )
+ allocate( bc%rsf%b(bc%nglob) )
+ 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.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
+
+ nV0 =0
+ nf0 =0
+ na =0
+ nb =0
+ nL =0
+ 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
+ bc%rsf%a = a
+ bc%rsf%b = b
+ 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)
+ call init_2d_distribution(bc%rsf%b,bc%coord,IIN_PAR,nb)
+ 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)
+
+
+ W1=15000._CUSTOM_REAL
+ W2=7500._CUSTOM_REAL
+ w=3000._CUSTOM_REAL
+ hypo_z = -7500._CUSTOM_REAL
+ do iglob=1,bc%nglob
+ x=bc%coord(1,iglob)
+ z=bc%coord(3,iglob)
+ c1=abs(x)<W1+w
+ c2=abs(x)>W1
+ c3=abs(z-hypo_z)<W2+w
+ c4=abs(z-hypo_z)>W2
+ if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
+
+ if (c1 .and. c2) then
+ b11 = w/(abs(x)-W1-w)
+ b12 = w/(abs(x)-W1)
+ B1 = 0.5 * (ONE + tanh(b11 + b12))
+ elseif(abs(x)<=W1) then
+ B1 = 1._CUSTOM_REAL
+ else
+ B1 = 0._CUSTOM_REAL
+ endif
+
+ if (c3 .and. c4) then
+ b21 = w/(abs(z-hypo_z)-W2-w)
+ b22 = w/(abs(z-hypo_z)-W2)
+ B2 = 0.5 * (ONE + tanh(b21 + b22))
+ elseif(abs(z-hypo_z)<=W2) then
+ B2 = 1._CUSTOM_REAL
+ else
+ B2 = 0._CUSTOM_REAL
+ endif
+
+
+ bc%rsf%a(iglob) = 0.008 + 0.008 * (ONE - B1*B2)
+ elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
+ bc%rsf%a(iglob) = 0.008
+ else
+ bc%rsf%a(iglob) = 0.016
+ endif
+
+ bc%rsf%theta(iglob) = L/V0 * exp( (bc%rsf%a(iglob) * log(2.0*sinh(-S1/S3/bc%rsf%a(iglob))) - f0 - bc%rsf%a(iglob)*log(V_init/V0) )/b )
+
+
+ enddo
+
+
+
+ allocate(bc%MU(bc%nglob))
+
+ bc%V(1,:) = bc%rsf%V_init
+ allocate( bc%asp )
+ allocate( bc%asp%Fload(bc%nglob) )
+
+ Fload =0.e0_CUSTOM_REAL
+ nFload =0
+
+ read(IIN_PAR, nml=ASP)
+ bc%asp%Fload =Fload
+ call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
+
+ endif
+
+
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc,bc%nglob)
+
+end subroutine init_one_fault
+
+!---------------------------------------------------------------------
+subroutine compute_R(R,nglob,nx,ny,nz)
+
+ integer :: nglob
+ real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
+ real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
+
+ real(kind=CUSTOM_REAL), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
+
+ ! Percy , defining fault directions (in concordance with SCEC conventions) .
+ ! fault coordinates (s,d,n) = (1,2,3)
+ ! s = strike , d = dip , n = n.
+ ! 1 = strike , 2 = dip , 3 = n.
+ norm = sqrt(nx*nx+ny*ny)
+ sx = ny/norm
+ sy = -nx/norm
+ sz = 0.e0_CUSTOM_REAL
+
+ norm = sqrt(sy*sy*nz*nz+sx*sx*nz*nz+(sy*nx-ny*sx)*(nx*sy-ny*sx))
+ dx = -sy*nz/norm
+ dy = sx*nz/norm
+ dz = (sy*nx-ny*sx)/norm
+ !Percy, dz is always dipwards = -1/norm , because (nx*sy-ny*sx)= - 1
+
+ R(1,1,:)=sx
+ R(1,2,:)=sy
+ R(1,3,:)=sz
+ R(2,1,:)=dx
+ R(2,2,:)=dy
+ R(2,3,:)=dz
+ R(3,1,:)=nx
+ R(3,2,:)=ny
+ R(3,3,:)=nz
+
+end subroutine compute_R
+
+!---------------------------------------------------------------------
+! adds a value to a fault parameter inside an area with prescribed shape
+subroutine init_2d_distribution(a,coord,iin,n)
+
+ real(kind=CUSTOM_REAL), intent(inout) :: a(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: iin,n
+
+ 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
+
+ if (n==0) return
+
+ do i=1,n
+ shape = ''
+ xc = 0e0_CUSTOM_REAL
+ yc = 0e0_CUSTOM_REAL
+ zc = 0e0_CUSTOM_REAL
+ r = 0e0_CUSTOM_REAL
+ l = 0e0_CUSTOM_REAL
+ lx = 0e0_CUSTOM_REAL
+ ly = 0e0_CUSTOM_REAL
+ lz = 0e0_CUSTOM_REAL
+ valh = 0e0_CUSTOM_REAL
+
+ read(iin,DIST2D)
+ 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')
+ b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle-taper')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
+ case default
+ stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
+ end select
+
+ where (b /= 0e0_CUSTOM_REAL) a = b
+ enddo
+
+end subroutine init_2d_distribution
+
+!---------------------------------------------------------------------
+elemental function heaviside(x)
+
+ real(kind=CUSTOM_REAL), intent(in) :: x
+ real(kind=CUSTOM_REAL) :: heaviside
+
+ if (x>=0e0_CUSTOM_REAL) then
+ heaviside = 1e0_CUSTOM_REAL
+ else
+ heaviside = 0e0_CUSTOM_REAL
+ endif
+
+end function heaviside
+
+!=====================================================================
+! adds boundary term Bt into Force array for each fault.
+!
+subroutine bc_dynflt_set3d_all(F,Vel,Dis)
+
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
+
+ integer :: iflt
+
+ if (.not. allocated(faults)) return
+ do iflt=1,size(faults)
+ if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
+ enddo
+
+end subroutine bc_dynflt_set3d_all
+
+!---------------------------------------------------------------------
+subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
+
+ use specfem_par, only:it,NSTEP
+
+ real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
+ type(bc_dynflt_type), intent(inout) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
+ integer,intent(in) :: iflt
+
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,tnew
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, theta_new, Vnorm, Vnorm_old, dc
+ real(kind=CUSTOM_REAL) :: half_dt
+
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: ta, Vf,Vf1,tau1,Vf2,tau2,Vfavg
+ integer :: ierr,iNode
+ 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,:))
+
+ ! get predicted values
+ dD = get_jump(bc,D) ! dD_predictor
+ dV = get_jump(bc,V) ! dV_predictor
+ dA = get_weighted_jump(bc,MxA) ! dA_free
+
+ ! rotate to fault frame (tangent,normal)
+ ! component 3 is normal to the fault
+ dD = rotate(bc,dD,1)
+ dV = rotate(bc,dV,1)
+ dA = rotate(bc,dA,1)
+
+ ! T_stick
+ T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
+ T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
+ T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
+
+ !Warning : dirty particular free surface condition z = 0.
+ ! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
+ ! do k=1,bc%nglob
+ ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
+ ! end do
+
+ ! add initial stress
+ T = T + bc%T0
+
+ ! Solve for normal stress (negative is compressive)
+ ! Opening implies free stress
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+
+ if(.not. Rate_AND_State) then ! Update slip weakening friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
+ theta_old = bc%swf%theta
+ call swf_update_state(bc%D,dD,bc%V,bc%swf)
+
+ ! Update friction coeficient
+ bc%MU = swf_mu(bc%swf)
+
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ if (TPV16) then
+ where (bc%swf%T <= it*bc%dt) bc%MU = bc%swf%mud
+ endif
+
+ ! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
+
+ ! Solve for shear stress
+ tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnew = min(tStick,strength)
+ tStick = max(tStick,1e0_CUSTOM_REAL)
+ T(1,:) = tnew * T(1,:)/tStick
+ T(2,:) = tnew * T(2,:)/tStick
+
+ else ! Update rate and state friction:
+
+ ! smooth loading within nucleation patch
+ !WARNING : ad hoc for SCEC benchmark TPV10x
+ TxExt = 0._CUSTOM_REAL
+ TLoad = 1.0_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 = 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)
+
+ tStick = sqrt(T(1,:)**2 + T(2,:)**2)
+ do iNode=1,bc%nglob
+ Vf1(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
+ tau1(iNode) = tStick(iNode) - bc%Z(iNode)*Vf1(iNode)
+ enddo
+
+ ! Updating state variable: 2nd loop
+ Vfavg = 0.5e0_CUSTOM_REAL*(Vf + Vf1)
+ bc%rsf%theta = theta_old
+ call rsf_update_state(Vfavg,bc%dt,bc%rsf)
+
+ ! NR search 2nd loop
+ do iNode=1,bc%nglob
+ Vf2(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
+ tau2(iNode) = tStick(iNode) - bc%Z(iNode)*Vf2(iNode)
+ enddo
+
+ tStick = max(tStick,1e0_CUSTOM_REAL)
+ T(1,:) = tau2 * T(1,:)/tStick
+ T(2,:) = tau2 * T(2,:)/tStick
+
+ endif
+
+ ! Save total tractions
+ bc%T = T
+
+ ! Subtract initial stress
+ T = T - bc%T0
+
+ ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
+ dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
+ dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
+ dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
+
+ ! Update slip and slip rate, in fault frame
+ bc%D = dD
+ bc%V = dV + half_dt*dA
+
+ ! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
+
+ ! Add boundary term B*T to M*a
+ MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
+ MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
+ MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
+
+ MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
+ MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
+ MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
+
+ !-- intermediate storage of outputs --
+ Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ if(.not. Rate_AND_State) then
+ theta_new = bc%swf%theta
+ dc = bc%swf%dc
+ else
+ theta_new = bc%rsf%theta
+ dc = bc%rsf%L
+ endif
+
+ call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
+ Vnorm_old, Vnorm, it*bc%dt,bc%dt)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
+
+ !-- outputs --
+ ! write dataT every NTOUT time step or at the end of simulation
+ if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
+ ! write dataXZ every NSNAP time step
+ if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
+ if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
+
+end subroutine BC_DYNFLT_set3d
+
+!===============================================================
+function get_jump (bc,v) result(dv)
+
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
+ real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
+
+ ! diference between side 2 and side 1 of fault nodes. dv
+ dv(1,:) = v(1,bc%ibulk2)-v(1,bc%ibulk1)
+ dv(2,:) = v(2,bc%ibulk2)-v(2,bc%ibulk1)
+ dv(3,:) = v(3,bc%ibulk2)-v(3,bc%ibulk1)
+
+end function get_jump
+
+!---------------------------------------------------------------------
+function get_weighted_jump (bc,f) result(da)
+
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
+
+ real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
+
+ ! diference between side 2 and side 1 of fault nodes. M-1 * F
+ da(1,:) = bc%invM2*f(1,bc%ibulk2)-bc%invM1*f(1,bc%ibulk1)
+ da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
+ da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
+
+end function get_weighted_jump
+
+!----------------------------------------------------------------------
+function rotate(bc,v,fb) result(vr)
+
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: v(3,bc%nglob)
+ integer, intent(in) :: fb
+ real(kind=CUSTOM_REAL) :: vr(3,bc%nglob)
+
+ ! Percy, tangential direction Vt, equation 7 of Pablo's notes in agreement with SPECFEM3D
+
+ ! forward rotation
+ if (fb==1) then
+ vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(1,2,:)+v(3,:)*bc%R(1,3,:) ! vs
+ vr(2,:) = v(1,:)*bc%R(2,1,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(2,3,:) ! vd
+ vr(3,:) = v(1,:)*bc%R(3,1,:)+v(2,:)*bc%R(3,2,:)+v(3,:)*bc%R(3,3,:) ! vn
+
+ ! backward rotation
+ else
+ vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(2,1,:)+v(3,:)*bc%R(3,1,:) !vx
+ vr(2,:) = v(1,:)*bc%R(1,2,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(3,2,:) !vy
+ vr(3,:) = v(1,:)*bc%R(1,3,:)+v(2,:)*bc%R(2,3,:)+v(3,:)*bc%R(3,3,:) !vz
+
+ endif
+
+end function rotate
+
+
+!=====================================================================
+subroutine swf_update_state(dold,dnew,vold,f)
+
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
+ type(swf_type), intent(inout) :: f
+
+ real(kind=CUSTOM_REAL) :: vnorm
+ integer :: k,npoin
+
+ f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
+
+ if (f%healing) then
+ npoin = size(vold,2)
+ do k=1,npoin
+ vnorm = sqrt(vold(1,k)**2 + vold(2,k)**2)
+ if (vnorm<V_HEALING) f%theta(k) = 0e0_CUSTOM_REAL
+ enddo
+ endif
+end subroutine swf_update_state
+
+!---------------------------------------------------------------------
+function swf_mu(f) result(mu)
+
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ mu = max( mu, f%mud)
+
+end function swf_mu
+
+!=====================================================================
+! Rate and state friction coefficient
+function rsf_mu(f,V) result(mu)
+
+ type(rsf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ real(kind=CUSTOM_REAL) :: mu(size(V))
+ mu = f%a * asinh( V/2.0/f%V0 * exp((f%f0 + f%b*log(f%theta*f%V0/f%L))/f%a ) ) ! Regularized
+
+end function rsf_mu
+
+!---------------------------------------------------------------------
+subroutine funcd(x,fn,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+
+ real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
+ double precision :: arg,fn,df,x
+
+ arg = exp((f0+dble(b)*log(V0*theta/L))/a)/2._CUSTOM_REAL/V0
+ fn = tStick - Z*x - a*Seff*asinh(x*arg)
+ df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
+end subroutine funcd
+
+!---------------------------------------------------------------------
+function rtsafe(funcd,x1,x2,xacc,tStick,Seff,Z,f0,V0,a,b,L,theta)
+
+ integer, parameter :: MAXIT=200
+ real(kind=CUSTOM_REAL) :: x1,x2,xacc
+ EXTERNAL funcd
+ integer :: j
+ !real(kind=CUSTOM_REAL) :: df,dx,dxold,f,fh,fl,temp,xh,xl
+ double precision :: df,dx,dxold,f,fh,fl,temp,xh,xl,rtsafe
+ real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
+
+ call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ if( (fl>0 .and. fh>0) .or. (fl<0 .and. fh<0) ) stop 'root must be bracketed in rtsafe'
+ if(fl==0.) then
+ rtsafe=x2
+ return
+ elseif(fh==0.) then
+ rtsafe=x2
+ return
+ elseif(fl<0) then
+ xl=x1
+ xh=x2
+ else
+ xh=x1
+ xl=x2
+ endif
+
+ rtsafe=0.5d0*(x1+x2)
+ dxold=abs(x2-x1)
+ dx=dxold
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ do j=1,MAXIT
+ if( ((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f)>0 .or. abs(2.*f)>abs(dxold*df) ) then
+ dxold=dx
+ dx=0.5d0*(xh-xl)
+ rtsafe=xl+dx
+ if(xl==rtsafe) return
+ else
+ dxold=dx
+ dx=f/df
+ temp=rtsafe
+ rtsafe=rtsafe-dx
+ if(temp==rtsafe) return
+ endif
+ if(abs(dx)<xacc) return
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ if(f<0.) then
+ xl=rtsafe
+ else
+ xh=rtsafe
+ endif
+ enddo
+ stop 'rtsafe exceeding maximum iterations'
+ return
+
+end function rtsafe
+
+!---------------------------------------------------------------------
+subroutine rsf_update_state(V,dt,f)
+
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ type(rsf_type), intent(inout) :: f
+ real(kind=CUSTOM_REAL), intent(in) :: dt
+
+ real(kind=CUSTOM_REAL) :: vDtL(size(V))
+
+ vDtL = V * dt / f%L
+
+ ! ageing law
+ if (f%StateLaw == 1) then
+ where(vDtL > 1.e-5_CUSTOM_REAL)
+ f%theta = f%theta * exp(-vDtL) + f%L/V * (ONE - exp(-vDtL))
+ elsewhere
+ f%theta = f%theta * exp(-vDtL) + dt*( ONE - HALF*vDtL )
+ endwhere
+
+ ! slip law
+ else
+ f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
+ endif
+
+end subroutine rsf_update_state
+
+!===============================================================
+! OUTPUTS
+subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+ ! NT = total number of time steps
+
+ integer, intent(in) :: nglob,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
+ type (dataT_type), intent(out) :: DataT
+
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
+
+ ! 1. read fault output coordinates from user file,
+ ! 2. define iglob: the fault global index of the node nearest to user
+ ! requested coordinate
+
+ IIN = 251
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ read(IIN,*) np
+ DataT%npoin =0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt==iflt) DataT%npoin = DataT%npoin +1
+ enddo
+ close(IIN)
+
+ if (DataT%npoin == 0) return
+
+ allocate(DataT%iglob(DataT%npoin))
+ allocate(DataT%name(DataT%npoin))
+
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
+ read(IIN,*) np
+ k = 0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt/=iflt) cycle
+ k = k+1
+ DataT%name(k) = tmpname
+ !search nearest node
+ distkeep = huge(distkeep)
+
+ do iglob=1,nglob
+ dist = sqrt((coord(1,iglob)-xtarget)**2 &
+ + (coord(2,iglob)-ytarget)**2 &
+ + (coord(3,iglob)-ztarget)**2)
+ if (dist < distkeep) then
+ distkeep = dist
+ DataT%iglob(k) = iglob
+ endif
+ enddo
+ enddo
+
+ ! 3. allocate arrays and set to zero
+ allocate(DataT%d1(NT,DataT%npoin))
+ allocate(DataT%v1(NT,DataT%npoin))
+ allocate(DataT%t1(NT,DataT%npoin))
+ allocate(DataT%d2(NT,DataT%npoin))
+ allocate(DataT%v2(NT,DataT%npoin))
+ allocate(DataT%t2(NT,DataT%npoin))
+ allocate(DataT%t3(NT,DataT%npoin))
+ allocate(DataT%theta(NT,DataT%npoin))
+ DataT%d1 = 0e0_CUSTOM_REAL
+ DataT%v1 = 0e0_CUSTOM_REAL
+ DataT%t1 = 0e0_CUSTOM_REAL
+ DataT%d2 = 0e0_CUSTOM_REAL
+ DataT%v2 = 0e0_CUSTOM_REAL
+ DataT%t2 = 0e0_CUSTOM_REAL
+ DataT%t3 = 0e0_CUSTOM_REAL
+ DataT%theta = 0e0_CUSTOM_REAL
+
+ close(IIN)
+
+end subroutine init_dataT
+
+!---------------------------------------------------------------
+subroutine store_dataT(dataT,d,v,t,theta,itime)
+
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
+ integer, intent(in) :: itime
+
+ integer :: i,k
+
+ do i=1,dataT%npoin
+ k = dataT%iglob(i)
+ dataT%d1(itime,i) = d(1,k)
+ dataT%d2(itime,i) = d(2,k)
+ dataT%v1(itime,i) = v(1,k)
+ dataT%v2(itime,i) = v(2,k)
+ dataT%t1(itime,i) = t(1,k)
+ dataT%t2(itime,i) = t(2,k)
+ dataT%t3(itime,i) = t(3,k)
+ dataT%theta(itime,i) = theta(k)
+ enddo
+
+end subroutine store_dataT
+
+!-----------------------------------------------------------------
+subroutine write_dataT_all(nt)
+
+ integer, intent(in) :: nt
+
+ integer :: i
+
+ if (.not.allocated(faults)) return
+ do i = 1,size(faults)
+ call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
+ enddo
+
+end subroutine write_dataT_all
+
+!------------------------------------------------------------------------
+subroutine SCEC_write_dataT(dataT,DT,NT)
+
+ type(dataT_type), intent(in) :: dataT
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT
+
+ integer :: i,k,IOUT
+ character :: NTchar*5
+ integer :: today(3), now(3)
+
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+
+ write(NTchar,1) NT
+ NTchar = adjustl(NTchar)
+
+1 format(I5)
+ do i=1,dataT%npoin
+
+ open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
+ write(IOUT,*) "# problem=TPV102"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ write(IOUT,1000) today(2), today(1), today(3), now
+ write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
+ write(IOUT,*) "# code_version=1.1"
+ write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
+ write(IOUT,*) "# time_step=",DT
+ write(IOUT,*) "# location=",trim(dataT%name(i))
+ write(IOUT,*) "# Column #1 = Time (s)"
+ write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
+ write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
+ write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
+ write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
+ write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
+ write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
+ write(IOUT,*) "# Column #8 = normal stress (MPa)"
+ if(Rate_AND_State) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
+ write(IOUT,*) "#"
+ write(IOUT,*) "# The line below lists the names of the data fields:"
+ if(.not. Rate_AND_State) then
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(8(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL
+ enddo
+ else
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
+ enddo
+ endif
+ close(IOUT)
+
+
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+
+
+ enddo
+
+end subroutine SCEC_write_dataT
+
+!-------------------------------------------------------------------------------------------------
+subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
+
+ type(dataXZ_type), intent(in) :: dataXZ
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT,iflt
+
+ integer :: i,IOUT
+ character(len=70) :: filename
+ integer*4 today(3), now(3)
+
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+
+
+ write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
+
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+
+ open(IOUT,file=trim(filename),status='replace')
+ write(IOUT,*) "# problem=TPV102"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ write(IOUT,1000) today(2), today(1), today(3), now
+ write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
+ write(IOUT,*) "# code_version=1.1"
+ write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
+ write(IOUT,*) "# Column #1 = horizontal coordinate, distance along strike (m)"
+ write(IOUT,*) "# Column #2 = vertical coordinate, distance down-dip (m)"
+ write(IOUT,*) "# Column #3 = rupture time (s)"
+ write(IOUT,*) "# "
+ write(IOUT,*) "j k t"
+ do i = 1,size(dataXZ%tRUP)
+ write(IOUT,'(3(E15.7))') dataXZ%xcoord(i), -dataXZ%zcoord(i), dataXZ%tRUP(i)
+ end do
+
+ close(IOUT)
+
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+
+end subroutine SCEC_Write_RuptureTime
+
+!-------------------------------------------------------------------------------------------------
+subroutine init_dataXZ(DataXZ,bc,nglob)
+
+ type(dataXZ_type), intent(inout) :: DataXZ
+ type(bc_dynflt_type) :: bc
+ integer, intent(in) :: nglob
+
+ allocate(DataXZ%stg(nglob))
+ if(.not. Rate_AND_State) then
+ DataXZ%sta => bc%swf%theta
+ else
+ DataXZ%sta => bc%rsf%theta
+ endif
+ DataXZ%d1 => bc%d(1,:)
+ DataXZ%d2 => bc%d(2,:)
+ DataXZ%v1 => bc%v(1,:)
+ DataXZ%v2 => bc%v(2,:)
+ DataXZ%t1 => bc%t(1,:)
+ DataXZ%t2 => bc%t(2,:)
+ DataXZ%t3 => bc%t(3,:)
+ DataXZ%xcoord => bc%coord(1,:)
+ DataXZ%ycoord => bc%coord(2,:)
+ DataXZ%zcoord => bc%coord(3,:)
+ allocate(DataXZ%tRUP(nglob))
+ allocate(DataXZ%tPZ(nglob))
+
+ !Percy , setting up initial rupture time null for all faults.
+ DataXZ%tRUP = 0e0_CUSTOM_REAL
+ DataXZ%tPZ = 0e0_CUSTOM_REAL
+
+end subroutine init_dataXZ
+
+!---------------------------------------------------------------
+subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
+
+ type(dataXZ_type), intent(inout) :: dataXZ
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: stg,dold,dnew,dc,vold,vnew
+ real(kind=CUSTOM_REAL), intent(in) :: time,dt
+
+ integer :: i
+
+ ! "stg" : strength .
+
+ dataXZ%stg = stg
+
+ do i = 1,size(stg)
+ ! process zone time = first time when slip = dc (break down process).
+ ! with linear time interpolation
+ if (dataXZ%tPZ(i)==0e0_CUSTOM_REAL) then
+ if (dold(i)<=dc(i) .and. dnew(i) >= dc(i)) then
+ dataXZ%tPZ(i) = time-dt*(dnew(i)-dc(i))/(dnew(i)-dold(i))
+ endif
+ endif
+ ! rupture time = first time when slip velocity = vc
+ ! with linear time interpolation
+ ! vc should be pre-defined as input data .
+
+ if (dataXZ%tRUP(i)==0e0_CUSTOM_REAL) then
+ if (vold(i)<=V_RUPT .and. vnew(i)>=V_RUPT) dataXZ%tRUP(i)= time-dt*(vnew(i)-V_RUPT)/(vnew(i)-vold(i))
+ endif
+ enddo
+
+ ! To do : add stress criteria (firs time strength is reached).
+
+ ! note: the other arrays in dataXZ are pointers to arrays in bc
+ ! they do not need to be updated here
+
+end subroutine store_dataXZ
+
+!---------------------------------------------------------------
+subroutine write_dataXZ(dataXZ,itime,iflt)
+
+ type(dataXZ_type), intent(in) :: dataXZ
+ integer, intent(in) :: itime,iflt
+
+ character(len=70) :: filename
+
+ write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
+ ! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
+ ! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
+ ! compilers.
+
+ ! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
+
+ open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
+
+ write(IOUT) dataXZ%xcoord
+ write(IOUT) dataXZ%ycoord
+ write(IOUT) dataXZ%zcoord
+ write(IOUT) dataXZ%d1
+ write(IOUT) dataXZ%d2
+ write(IOUT) dataXZ%v1
+ write(IOUT) dataXZ%v2
+ write(IOUT) dataXZ%t1
+ write(IOUT) dataXZ%t2
+ write(IOUT) dataXZ%t3
+ write(IOUT) dataXZ%sta
+ write(IOUT) dataXZ%stg
+ write(IOUT) dataXZ%tRUP
+ write(IOUT) dataXZ%tPZ
+ close(IOUT)
+
+end subroutine write_dataXZ
+
+
+end module fault_solver
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-04-12 18:49:24 UTC (rev 19939)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-04-12 21:26:05 UTC (rev 19940)
@@ -223,9 +223,9 @@
! NOTE: same Bi on both sides, see note above
allocate(bc%Z(bc%nglob))
bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
- ! WARNING: In some fault edge nodes (non-split) Minv is assembled across the fault
+ ! WARNING: In non-split nodes at fault edges Minv is assembled across the fault
! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
- ! In a symmetric mesh (M1=M2) Z will be twice its intended value
+ ! For instance, in a symmetric mesh (M1=M2) Z is twice its intended value
allocate(bc%T(3,bc%nglob))
allocate(bc%D(3,bc%nglob))
@@ -306,6 +306,12 @@
end subroutine BC_KINFLT_set_all
!---------------------------------------------------------------------
+!
+!NOTE: On non-split nodes at fault edges, dD=dV=dA=0 but bc%T is corrupted.
+! That does not affect computations: the net contribution of B*T is =0.
+! However, the output T in these nodes should be ignored.
+! It is =0 if the user sets bc%V=0 there in the input slip rates.
+!
subroutine BC_KINFLT_set_single(bc,MxA,V,D,iflt)
use specfem_par, only:it,NSTEP
@@ -326,7 +332,6 @@
dD = get_jump(bc,D) ! dD_predictor
dV = get_jump(bc,V) ! dV_predictor
dA = get_weighted_jump(bc,MxA) ! dA_free
- ! NOTE: In non-split nodes at fault edges, dD=dV=dA=0
! rotate to fault frame (tangent,normal)
! component 3 is normal to the fault
@@ -374,7 +379,7 @@
t2 = bc%kin_it * bc%kin_dt
! Kinematic velocity_rate
- ! bc%V : Imposed apriori and read from slip rate snapshots (from time reversal)
+ ! bc%V : Imposed a priori and read from slip rate snapshots (from time reversal)
! Linear interpolation between consecutive kinematic time steps.
! V will be given at each time step.
bc%V(1,:) = ( (t2 - time)*bc%v_kin_t1(1,:) + (time - t1)*bc%v_kin_t2(1,:) )/ bc%kin_dt
@@ -393,8 +398,6 @@
! Save tractions
bc%T = T
- !NOTE: On non-split nodes at fault edges, bc%T is corrupted (non zero)
- ! unless the user makes sure that bc%V=0 there.
! Update slip in fault frame
bc%D = dD
@@ -410,7 +413,6 @@
MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
- !NOTE: On non-split nodes at fault edges, the net contribution of B*T is =0
!-- intermediate storage of outputs --
call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
@@ -449,8 +451,8 @@
da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
- ! NOTE: In some fault edge nodes (non-split) Minv and f are assembled across the fault.
- ! Hence, f1=f2, invM1=invM2=1/(M1+M2) instead of invMi=1/Mi, and da=0.
+ ! NOTE: In non-split nodes at fault edges da=0
+ ! (invM and f are assembled across the fault, hence f1=f2 and invM1=invM2)
end function get_weighted_jump
More information about the CIG-COMMITS
mailing list