[cig-commits] r19867 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
ampuero at geodynamics.org
ampuero at geodynamics.org
Mon Mar 26 00:46:40 PDT 2012
Author: ampuero
Date: 2012-03-26 00:46:39 -0700 (Mon, 26 Mar 2012)
New Revision: 19867
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
many fixes in rsf solver
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-26 06:05:51 UTC (rev 19866)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-26 07:46:39 UTC (rev 19867)
@@ -30,93 +30,96 @@
module fault_solver
- implicit none
+ implicit none
- include 'constants.h'
+ include 'constants.h'
- private
+ 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
- character(len=70), dimension(:), pointer :: name
- end type dataT_type
+ ! 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
+ 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
+ ! 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 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 rs_type
- 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 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 :: AspTx=>null(), Fload=>null()
- end type asperity_type
+ type asperity_type
+ private
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: AspTx=>null(), 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(rs_type), pointer :: rs => 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
+ 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 :: rs => 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(:)
+ 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 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
+ !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
+ !Number of time steps defined by the user : NTOUT
+ integer, save :: NTOUT,NSNAP
- logical, save :: SIMULATION_TYPE_DYN = .false.
+ logical, save :: SIMULATION_TYPE_DYN = .false.
- logical, save :: TPV16 = .false.
+ logical, save :: TPV16 = .false.
- logical, save :: Rate_AND_State = .true.
+ logical, save :: Rate_AND_State = .true.
- logical, save :: SlipLaw = .false.
+ real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
- real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
+ public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
+ SIMULATION_TYPE_DYN
- public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
- SIMULATION_TYPE_DYN
-
contains
!=====================================================================
@@ -128,86 +131,86 @@
!
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
+ 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_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)
+ 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)
- 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
+ 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
+ ! WARNING TO DO: should be an MPI abort
end subroutine BC_DYNFLT_init
@@ -216,217 +219,217 @@
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)
- enddo
- enddo
+ 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
- !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
+ 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
- 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(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)
+ 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)
- 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)
+ 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)
- 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) )
+ 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)
+ ! 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
- do iLoc=1,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)
+ 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%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)
+ 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)
- ! WARNING: if V_HEALING is negative we turn off healing
- bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+ ! 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)
+ bc%swf%theta = 0e0_CUSTOM_REAL
+ allocate(bc%MU(bc%nglob))
+ bc%MU = swf_mu_TPV16(bc%swf,0.0,bc%nglob)
- 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)
- 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,:)
+ 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,:)
end subroutine init_one_fault_TPV16
@@ -435,252 +438,251 @@
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
+ 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), 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
- 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
- NAMELIST / ASP / AspTx,Fload,nAsp,nFload
+ NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
+ NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
+ NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init
+ NAMELIST / ASP / AspTx,Fload,nAsp,nFload
- 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)
- 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%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)
+ 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
+ else !Rate AND State
+ ! Set friction parameters and initialize friction variables
+ 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) )
+ ! WARNING: if V_HEALING is negative we turn off healing
- V0 =1.e-6_CUSTOM_REAL
- f0 =0.6_CUSTOM_REAL
- a =0.0220_CUSTOM_REAL !0.0080_CUSTOM_REAL
- b =0.0050_CUSTOM_REAL !0.0120_CUSTOM_REAL
- 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.0220_CUSTOM_REAL !0.0080_CUSTOM_REAL
+ b =0.0050_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
+ 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=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
+ 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)
- bc%rs%theta = bc%rs%theta_init
- allocate(bc%MU(bc%nglob)) !Not necessary?
- bc%MU = rs_mu(bc%rs,bc%rs%V_init)
- bc%T0(3,:) = bc%T0(3,:) + 1.e6_CUSTOM_REAL
- bc%T0(1,:) = bc%MU * bc%T0(3,:)
- bc%T0(2,:) = 0._CUSTOM_REAL
+ allocate(bc%MU(bc%nglob)) !Not necessary?
+ bc%MU = rsf_mu(bc%rsf,bc%rsf%V_init)
+! bc%T0(3,:) = bc%T0(3,:) + 1.e6_CUSTOM_REAL
+ bc%T0(1,:) = bc%MU * bc%T0(3,:)
+ bc%T0(2,:) = 0._CUSTOM_REAL
- allocate( bc%asp )
- allocate( bc%asp%AspTx(bc%nglob) )
- allocate( bc%asp%Fload(bc%nglob) )
+ allocate( bc%asp )
+ allocate( bc%asp%AspTx(bc%nglob) )
+ allocate( bc%asp%Fload(bc%nglob) )
- AspTx =0.e0_CUSTOM_REAL
- Fload =0.e0_CUSTOM_REAL
- nAsp =0
- nFload =0
+ AspTx =0.e0_CUSTOM_REAL
+ Fload =0.e0_CUSTOM_REAL
+ nAsp =0
+ nFload =0
- read(IIN_PAR, nml=ASP)
- bc%asp%AspTx =AspTx
- bc%asp%Fload =Fload
- call init_2d_distribution(bc%asp%AspTx,bc%coord,IIN_PAR,nAsp)
- call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
+ read(IIN_PAR, nml=ASP)
+ bc%asp%AspTx =AspTx
+ bc%asp%Fload =Fload
+ call init_2d_distribution(bc%asp%AspTx,bc%coord,IIN_PAR,nAsp)
+ call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
- endif
+ 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
!---------------------------------------------------------------------
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
+ 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
+ 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
+ ! 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
+ 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
+ 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
@@ -688,72 +690,72 @@
! 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), 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
- integer :: i
+ 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
- NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
- if (n==0) return
+ 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
+ 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 ('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
+ 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
+ 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
+ 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
+ if (x>=0e0_CUSTOM_REAL) then
+ heaviside = 1e0_CUSTOM_REAL
+ else
+ heaviside = 0e0_CUSTOM_REAL
+ endif
end function heaviside
@@ -762,388 +764,384 @@
!
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
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
- integer :: iflt
+ 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)
- else
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d_TPV16(faults(iflt),F,Vel,Dis,iflt)
- endif
- enddo
+ 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
!---------------------------------------------------------------------
subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
- use specfem_par, only:it,NSTEP
+ 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), 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) :: 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
+ 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,:))
+ 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,:))
- ! 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
+ ! 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)
+ ! 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,:) )
+ ! 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
+ !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
+ ! 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)
+ ! 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)
+ ! 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_TPV16(bc%swf,it*bc%dt,bc%nglob)
+ ! 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) )
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
- ! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
+ ! 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
+ ! 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
+ tnew = min(tnorm,strength)
+ T(1,:) = tnew * t1
+ T(2,:) = tnew * t2
- ! Save total tractions
- bc%T = T
+ ! Save total tractions
+ bc%T = T
- ! Subtract initial stress
- T = T - bc%T0
+ ! 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 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
+ ! 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)
+ ! 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,:)
+ ! 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,:)
+ 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,:))
- 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)
+ !-- 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)
+ !-- 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_TPV16
!---------------------------------------------------------------------
subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
- use specfem_par, only:it,NSTEP
+ 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), 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) :: 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
+ 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
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,ta,tau_ratio
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: Vf,Vf1,tau1,Vf2,tau2,Vfavg
- integer, dimension(bc%nglob) :: ierr
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: TxExt
- real(kind=CUSTOM_REAL) :: TLoad,GLoad
- real(kind=CUSTOM_REAL) :: theta1,thetan
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,ta
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: Vf,Vf1,tau1,Vf2,tau2,Vfavg
+ integer, dimension(bc%nglob) :: ierr
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: TxExt
+ real(kind=CUSTOM_REAL) :: TLoad,DTau0,GLoad
+ real(kind=CUSTOM_REAL) :: thetan,time
- 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,:))
+ 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,:))
- ! 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
+ ! 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)
+ ! 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,:) )
+ ! 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
+ !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
+ ! 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)
+ ! 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
- call swf_update_state(bc%D,dD,bc%V,bc%swf)
+ 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)
+ ! 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) )
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ ! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
- ! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
+ ! 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
- ! 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
- tnew = min(tnorm,strength)
- T(1,:) = tnew * t1
- T(2,:) = tnew * t2
+ else ! Update rate and state friction:
- 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)
- thetan = bc%rs%theta
- call rs_update_state(Vf,bc%dt,bc%rs)
+ ! smooth loading within nucleation patch
+ !WARNING : ad hoc for SCEC becnhmark TPV10x
+ TxExt = 0._CUSTOM_REAL
+ TLoad = 1.0_CUSTOM_REAL
+ DTau0 = 25e6_CUSTOM_REAL
+ time = it*bc%dt
+ if (time==0._CUSTOM_REAL) then
+ GLoad = 0._CUSTOM_REAL
+ elseif (time <= TLoad) then
+ 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
- ! smooth loading within nucleation patch
- TxExt = 0._CUSTOM_REAL
- TLoad = 0.1_CUSTOM_REAL
- if(it*bc%dt <= TLoad) then
- GLoad = exp(it*bc%dt - TLoad)**2 / (it*bc%dt*(it*bc%dt-2.0*TLoad))
- TxExt = (bc%asp%AspTx - bc%T0(1,:))*bc%asp%Fload*GLoad*bc%T0(1,:)
- endif
+ tStick = sqrt(T(1,:)**2 + T(2,:)**2)
- tStick = sqrt( (T(1,:) - bc%T0(1,:) + TxExt)**2 + T(2,:)**2)
- ! tStick = sqrt(T(1,:)**2 + T(2,:)**2)
+ Vf = sqrt(bc%V(1,:)**2 + bc%V(2,:)**2)
+ thetan = bc%rsf%theta
+ call rsf_update_state(Vf,bc%dt,bc%rsf)
- tau_ratio = T(2,:)/T(1,:)
- ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
- call NRsearch(Vf1,tau1,bc%rs%f0,bc%rs%V0,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr,bc%nglob)
+ ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
+ call NRsearch(Vf1,tau1,bc%rsf%f0,bc%rsf%V0,bc%rsf%a,bc%rsf%b,-T(3,:),log(bc%rsf%theta),bc%Z,ta,tStick,ierr,bc%nglob)
- ! 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
+ ! 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
- if (any(ierr == 1) ) then
- print*, 'NR search 1st loop failed!'
- print*, 'iteration = ', it
- !call exit_MPI(myrank,'NR search 2nd loop failed!')
- stop
- endif
- theta1 = bc%rs%theta
- bc%rs%theta = thetan
- ! Updating state variable: 2nd loop
- Vfavg = 0.5e0_CUSTOM_REAL**abs(Vf + Vf1)
- call rs_update_state(Vfavg,bc%dt,bc%rs)
+ if (any(ierr == 1) ) then
+ print*, 'NR search 1st loop failed!'
+ print*, 'iteration = ', it
+ !call exit_MPI(myrank,'NR search 2nd loop failed!')
+ stop
+ endif
+ ! Updating state variable: 2nd loop
+ Vfavg = 0.5e0_CUSTOM_REAL**abs(Vf + Vf1)
+ bc%rsf%theta = thetan
+ call rsf_update_state(Vfavg,bc%dt,bc%rsf)
- ! NR search 2nd loop
- call NRsearch(Vf2,tau2,bc%rs%f0,bc%rs%V0,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr,bc%nglob)
+ ! NR search 2nd loop
+ call NRsearch(Vf2,tau2,bc%rsf%f0,bc%rsf%V0,bc%rsf%a,bc%rsf%b,-T(3,:),log(bc%rsf%theta),bc%Z,ta,tStick,ierr,bc%nglob)
- where (Vf2 < 0._CUSTOM_REAL .or. tau2 < 0._CUSTOM_REAL .or. ierr == 1)
- tau2 = tStick
- Vf2 = 0._CUSTOM_REAL
- endwhere
+ where (Vf2 < 0._CUSTOM_REAL .or. tau2 < 0._CUSTOM_REAL .or. ierr == 1)
+ tau2 = tStick
+ Vf2 = 0._CUSTOM_REAL
+ endwhere
- ! 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
+ if (any(ierr == 1) ) then
+ print*, 'NR search 2nd loop failed!'
+ print*, 'iteration = ', it
+ !call exit_MPI(myrank,'NR search 2nd loop failed!')
+ stop
+ endif
- 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
+ tStick = max(tStick,1e0_CUSTOM_REAL)
+ t1 = T(1,:)/tnorm
+ t2 = T(2,:)/tnorm
+ T(1,:) = tau2 * t1
+ T(2,:) = tau2 * t2
- if (any(ierr == 1) ) then
- print*, 'NR search 2nd loop failed!'
- print*, 'iteration = ', it
- !call exit_MPI(myrank,'NR search 2nd loop failed!')
- stop
- endif
+ endif
- endif
+ ! Save total tractions
+ bc%T = T
- ! Save total tractions
- bc%T = T
+ ! Subtract initial stress
+ T = T - bc%T0
- ! 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 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
- ! 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)
- ! 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,:)
- ! 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,:)
- 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,:))
+ 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)
- !-- 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)
- !-- 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)
+ 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)
+ ! 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(:,:)
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
- real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
+ 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)
+ ! 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)
+ 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
+ ! 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
+ ! 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
+ ! 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
+ endif
end function rotate
@@ -1151,63 +1149,66 @@
!=====================================================================
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), dimension(:,:), intent(in) :: vold,dold,dnew
+ type(swf_type), intent(inout) :: f
- real(kind=CUSTOM_REAL) :: vnorm
- integer :: k,npoin
+ real(kind=CUSTOM_REAL) :: vnorm
+ integer :: k,npoin
- f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
+ 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
+ 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
!=====================================================================
-subroutine rs_update_state(V,dt,f)
+subroutine rsf_update_state(V,dt,f)
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
- type(rs_type), intent(inout) :: f
- real(kind=CUSTOM_REAL), intent(in) :: dt
+ 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))
+ real(kind=CUSTOM_REAL) :: vDtL(size(V))
- vDtL = V * dt * f%L
+ vDtL = V * dt * f%L
- if(.not. SlipLaw) 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 - HALF*V/f%L*dt**2
- endwhere
- else
- f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
- endif
+ ! 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 - HALF*V/f%L*dt**2
+ endwhere
-end subroutine rs_update_state
+ ! slip law
+ else
+ f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
+ endif
+end subroutine rsf_update_state
+
!=====================================================================
! Friction coefficient
function swf_mu_TPV16(f,time,nglob) result(mu)
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- real(kind=CUSTOM_REAL) :: time
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+ real(kind=CUSTOM_REAL) :: time
- integer :: i,nglob
+ 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
+ !-- 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
end function swf_mu_TPV16
@@ -1215,236 +1216,327 @@
! Friction coefficient
function swf_mu(f) result(mu)
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+ 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)
+ !-- linear slip weakening:
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ mu = max( mu, f%mud)
end function swf_mu
!=====================================================================
-! Friction coefficient
-function rs_mu(f,V) result(mu)
+! Rate and state friction coefficient
+function rsf_mu(f,V) result(mu)
- type(rs_type), intent(in) :: f
- !real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
- !real(kind=CUSTOM_REAL) :: mu(size(V,2))
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
- real(kind=CUSTOM_REAL) :: mu(size(V))
+ type(rsf_type), intent(in) :: f
+ !real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
+ !real(kind=CUSTOM_REAL) :: mu(size(V,2))
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ real(kind=CUSTOM_REAL) :: mu(size(V))
- !-- 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))
- mu = f%f0 + (f%a * log(V/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
+ !mu = f%f0 + (f%a * log(sqrt(V(1,:)**2 + V(2,:)**2)/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
+ mu = f%f0 + (f%a * log(V/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
-end function rs_mu
+end function rsf_mu
!===============================================================
+! Newton-Raphson search for rate-and-state solver
+!
+! 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
+! Surendra : Modified to loop over fault nodes
+subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror,nglob)
+
+ !implicit none
+ !include "constants.h"
+
+ integer :: iglob,nglob
+
+ real(kind=CUSTOM_REAL), dimension(nglob) :: cA,eps,delta,fa,help,help1,help2,Vfprime
+ real(kind=CUSTOM_REAL), dimension(nglob) :: Vfout,tauout,fo,Vo,cca,ccb,Seff,tau,psi,FaultZ,TauStick
+ real(kind=CUSTOM_REAL), dimension(nglob) :: f_bisec,f_bisec_Vflow,Vfbisec,Vfhigh,Vflow
+ integer :: j,k
+ integer, dimension(nglob) :: ierror
+
+
+ do iglob=1,nglob
+ cA(iglob) = cca(iglob)*Seff(iglob)
+ eps(iglob) = 0.001*cA(iglob)
+ k=0
+ delta(iglob) = HUGEVAL
+
+ DO WHILE ( abs(delta(iglob)) > eps(iglob) )
+ fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
+ help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
+ help1(iglob) = exp(help(iglob)+fa(iglob))
+ help2(iglob) = exp(help(iglob)-fa(iglob))
+ Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
+ Vfprime(iglob) = (Vo(iglob)/cA(iglob))*(help1(iglob) + help2(iglob))
+ delta(iglob) = (TauStick(iglob) - HALF*FaultZ(iglob)*Vfout(iglob) - tau(iglob))/(ONE + HALF*FaultZ(iglob)*Vfprime(iglob))
+ tau(iglob) = tau(iglob) + delta(iglob)
+ k = k + 1
+ if ( abs(delta(iglob)) > 1.e10 .OR. k > 100 ) then
+ ! try bisection if NR search fails
+ Vflow(iglob) = 0._CUSTOM_REAL ! lower bound of Vf for bisection
+ Vfhigh(iglob) = 5._CUSTOM_REAL ! upper bound of Vf for bisection
+ do j = 1,200
+ Vfbisec(iglob) = (Vflow(iglob) + Vfhigh(iglob))/TWO
+ f_bisec(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vfbisec(iglob) - cA(iglob)*asinh(Vfbisec(iglob)/(TWO*Vo(iglob))*exp(-help(iglob)))
+ if ( abs(f_bisec(iglob)) < eps(iglob) ) then
+ Vfout(iglob) = Vfbisec(iglob)
+ tauout(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vfout(iglob)
+ return
+ endif
+ f_bisec_Vflow(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vflow(iglob) - cA(iglob)*asinh(Vflow(iglob)/(TWO*Vo(iglob))*exp(-help(iglob)))
+ if ( f_bisec(iglob)*f_bisec_Vflow(iglob) < 0._CUSTOM_REAL ) then
+ Vfhigh(iglob) = Vfbisec(iglob)
+ else
+ Vflow(iglob) = Vfbisec(iglob)
+ endif
+ enddo
+
+ ierror(iglob) = 1
+ !compute updated Vf before exiting NR search
+ fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
+ help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
+ help1(iglob) = exp(help(iglob)+fa(iglob))
+ help2(iglob) = exp(help(iglob)-fa(iglob))
+ Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
+ tauout(iglob) = tau(iglob)
+ return
+ endif
+
+ END DO
+
+ !compute the updated Vf before exiting NR search
+ fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
+ help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
+ help1(iglob) = exp(help(iglob)+fa(iglob))
+ help2(iglob) = exp(help(iglob)-fa(iglob))
+ Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
+ tauout(iglob) = tau(iglob)
+ enddo
+end subroutine NRsearch
+
+!===============================================================
! OUTPUTS
subroutine init_dataT(DataT,coord,nglob,NT,iflt)
- ! NT = total number of time steps
+ ! 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
+ 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
+ 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
+ ! 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)
+ 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
+ if (DataT%npoin == 0) return
- allocate(DataT%iglob(DataT%npoin))
- allocate(DataT%name(DataT%npoin))
+ 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)
+ 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
+ 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))
- 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
+ ! 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)
+ close(IIN)
end subroutine init_dataT
!---------------------------------------------------------------
subroutine store_dataT(dataT,d,v,t,itime)
- type(dataT_type), intent(inout) :: dataT
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
- integer, intent(in) :: itime
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ integer, intent(in) :: itime
- integer :: i,k
+ 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
+ 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 subroutine store_dataT
!-----------------------------------------------------------------
subroutine write_dataT_all(nt)
- integer, intent(in) :: nt
+ integer, intent(in) :: nt
- integer :: i
+ integer :: i
- if (.not.allocated(faults)) return
- do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
- enddo
+ 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
+ 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)
+ 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
+ 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
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
- write(NTchar,1) NT
- NTchar = adjustl(NTchar)
+ write(NTchar,1) NT
+ NTchar = adjustl(NTchar)
1 format(I5)
- do i=1,dataT%npoin
+ 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)
+ 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
+ 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
+ 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)
+ 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
+ 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
+ write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+ 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
+ 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)
+ close(IOUT)
1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
@@ -1453,195 +1545,102 @@
!-------------------------------------------------------------------------------------------------
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
!---------------------------------------------------------------
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 .
+ ! "stg" : strength .
- dataXZ%stg = stg
+ 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 .
+ 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
+ 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).
+ ! 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
+ ! 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
+ type(dataXZ_type), intent(in) :: dataXZ
+ integer, intent(in) :: itime,iflt
- character(len=70) :: filename
+ 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(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
+ ! 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')
+ 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)
+ 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
-! Surendra : Modified to loop over fault nodes
-subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror,nglob)
-
- !implicit none
- !include "constants.h"
-
- integer :: iglob,nglob
-
- real(kind=CUSTOM_REAL), dimension(nglob) :: cA,eps,delta,fa,help,help1,help2,Vfprime
- real(kind=CUSTOM_REAL), dimension(nglob) :: Vfout,tauout,fo,Vo,cca,ccb,Seff,tau,psi,FaultZ,TauStick
- real(kind=CUSTOM_REAL), dimension(nglob) :: f_bisec,f_bisec_Vflow,Vfbisec,Vfhigh,Vflow
- integer :: j,k
- integer, dimension(nglob) :: ierror
-
-
- do iglob=1,nglob
- cA(iglob) = cca(iglob)*Seff(iglob)
- eps(iglob) = 0.001*cA(iglob)
- k=0
- delta(iglob) = HUGEVAL
-
- DO WHILE ( abs(delta(iglob)) > eps(iglob) )
- fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
- help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
- help1(iglob) = exp(help(iglob)+fa(iglob))
- help2(iglob) = exp(help(iglob)-fa(iglob))
- Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
- Vfprime(iglob) = (Vo(iglob)/cA(iglob))*(help1(iglob) + help2(iglob))
- delta(iglob) = (TauStick(iglob) - HALF*FaultZ(iglob)*Vfout(iglob) - tau(iglob))/(ONE + HALF*FaultZ(iglob)*Vfprime(iglob))
- tau(iglob) = tau(iglob) + delta(iglob)
- k = k + 1
- if ( abs(delta(iglob)) > 1.e10 .OR. k > 100 ) then
- ! try bisection if NR search fails
- Vflow(iglob) = 0._CUSTOM_REAL ! lower bound of Vf for bisection
- Vfhigh(iglob) = 5._CUSTOM_REAL ! upper bound of Vf for bisection
- do j = 1,200
- Vfbisec(iglob) = (Vflow(iglob) + Vfhigh(iglob))/TWO
- f_bisec(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vfbisec(iglob) - cA(iglob)*asinh(Vfbisec(iglob)/(TWO*Vo(iglob))*exp(-help(iglob)))
- if ( abs(f_bisec(iglob)) < eps(iglob) ) then
- Vfout(iglob) = Vfbisec(iglob)
- tauout(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vfout(iglob)
- return
- endif
- f_bisec_Vflow(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vflow(iglob) - cA(iglob)*asinh(Vflow(iglob)/(TWO*Vo(iglob))*exp(-help(iglob)))
- if ( f_bisec(iglob)*f_bisec_Vflow(iglob) < 0._CUSTOM_REAL ) then
- Vfhigh(iglob) = Vfbisec(iglob)
- else
- Vflow(iglob) = Vfbisec(iglob)
- endif
- enddo
-
- ierror(iglob) = 1
- !compute updated Vf before exiting NR search
- fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
- help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
- help1(iglob) = exp(help(iglob)+fa(iglob))
- help2(iglob) = exp(help(iglob)-fa(iglob))
- Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
- tauout(iglob) = tau(iglob)
- return
- endif
-
- END DO
-
- !compute the updated Vf before exiting NR search
- fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
- help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
- help1(iglob) = exp(help(iglob)+fa(iglob))
- help2(iglob) = exp(help(iglob)-fa(iglob))
- Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
- tauout(iglob) = tau(iglob)
- enddo
-end subroutine NRsearch
-
-!---------------------------------------------------------------
-
-
end module fault_solver
More information about the CIG-COMMITS
mailing list