[cig-commits] r19862 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Sat Mar 24 05:04:54 PDT 2012
Author: surendra
Date: 2012-03-24 05:04:53 -0700 (Sat, 24 Mar 2012)
New Revision: 19862
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Finished rate & state friction
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-24 02:37:29 UTC (rev 19861)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-24 12:04:53 UTC (rev 19862)
@@ -36,66 +36,65 @@
private
- ! outputs on selected fault nodes at every time step:
- ! slip, slip velocity, fault stresses
+ ! 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
- character(len=70), dimension(:), pointer :: name
+ 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
+ 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
+
+ ! 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
+ 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()
+ 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 rs_type
- private
- logical :: healing = .false.
- real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), theta_init=>null(), V_init=>null(), a=>null(), b=>null(), theta=>null(), T=>null(), C=>null()
+ private
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), theta_init=>null(), V_init=>null(), a=>null(), b=>null(), theta=>null(), T=>null(), C=>null()
end type rs_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(rs_type), pointer :: rs => null()
- logical :: allow_opening = .false. ! default : do not allow opening
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ
+ 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(rs_type), pointer :: rs => 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
+ !slip velocity threshold for healing
+ !WARNING: not very robust
real(kind=CUSTOM_REAL), save :: V_HEALING
- !slip velocity threshold for definition of rupture front
+ !slip velocity threshold for definition of rupture front
real(kind=CUSTOM_REAL), save :: V_RUPT
- !Number of time steps defined by the user : NTOUT
+ !Number of time steps defined by the user : NTOUT
integer, save :: NTOUT,NSNAP
logical, save :: SIMULATION_TYPE_DYN = .false.
@@ -105,534 +104,533 @@
logical, save :: Rate_AND_State = .false.
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
-
+
public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
- SIMULATION_TYPE_DYN
+ 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)
+ !=====================================================================
+ ! 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
+ 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
+ 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
+ NAMELIST / BEGIN_FAULT / dummy_idfault
- dummy_idfault = 0
+ 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
+ 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
+ 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)
+ if(.not.TPV16) then
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
+ else
+ call init_one_fault_TPV16(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
+ endif
+ enddo
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)
- if(.not.TPV16) then
- call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
- else
- call init_one_fault_TPV16(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
+
+ 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
- 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
+ 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
+ ! WARNING TO DO: should be an MPI abort
-end subroutine BC_DYNFLT_init
+ end subroutine BC_DYNFLT_init
-!---------------------------------------------------------------------
+ !---------------------------------------------------------------------
-subroutine init_one_fault_TPV16(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+ subroutine init_one_fault_TPV16(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
+ 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), 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
- integer, parameter :: IIN_NUC =270
- integer :: ipar
- integer :: ier
+ 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 :: 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, 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, dimension(:), allocatable :: iLoc
+ !integer, dimension(:), allocatable :: iLoc
- integer :: iLoc
+ integer :: iLoc
- NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
+ NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
+ NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
- 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))
+ 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%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)
+ 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
- enddo
-!JPA assemble bc%B across processors
-!JPA 1. put bc%B in the vector 'accel', which is defined in every GLL node including in the bulk (contains zeros off the fault)
-!JPA 2. call assemble_MPI_vector_* (or maybe assemble_MPI_scalar_*)
-!JPA do the same for bc%n
+ !JPA assemble bc%B across processors
+ !JPA 1. put bc%B in the vector 'accel', which is defined in every GLL node including in the bulk (contains zeros off the fault)
+ !JPA 2. call assemble_MPI_vector_* (or maybe assemble_MPI_scalar_*)
+ !JPA do the same for bc%n
- 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)
+ 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
-! 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)
+ allocate( bc%R(3,3,bc%nglob) )
+ call compute_R(bc%R,bc%nglob,nx,ny,nz)
-! 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
+ ! 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(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))
+ 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
- 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(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))
-print*,minval(inp_nx),maxval(inp_nx)
-print*,minval(inp_nz),maxval(inp_nz)
-print*,minval(loc_str),maxval(loc_str)
-print*,minval(loc_dip),maxval(loc_dip)
-print*,minval(sigma0),maxval(sigma0)
-print*,minval(tau0_str),maxval(tau0_str)
-print*,minval(tau0_dip),maxval(tau0_dip)
-print*,minval(Rstress_str),maxval(Rstress_str)
-print*,minval(Rstress_dip),maxval(Rstress_dip)
-print*,minval(static_fc),maxval(static_fc)
-print*,minval(dyn_fc),maxval(dyn_fc)
-print*,minval(swcd),maxval(swcd)
-print*,minval(cohes),maxval(cohes)
-print*,minval(tim_forcedRup),maxval(tim_forcedRup)
+ 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) )
-! do ipar=1,bc%nglob
-! iLoc = minloc( (loc_str(ipar)-bc%coord(1,:))**2 + (-loc_dip(ipar)-bc%coord(3,:))**2 , 1) !iloc_dip is negative of Z-coord
-! bc%T0(1,iLoc) = sigma0(ipar)
-! bc%T0(2,iLoc) = tau0_str(ipar)
-! bc%T0(3,iLoc) = tau0_dip(ipar)
+ print*,minval(inp_nx),maxval(inp_nx)
+ print*,minval(inp_nz),maxval(inp_nz)
+ print*,minval(loc_str),maxval(loc_str)
+ print*,minval(loc_dip),maxval(loc_dip)
+ print*,minval(sigma0),maxval(sigma0)
+ print*,minval(tau0_str),maxval(tau0_str)
+ print*,minval(tau0_dip),maxval(tau0_dip)
+ print*,minval(Rstress_str),maxval(Rstress_str)
+ print*,minval(Rstress_dip),maxval(Rstress_dip)
+ print*,minval(static_fc),maxval(static_fc)
+ print*,minval(dyn_fc),maxval(dyn_fc)
+ print*,minval(swcd),maxval(swcd)
+ print*,minval(cohes),maxval(cohes)
+ print*,minval(tim_forcedRup),maxval(tim_forcedRup)
-! 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
+ 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) )
- do iLoc=1,bc%nglob
-
- ipar = minloc( (-24000.0+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)
+ ! do ipar=1,bc%nglob
+ ! iLoc = minloc( (loc_str(ipar)-bc%coord(1,:))**2 + (-loc_dip(ipar)-bc%coord(3,:))**2 , 1) !iloc_dip is negative of Z-coord
+ ! bc%T0(1,iLoc) = sigma0(ipar)
+ ! bc%T0(2,iLoc) = tau0_str(ipar)
+ ! bc%T0(3,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
+ ! 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
-print*,' '
-print*,minval(bc%T0(1,:)),maxval(bc%T0(1,:))
-print*,minval(bc%T0(2,:)),maxval(bc%T0(2,:))
-print*,minval(bc%T0(3,:)),maxval(bc%T0(3,:))
-print*,minval(bc%swf%mus),maxval(bc%swf%mus)
-print*,minval(bc%swf%mud),maxval(bc%swf%mud)
-print*,minval(bc%swf%Dc),maxval(bc%swf%Dc)
-print*,minval(bc%swf%C),maxval(bc%swf%C)
-print*,minval(bc%swf%T),maxval(bc%swf%T)
+ do iLoc=1,bc%nglob
- ! WARNING: if V_HEALING is negative we turn off healing
- bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+ ipar = minloc( (-24000.0+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%theta = 0e0_CUSTOM_REAL
- allocate(bc%MU(bc%nglob))
- bc%MU = swf_mu_TPV16(bc%swf,0.0,bc%nglob)
+ 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
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc,bc%nglob)
+ print*,' '
+ print*,minval(bc%T0(1,:)),maxval(bc%T0(1,:))
+ print*,minval(bc%T0(2,:)),maxval(bc%T0(2,:))
+ print*,minval(bc%T0(3,:)),maxval(bc%T0(3,:))
+ print*,minval(bc%swf%mus),maxval(bc%swf%mus)
+ print*,minval(bc%swf%mud),maxval(bc%swf%mud)
+ print*,minval(bc%swf%Dc),maxval(bc%swf%Dc)
+ print*,minval(bc%swf%C),maxval(bc%swf%C)
+ print*,minval(bc%swf%T),maxval(bc%swf%T)
- bc%DataXZ%t1 => bc%T0(1,:)
- bc%DataXZ%t2 => bc%T0(2,:)
- bc%DataXZ%t3 => bc%T0(3,:)
- call write_dataXZ(bc%dataXZ,0,iflt)
- bc%DataXZ%t1 => bc%t(1,:)
- bc%DataXZ%t2 => bc%t(2,:)
- bc%DataXZ%t3 => bc%t(3,:)
+ ! WARNING: if V_HEALING is negative we turn off healing
+ bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+ bc%swf%theta = 0e0_CUSTOM_REAL
+ allocate(bc%MU(bc%nglob))
+ bc%MU = swf_mu_TPV16(bc%swf,0.0,bc%nglob)
-end subroutine init_one_fault_TPV16
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc,bc%nglob)
-!---------------------------------------------------------------------
+ bc%DataXZ%t1 => bc%T0(1,:)
+ bc%DataXZ%t2 => bc%T0(2,:)
+ bc%DataXZ%t3 => bc%T0(3,:)
+ call write_dataXZ(bc%dataXZ,0,iflt)
+ bc%DataXZ%t1 => bc%t(1,:)
+ bc%DataXZ%t2 => bc%t(2,:)
+ bc%DataXZ%t3 => bc%t(3,:)
-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
+ end subroutine init_one_fault_TPV16
- NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
- NAMELIST / RS / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init
+ !---------------------------------------------------------------------
- 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)
+ 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
+
+ NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
+ NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
+ NAMELIST / RS / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init
+
+ 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
- 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
+ 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)
+ 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)
+ ! 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 ))
+ ! 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
+ 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
+ ! 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)
+ 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
+ !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(.not. Rate_AND_State) then
-! Set friction parameters and initialize friction variables
- 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) )
- ! WARNING: if V_HEALING is negative we turn off healing
- bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+ if(.not. Rate_AND_State) then
+ ! Set friction parameters and initialize friction variables
+ 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) )
+ ! 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
+ mus = 0.6e0_CUSTOM_REAL
+ mud = 0.1e0_CUSTOM_REAL
+ dc = 1e0_CUSTOM_REAL
+ nmus = 0
+ nmud = 0
+ ndc = 0
- read(IIN_PAR, nml=SWF)
- bc%swf%mus = mus
- bc%swf%mud = mud
- bc%swf%Dc = dc
- 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)
+ read(IIN_PAR, nml=SWF)
+ bc%swf%mus = mus
+ bc%swf%mud = mud
+ bc%swf%Dc = dc
+ 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)
- bc%swf%theta = 0e0_CUSTOM_REAL
- allocate(bc%MU(bc%nglob))
- bc%MU = swf_mu(bc%swf)
+ bc%swf%theta = 0e0_CUSTOM_REAL
+ allocate(bc%MU(bc%nglob))
+ bc%MU = swf_mu(bc%swf)
-else !Rate AND State (Ageing Law)
-! Set friction parameters and initialize friction variables
- allocate( bc%rs )
- allocate( bc%rs%V0(bc%nglob) )
- allocate( bc%rs%f0(bc%nglob) )
- allocate( bc%rs%a(bc%nglob) )
- allocate( bc%rs%b(bc%nglob) )
- allocate( bc%rs%L(bc%nglob) )
- allocate( bc%rs%V_init(bc%nglob) )
- allocate( bc%rs%theta_init(bc%nglob) )
- ! WARNING: if V_HEALING is negative we turn off healing
- bc%rs%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+ else !Rate AND State (Ageing Law)
+ ! Set friction parameters and initialize friction variables
+ allocate( bc%rs )
+ allocate( bc%rs%V0(bc%nglob) )
+ allocate( bc%rs%f0(bc%nglob) )
+ allocate( bc%rs%a(bc%nglob) )
+ allocate( bc%rs%b(bc%nglob) )
+ allocate( bc%rs%L(bc%nglob) )
+ allocate( bc%rs%V_init(bc%nglob) )
+ allocate( bc%rs%theta_init(bc%nglob) )
+ ! WARNING: if V_HEALING is negative we turn off healing
-V0 =1.e-6_CUSTOM_REAL
-f0 =0.6_CUSTOM_REAL
-a =0.0080_CUSTOM_REAL
-b =0.0120_CUSTOM_REAL
-L =0.0135_CUSTOM_REAL
-V_init =1.e-12_CUSTOM_REAL
-theta_init =1.084207680000000e+09_CUSTOM_REAL
+ V0 =1.e-6_CUSTOM_REAL
+ f0 =0.6_CUSTOM_REAL
+ a =0.0080_CUSTOM_REAL
+ b =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
+ nV0 =0
+ nf0 =0
+ na =0
+ nb =0
+ nL =0
+ nV_init =0
+ ntheta_init =0
- read(IIN_PAR, nml=RS)
- bc%rs%V0 = V0
- bc%rs%f0 = f0
- bc%rs%a = a
- bc%rs%b = b
- bc%rs%L = L
- bc%rs%V_init = V_init
- bc%rs%theta_init = theta_init
- call init_2d_distribution(bc%rs%V0,bc%coord,IIN_PAR,nV0)
- call init_2d_distribution(bc%rs%f0,bc%coord,IIN_PAR,nf0)
- call init_2d_distribution(bc%rs%a,bc%coord,IIN_PAR,na)
- call init_2d_distribution(bc%rs%b,bc%coord,IIN_PAR,nb)
- call init_2d_distribution(bc%rs%L,bc%coord,IIN_PAR,nL)
- call init_2d_distribution(bc%rs%V_init,bc%coord,IIN_PAR,nV_init)
- call init_2d_distribution(bc%rs%theta_init,bc%coord,IIN_PAR,ntheta_init)
+ read(IIN_PAR, nml=RS)
+ bc%rs%V0 = V0
+ bc%rs%f0 = f0
+ bc%rs%a = a
+ bc%rs%b = b
+ bc%rs%L = L
+ bc%rs%V_init = V_init
+ bc%rs%theta_init = theta_init
+ call init_2d_distribution(bc%rs%V0,bc%coord,IIN_PAR,nV0)
+ call init_2d_distribution(bc%rs%f0,bc%coord,IIN_PAR,nf0)
+ call init_2d_distribution(bc%rs%a,bc%coord,IIN_PAR,na)
+ call init_2d_distribution(bc%rs%b,bc%coord,IIN_PAR,nb)
+ call init_2d_distribution(bc%rs%L,bc%coord,IIN_PAR,nL)
+ call init_2d_distribution(bc%rs%V_init,bc%coord,IIN_PAR,nV_init)
+ call init_2d_distribution(bc%rs%theta_init,bc%coord,IIN_PAR,ntheta_init)
- bc%rs%theta = 0e0_CUSTOM_REAL
- allocate(bc%MU(bc%nglob))
- bc%MU = rs_mu(bc%rs)
-endif
+ bc%rs%theta = bc%rs%theta_init
+ allocate(bc%MU(bc%nglob)) !Not necessary?
+ !bc%MU = rs_mu(bc%rs)
+ endif
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc,bc%nglob)
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc,bc%nglob)
-end subroutine init_one_fault
+ 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
+ !---------------------------------------------------------------------
+ subroutine compute_R(R,nglob,nx,ny,nz)
- real(kind=CUSTOM_REAL), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
+ integer :: nglob
+ real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
+ real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
-! 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.
+ 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
@@ -642,7 +640,7 @@
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
+ !Percy, dz is always dipwards = -1/norm , because (nx*sy-ny*sx)= - 1
R(1,1,:)=sx
R(1,2,:)=sy
@@ -653,801 +651,938 @@
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)
+ end subroutine compute_R
- real(kind=CUSTOM_REAL), intent(inout) :: a(:)
- real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
- integer, intent(in) :: iin,n
+ !---------------------------------------------------------------------
+ ! 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) :: b(size(a))
- character(len=20) :: shape
- real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
- integer :: i
+ real(kind=CUSTOM_REAL), intent(inout) :: a(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: iin,n
- NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
+ 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
+ integer :: i
- 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
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
- 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 ('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
+ if (n==0) return
-!---------------------------------------------------------------------
-elemental function heaviside(x)
+ 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
- real(kind=CUSTOM_REAL), intent(in) :: x
- real(kind=CUSTOM_REAL) :: heaviside
+ 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 ('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
- if (x>=0e0_CUSTOM_REAL) then
- heaviside = 1e0_CUSTOM_REAL
- else
- heaviside = 0e0_CUSTOM_REAL
- endif
+ where (b /= 0e0_CUSTOM_REAL) a = b
+ enddo
-end function heaviside
+ end subroutine init_2d_distribution
-!=====================================================================
-! adds boundary term Bt into Force array for each fault.
-!
-subroutine bc_dynflt_set3d_all(F,Vel,Dis)
+ !---------------------------------------------------------------------
+ elemental function heaviside(x)
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
- real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
+ real(kind=CUSTOM_REAL), intent(in) :: x
+ real(kind=CUSTOM_REAL) :: heaviside
- integer :: iflt
-
- if (.not. allocated(faults)) return
- do iflt=1,size(faults)
- if(.not.TPV16) then
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
+ if (x>=0e0_CUSTOM_REAL) then
+ heaviside = 1e0_CUSTOM_REAL
else
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d_TPV16(faults(iflt),F,Vel,Dis,iflt)
+ heaviside = 0e0_CUSTOM_REAL
endif
- enddo
-
-end subroutine bc_dynflt_set3d_all
-!---------------------------------------------------------------------
-subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
-
- use specfem_par, only:it,NSTEP
+ end function heaviside
- 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
+ !=====================================================================
+ ! adds boundary term Bt into Force array for each fault.
+ !
+ subroutine bc_dynflt_set3d_all(F,Vel,Dis)
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
- real(kind=CUSTOM_REAL) :: half_dt
-! integer :: k
-
- half_dt = 0.5e0_CUSTOM_REAL*bc%dt
- theta_old = bc%swf%theta
- Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
-! 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
+ integer :: iflt
-! 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)
+ if (.not. allocated(faults)) return
+ do iflt=1,size(faults)
+ if(.not.TPV16) then
+ if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
+ else
+ if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d_TPV16(faults(iflt),F,Vel,Dis,iflt)
+ endif
+ enddo
+ end subroutine bc_dynflt_set3d_all
-! 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,:) )
+ !---------------------------------------------------------------------
+ subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
-!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)
+ use specfem_par, only:it,NSTEP
-! Update slip weakening friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
- call swf_update_state(bc%D,dD,bc%V,bc%swf)
+ 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
- ! Update friction coeficient
- bc%MU = swf_mu_TPV16(bc%swf,it*bc%dt,bc%nglob)
-
-! combined with time-weakening for nucleation
-! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
+ real(kind=CUSTOM_REAL) :: half_dt
+ ! integer :: k
-! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
-
-! Solve for shear stress
- tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
- tnorm = max(tnorm,1e0_CUSTOM_REAL)
- t1 = T(1,:)/tnorm
- t2 = T(2,:)/tnorm
-
- tnew = min(tnorm,strength)
- T(1,:) = tnew * t1
- T(2,:) = tnew * t2
+ half_dt = 0.5e0_CUSTOM_REAL*bc%dt
+ theta_old = bc%swf%theta
+ Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
-! Save total tractions
- bc%T = T
+ ! 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
-! Subtract initial stress
- T = T - bc%T0
+ ! 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)
-! 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
+ ! 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,:) )
-! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
+ !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 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,:)
+ ! add initial stress
+ T = T + bc%T0
- 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,:)
+ ! Solve for normal stress (negative is compressive)
+ ! Opening implies free stress
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+ ! Update slip weakening friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
+ call swf_update_state(bc%D,dD,bc%V,bc%swf)
-!-- intermediate storage of outputs --
- Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
- Vnorm_old, Vnorm, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
+ ! Update friction coeficient
+ bc%MU = swf_mu_TPV16(bc%swf,it*bc%dt,bc%nglob)
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
-!-- 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)
+ ! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
-end subroutine BC_DYNFLT_set3d_TPV16
+ ! Solve for shear stress
+ tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnorm = max(tnorm,1e0_CUSTOM_REAL)
+ t1 = T(1,:)/tnorm
+ t2 = T(2,:)/tnorm
-!---------------------------------------------------------------------
-subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
-
- use specfem_par, only:it,NSTEP
+ tnew = min(tnorm,strength)
+ T(1,:) = tnew * t1
+ T(2,:) = tnew * t2
- 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
+ ! Save total tractions
+ bc%T = T
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
- real(kind=CUSTOM_REAL) :: half_dt
-
- half_dt = 0.5e0_CUSTOM_REAL*bc%dt
- theta_old = bc%swf%theta
- Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ ! Subtract initial stress
+ T = T - bc%T0
- ! 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
+ ! 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)
- ! 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)
+ ! Update slip and slip rate, in fault frame
+ bc%D = dD
+ bc%V = dV + half_dt*dA
- ! 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,:) )
+ ! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
-!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 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,:)
- ! 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)
+ 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,:)
-if(.not. Rate_AND_State) then ! Update slip weakening friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
- call swf_update_state(bc%D,dD,bc%V,bc%swf)
- ! Update friction coeficient
- bc%MU = swf_mu(bc%swf)
-else ! Update rate and state friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
- call rs_update_state(bc%V,bc%dt,bc%rs)
+ !-- intermediate storage of outputs --
+ Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
+ Vnorm_old, Vnorm, it*bc%dt,bc%dt)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
- ! Update friction coeficient
- bc%MU = rs_mu(bc%rs,bc%V)
-endif
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ !-- 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)
- ! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
+ end subroutine BC_DYNFLT_set3d_TPV16
- ! Solve for shear stress
- tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
- tnorm = max(tnorm,1e0_CUSTOM_REAL)
- t1 = T(1,:)/tnorm
- t2 = T(2,:)/tnorm
+ !---------------------------------------------------------------------
+ subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
- tnew = min(tnorm,strength)
- T(1,:) = tnew * t1
- T(2,:) = tnew * t2
+ use specfem_par, only:it,NSTEP
- ! Save total tractions
- bc%T = T
+ 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
- ! Subtract initial stress
- T = T - bc%T0
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
+ real(kind=CUSTOM_REAL) :: half_dt
- ! 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
+ integer :: ierr
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,ta,tau_ratio
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: Vf,Vf1,tau1
+ real(kind=CUSTOM_REAL) :: ONE
+ logical, dimension(bc%nglob) :: ierr
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
+ ONE = 1.e0_CUSTOM_REAL
+ half_dt = 0.5e0_CUSTOM_REAL*bc%dt
+ theta_old = bc%swf%theta
+ Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- ! 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,:)
+ ! 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
- 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,:)
+ ! 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)
- !-- intermediate storage of outputs --
- Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
- Vnorm_old, Vnorm, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
+ ! 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,:) )
- !-- 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)
+ !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
-end subroutine BC_DYNFLT_set3d
+ ! add initial stress
+ T = T + bc%T0
-!===============================================================
-function get_jump (bc,v) result(dv)
+ ! Solve for normal stress (negative is compressive)
+ ! Opening implies free stress
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
- real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
+ if(.not. Rate_AND_State) then ! Update slip weakening friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
+ call swf_update_state(bc%D,dD,bc%V,bc%swf)
- ! 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
+ ! Update friction coeficient
+ bc%MU = swf_mu(bc%swf)
-!---------------------------------------------------------------------
-function get_weighted_jump (bc,f) result(da)
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
- real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
+ ! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
- ! 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
+ ! Solve for shear stress
+ tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnorm = max(tnorm,1e0_CUSTOM_REAL)
+ t1 = T(1,:)/tnorm
+ t2 = T(2,:)/tnorm
-!----------------------------------------------------------------------
-function rotate(bc,v,fb) result(vr)
+ tnew = min(tnorm,strength)
+ T(1,:) = tnew * t1
+ T(2,:) = tnew * t2
- 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
+ else ! Update rate and state friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
+ Vf = sqrt(bc%V(1,:)**2 + bc%V(2,:)**2)
+ call rs_update_state(Vf,bc%dt,bc%rs)
- ! 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
+ tau_ratio = T(2,:)/T(1,:)
+ tStick = sqrt(T(1,:)**2 + T(2,:)**2)
+ ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
+ call NRsearch(Vf1,tau1,bc%rs%fo,bc%rs%Vo,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr)
- endif
+ ! If Vf1 or tau1 is negative, set Vf1 = 0 and tau1 = TauStick
+ where (Vf1 < 0._CUSTOM_REAL .or. tau1 < 0._CUSTOM_REAL .or. ierr == 1)
+ tau1 = tStick
+ Vf1 = 0._CUSTOM_REAL
+ endwhere
-end function rotate
+ if (any(ierr) == 1) then
+ print*, 'NR search 1st loop failed!'
+ print*, 'iteration = ', it
+ endif
+ ! Updating state variable: 2nd loop
+ Vfavg = 0.5e0_CUSTOM_REAL**abs(Vf + Vf1)
+ call rs_update_state(Vfavg,bc%dt,bc%rs)
-!=====================================================================
-subroutine swf_update_state(dold,dnew,vold,f)
+ ! NR search 2nd loop
+ call NRsearch(Vf2,tau2,bc%rs%fo,bc%rs%Vo,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr)
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
- type(swf_type), intent(inout) :: f
+ where (Vf2 < 0._CUSTOM_REAL .or. tau2 < 0._CUSTOM_REAL .or. ierr == 1)
+ tau2 = tStick
+ Vf2 = 0._CUSTOM_REAL
+ endwhere
- real(kind=CUSTOM_REAL) :: vnorm
- integer :: k,npoin
+ ! back to the relative shear stresses
+ where(T(1,:) > -VERYSMALLVAL)
+ T(1,:) = tau2/sqrt(ONE+tau_ratio**2) - bc%T0(1,:)
+ elsewhere
+ T(1,:) = -tau2/sqrt(ONE+tau_ratio**2) - bc%T0(1,:)
+ endwhere
- f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
+ where(T(2,:) > -VERYSMALLVAL)
+ T(2,:) = tau2*abs(tau_ratio)/sqrt(ONE+tau_ratio**2) - bc%T0(2,:)
+ elsewhere
+ T(2,:) = -tau2*abs(tau_ratio)/sqrt(ONE+tau_ratio**2) - bc%T0(2,:)
+ endwhere
- 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
+ if (any(ierr) == 1 ) then
+ print*, 'NR search 2nd loop failed!'
+ print*, 'iteration = ', it
+ call exit_MPI(myrank,'NR search 2nd loop failed!')
+ endif
-!=====================================================================
-subroutine rs_update_state(V,dt,f)
+ endif
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
- type(swf_type), intent(inout) :: f
- double precision, intent(in) :: dt
+ ! Save total tractions
+ bc%T = T
-vDtL = sqrt(V(1,:)**2 + V(2,:)**2) * dt * L
+ ! Subtract initial stress
+ T = T - bc%T0
-where(vDtL > 1.e-20_CUSTOM_REAL)
- f%theta = f%theta * exp(vDtL) + L/V * (1 - exp(vDtL))
-elsewhere
- f%theta = f%theta * exp(vDtL) + dt - 0.5*V/L*dt**2
-endwhere
+ ! 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)
- if (f%healing) then
- where(sqrt(V(1,:)**2 + V(2,:)**2) < V_HEALING)
- f%theta = 0e0_CUSTOM_REAL
- endwhere
- endif
-end subroutine rs_update_state
+ ! Update slip and slip rate, in fault frame
+ bc%D = dD
+ bc%V = dV + half_dt*dA
-!=====================================================================
-! Friction coefficient
-function swf_mu_TPV16(f,time,nglob) result(mu)
+ ! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- real(kind=CUSTOM_REAL) :: time
+ ! 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,:)
- integer :: i,nglob
+ 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,:)
- !-- linear slip weakening:
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- do i=1,nglob
- if ( (f%theta(i) >= f%Dc(i)) .or. (time >= f%T(i)) ) then
- mu(i) = f%mud(i)
+ !-- intermediate storage of outputs --
+ Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
+ Vnorm_old, Vnorm, it*bc%dt,bc%dt)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,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
- enddo
-
-end function swf_mu_TPV16
-!=====================================================================
-! Friction coefficient
-function swf_mu(f) result(mu)
+ end function rotate
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- !-- linear slip weakening:
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- mu = max( mu, f%mud)
-
-end function swf_mu
+ !=====================================================================
+ subroutine swf_update_state(dold,dnew,vold,f)
-!=====================================================================
-! Friction coefficient
-function rs_mu(f,V) result(mu)
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
+ type(swf_type), intent(inout) :: f
- type(rs_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- real(kind=CUSTOM_REAL) :: V(size(f%theta))
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
+ real(kind=CUSTOM_REAL) :: vnorm
+ integer :: k,npoin
- !-- rate and state friction (ageing law):
- mu = f%f0 + (f%a * log(sqrt(V(1,:)**2 + V(2,:)**2))/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
-
-end function rs_mu
+ f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
-!===============================================================
-! OUTPUTS
-subroutine init_dataT(DataT,coord,nglob,NT,iflt)
- ! NT = total number of time steps
+ 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
- integer, intent(in) :: nglob,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
- type (dataT_type), intent(out) :: DataT
+ !=====================================================================
+ subroutine rs_update_state(V,dt,f)
- real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
- integer :: i, iglob , IIN, ier, jflt, np, k
- character(len=70) :: tmpname
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ type(swf_type), intent(inout) :: f
+ double precision, intent(in) :: dt
- ! 1. read fault output coordinates from user file,
- ! 2. define iglob: the fault global index of the node nearest to user
- ! requested coordinate
+ vDtL = V * dt * f%L
- 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
+ where(vDtL > 1.e-20_CUSTOM_REAL)
+ f%theta = f%theta * exp(vDtL) + f%L/V * (1 - exp(vDtL))
+ elsewhere
+ f%theta = f%theta * exp(vDtL) + dt - 0.5*V/f%L*dt**2
+ endwhere
- allocate(DataT%iglob(DataT%npoin))
- allocate(DataT%name(DataT%npoin))
+ end subroutine rs_update_state
- 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)
+ !=====================================================================
+ ! Friction coefficient
+ function swf_mu_TPV16(f,time,nglob) result(mu)
- 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
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+ real(kind=CUSTOM_REAL) :: time
+
+ integer :: i,nglob
+
+ !-- linear slip weakening:
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ do i=1,nglob
+ if ( (f%theta(i) >= f%Dc(i)) .or. (time >= f%T(i)) ) then
+ mu(i) = f%mud(i)
+ 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))
- 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
- close(IIN)
+ end function swf_mu_TPV16
-end subroutine init_dataT
+ !=====================================================================
+ ! Friction coefficient
+ function swf_mu(f) result(mu)
-!---------------------------------------------------------------
-subroutine store_dataT(dataT,d,v,t,itime)
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- type(dataT_type), intent(inout) :: dataT
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
- integer, intent(in) :: itime
-
- integer :: i,k
+ !-- linear slip weakening:
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ mu = max( mu, f%mud)
- 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)
- enddo
+ end function swf_mu
-end subroutine store_dataT
+ !=====================================================================
+ ! Friction coefficient
+ function rs_mu(f,V) result(mu)
-!-----------------------------------------------------------------
-subroutine write_dataT_all(nt)
+ type(rs_type), intent(in) :: f
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: mu
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
- integer, intent(in) :: nt
-
- integer :: i
+ !-- rate and state friction (ageing law):
+ mu = f%f0 + (f%a * log(sqrt(V(1,:)**2 + V(2,:)**2))/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
- if (.not.allocated(faults)) return
- do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
- enddo
+ end function rs_mu
-end subroutine write_dataT_all
+ !===============================================================
+ ! OUTPUTS
+ subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+ ! NT = total number of time steps
-!------------------------------------------------------------------------
-subroutine SCEC_write_dataT(dataT,DT,NT)
+ integer, intent(in) :: nglob,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
+ type (dataT_type), intent(out) :: DataT
- type(dataT_type), intent(in) :: dataT
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
- integer :: i,k,IOUT
- character :: NTchar*5
- integer :: today(3), now(3)
+ ! 1. read fault output coordinates from user file,
+ ! 2. define iglob: the fault global index of the node nearest to user
+ ! requested coordinate
- 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
+ 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)
- write(NTchar,1) NT
- NTchar = adjustl(NTchar)
+ if (DataT%npoin == 0) return
-1 format(I5)
- do i=1,dataT%npoin
+ allocate(DataT%iglob(DataT%npoin))
+ allocate(DataT%name(DataT%npoin))
- open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "# problem=TPV16"
- 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=75 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)"
- write(IOUT,*) "#"
- write(IOUT,*) "# The line below lists the names of the data fields:"
- 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
- close(IOUT)
+ 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
- 1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+ ! 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))
+ 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
+ close(IIN)
- enddo
+ end subroutine init_dataT
-end subroutine SCEC_write_dataT
+ !---------------------------------------------------------------
+ subroutine store_dataT(dataT,d,v,t,itime)
-!-------------------------------------------------------------------------------------------------
-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)
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ integer, intent(in) :: itime
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+ 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)
+ enddo
- write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
+ end subroutine store_dataT
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
-
- open(IOUT,file=trim(filename),status='replace')
- write(IOUT,*) "# problem=TPV16"
- 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=75 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
+ !-----------------------------------------------------------------
+ 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=TPV16"
+ 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=75 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)"
+ write(IOUT,*) "#"
+ write(IOUT,*) "# The line below lists the names of the data fields:"
+ 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
+ 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=TPV16"
+ 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=75 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 )
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
-end subroutine SCEC_Write_RuptureTime
+ end subroutine SCEC_Write_RuptureTime
-!-------------------------------------------------------------------------------------------------
-subroutine init_dataXZ(DataXZ,bc,nglob)
+ !-------------------------------------------------------------------------------------------------
+ subroutine init_dataXZ(DataXZ,bc,nglob)
- type(dataXZ_type), intent(inout) :: DataXZ
- type(bc_dynflt_type) :: bc
- integer, intent(in) :: nglob
+ type(dataXZ_type), intent(inout) :: DataXZ
+ type(bc_dynflt_type) :: bc
+ integer, intent(in) :: nglob
- allocate(DataXZ%stg(nglob))
- DataXZ%sta => bc%swf%theta
- 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))
+ allocate(DataXZ%stg(nglob))
+ DataXZ%sta => bc%swf%theta
+ 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
+ !Percy , setting up initial rupture time null for all faults.
+ DataXZ%tRUP = 0e0_CUSTOM_REAL
+ DataXZ%tPZ = 0e0_CUSTOM_REAL
-end subroutine init_dataXZ
+ end subroutine init_dataXZ
-!---------------------------------------------------------------
-subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
+ !---------------------------------------------------------------
+ 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
+ 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
+ integer :: i
-! "stg" : strength .
-
- dataXZ%stg = stg
+ ! "stg" : strength .
- 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).
+ dataXZ%stg = stg
- ! note: the other arrays in dataXZ are pointers to arrays in bc
- ! they do not need to be updated here
+ 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 .
-end subroutine store_dataXZ
+ 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
-!---------------------------------------------------------------
-subroutine write_dataXZ(dataXZ,itime,iflt)
+ ! To do : add stress criteria (firs time strength is reached).
- type(dataXZ_type), intent(in) :: dataXZ
- integer, intent(in) :: itime,iflt
-
- character(len=70) :: filename
+ ! note: the other arrays in dataXZ are pointers to arrays in bc
+ ! they do not need to be updated here
- 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.
+ end subroutine store_dataXZ
-! 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')
+ !---------------------------------------------------------------
+ subroutine write_dataXZ(dataXZ,itime,iflt)
- 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)
+ type(dataXZ_type), intent(in) :: dataXZ
+ integer, intent(in) :: itime,iflt
-end subroutine write_dataXZ
+ 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
+
+ !---------------------------------------------------------------
+
+ ! NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror)
+ ! Newton-Raphson search used to find slip rates
+ !
+ ! When Newton-Raphson search fails (happens near the free surface points),
+ ! use bisection method which will not fail to find the root if the specified
+ ! interval contains a root.
+ ! NOTE: When an interval contains more than one root, the bisection method
+ ! can find one of them.
+ !
+ !
+ ! INPUT VARIABLES fo,Vo,cca,ccb,Seff,tau,FaultVFree,FaultZ,ierror
+ !
+ ! OUTPUT VARIABLES Vfout,tauout
+ !
+ ! Yoshihiro Kaneko: ykaneko at gps.caltech.edu
+ !
+ subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror)
+
+ implicit none
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL) :: cA,eps,delta,fa,help,help1,help2,Vfprime
+ real(kind=CUSTOM_REAL) :: Vfout,tauout,fo,Vo,cca,ccb,Seff,tau,psi,FaultZ,TauStick
+ real(kind=CUSTOM_REAL) :: f_bisec,f_bisec_Vflow,Vfbisec,Vfhigh,Vflow
+ integer :: j,k,ierror
+
+ cA = cca*Seff
+ eps = 0.001*cA
+ k=0
+ delta = HUGEVAL
+
+ DO WHILE ( abs(delta) > eps )
+ fa = tau/(Seff*cca)
+ help = -(fo+ccb*psi)/cca
+ help1 = exp(help+fa)
+ help2 = exp(help-fa)
+ Vfout = Vo*(help1 - help2)
+ Vfprime = (Vo/cA)*(help1 + help2)
+ delta = (TauStick - HALF*FaultZ*Vfout - tau)/(ONE + HALF*FaultZ*Vfprime)
+ tau = tau + delta
+ k = k + 1
+ if ( abs(delta) > 1.e10 .OR. k > 100 ) then
+ ! try bisection if NR search fails
+ Vflow = 0._CUSTOM_REAL ! lower bound of Vf for bisection
+ Vfhigh = 5._CUSTOM_REAL ! upper bound of Vf for bisection
+ do j = 1,200
+ Vfbisec = (Vflow + Vfhigh)/TWO
+ f_bisec = TauStick - HALF*FaultZ*Vfbisec - cA*asinh(Vfbisec/(TWO*Vo)*exp(-help))
+ if ( abs(f_bisec) < eps ) then
+ Vfout = Vfbisec
+ tauout = TauStick - HALF*FaultZ*Vfout
+ return
+ endif
+ f_bisec_Vflow = TauStick - HALF*FaultZ*Vflow - cA*asinh(Vflow/(TWO*Vo)*exp(-help))
+ if ( f_bisec*f_bisec_Vflow < 0._CUSTOM_REAL ) then
+ Vfhigh = Vfbisec
+ else
+ Vflow = Vfbisec
+ endif
+ enddo
+
+ ierror = 1
+ !compute updated Vf before exiting NR search
+ fa = tau/(Seff*cca)
+ help = -(fo+ccb*psi)/cca
+ help1 = exp(help+fa)
+ help2 = exp(help-fa)
+ Vfout = Vo*(help1 - help2)
+ tauout = tau
+ return
+ endif
+
+ END DO
+
+ !compute the updated Vf before exiting NR search
+ fa = tau/(Seff*cca)
+ help = -(fo+ccb*psi)/cca
+ help1 = exp(help+fa)
+ help2 = exp(help-fa)
+ Vfout = Vo*(help1 - help2)
+ tauout = tau
+
+ end subroutine NRsearch
+
+ !---------------------------------------------------------------
+
+
end module fault_solver
More information about the CIG-COMMITS
mailing list