[cig-commits] r19864 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Sun Mar 25 08:54:40 PDT 2012
Author: surendra
Date: 2012-03-25 08:54:40 -0700 (Sun, 25 Mar 2012)
New Revision: 19864
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Added asperity (time-dependent nucleation) loading for RS friction. Corrected for indentation. Appropriate state variable is used during the call to rs_update_state
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-24 17:46:48 UTC (rev 19863)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-25 15:54:40 UTC (rev 19864)
@@ -30,1572 +30,1618 @@
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 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 bc_dynflt_type
- private
- integer :: nspec,nglob
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0,T,V,D
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: coord
- real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R
- real(kind=CUSTOM_REAL), dimension(:), pointer :: MU,B,invM1,invM2,Z
- real(kind=CUSTOM_REAL) :: dt
- integer, dimension(:), pointer :: ibulk1, ibulk2
- type(swf_type), pointer :: swf => null()
- type(rs_type), pointer :: rs => null()
- logical :: allow_opening = .false. ! default : do not allow opening
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ
- end type bc_dynflt_type
+ type asperity_type
+ private
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: AspTx=>null(), Fload=>null()
+ end type asperity_type
- type(bc_dynflt_type), allocatable, save :: faults(:)
+ 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
- !slip velocity threshold for healing
- !WARNING: not very robust
- real(kind=CUSTOM_REAL), save :: V_HEALING
+ type(bc_dynflt_type), allocatable, save :: faults(:)
- !slip velocity threshold for definition of rupture front
- real(kind=CUSTOM_REAL), save :: V_RUPT
+ !slip velocity threshold for healing
+ !WARNING: not very robust
+ real(kind=CUSTOM_REAL), save :: V_HEALING
- !Number of time steps defined by the user : NTOUT
- integer, save :: NTOUT,NSNAP
+ !slip velocity threshold for definition of rupture front
+ real(kind=CUSTOM_REAL), save :: V_RUPT
- logical, save :: SIMULATION_TYPE_DYN = .false.
+ !Number of time steps defined by the user : NTOUT
+ integer, save :: NTOUT,NSNAP
- logical, save :: TPV16 = .false.
+ logical, save :: SIMULATION_TYPE_DYN = .false.
- logical, save :: Rate_AND_State = .false.
+ logical, save :: TPV16 = .false.
- logical, save :: SlipLaw = .false.
+ logical, save :: Rate_AND_State = .true.
- real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
+ logical, save :: SlipLaw = .false.
- public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
- SIMULATION_TYPE_DYN
+ real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
+ public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
+ SIMULATION_TYPE_DYN
+
contains
+!=====================================================================
+! BC_DYNFLT_init initializes dynamic faults
+!
+! prname fault database is read from file prname_fault_db.bin
+! Minv inverse mass matrix
+! dt global time step
+!
+subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt)
- !=====================================================================
- ! BC_DYNFLT_init initializes dynamic faults
- !
- ! prname fault database is read from file prname_fault_db.bin
- ! Minv inverse mass matrix
- ! dt global time step
- !
- subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt)
+ character(len=256), intent(in) :: prname ! 'proc***'
+ real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
+ double precision, intent(in) :: DTglobal
+ integer, intent(in) :: nt
- character(len=256), intent(in) :: prname ! 'proc***'
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
- double precision, intent(in) :: DTglobal
- integer, intent(in) :: nt
+ real(kind=CUSTOM_REAL) :: dt
+ integer :: iflt,ier,dummy_idfault
+ integer :: nbfaults
+ integer :: size_Kelvin_Voigt
+ integer :: SIMULATION_TYPE
+ character(len=256) :: filename
+ integer, parameter :: IIN_PAR =151
+ integer, parameter :: IIN_BIN =170
- real(kind=CUSTOM_REAL) :: dt
- integer :: iflt,ier,dummy_idfault
- integer :: nbfaults
- integer :: size_Kelvin_Voigt
- integer :: SIMULATION_TYPE
- character(len=256) :: filename
- integer, parameter :: IIN_PAR =151
- integer, parameter :: IIN_BIN =170
+ NAMELIST / BEGIN_FAULT / dummy_idfault
- NAMELIST / BEGIN_FAULT / dummy_idfault
+ dummy_idfault = 0
- dummy_idfault = 0
+ open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults',status='old',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File Par_file_faults not found: assume no faults'
+ close(IIN_PAR)
+ return
+ endif
- open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults',status='old',iostat=ier)
- if( ier /= 0 ) then
- write(6,*) 'File Par_file_faults not found: assume no faults'
- close(IIN_PAR)
- return
- endif
+ dt = real(DTglobal)
+ filename = prname(1:len_trim(prname))//'fault_db.bin'
+ open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File ',trim(filename),' not found. Abort'
+ stop
+ endif
+ ! WARNING TO DO: should be an MPI abort
- 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
+end subroutine BC_DYNFLT_init
- !---------------------------------------------------------------------
+!---------------------------------------------------------------------
- subroutine init_one_fault_TPV16(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+subroutine init_one_fault_TPV16(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
- type(bc_dynflt_type), intent(inout) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
- integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: dt
+ type(bc_dynflt_type), intent(inout) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
+ integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: dt
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
- integer, dimension(:,:), allocatable :: ibool1
- real(kind=CUSTOM_REAL) :: norm
- real(kind=CUSTOM_REAL) :: S1,S2,S3
- integer :: n1,n2,n3
- real(kind=CUSTOM_REAL) :: mus,mud,dc
- integer :: nmus,nmud,ndc,ij,k,e
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
+ integer, dimension(:,:), allocatable :: ibool1
+ real(kind=CUSTOM_REAL) :: norm
+ real(kind=CUSTOM_REAL) :: S1,S2,S3
+ integer :: n1,n2,n3
+ real(kind=CUSTOM_REAL) :: mus,mud,dc
+ integer :: nmus,nmud,ndc,ij,k,e
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
- integer, parameter :: IIN_NUC =270
- integer :: ipar
- integer :: ier
+ integer, parameter :: IIN_NUC =270
+ integer :: ipar
+ integer :: ier
- integer :: relz_num,sub_relz_num
- integer :: num_cell_str,num_cell_dip
- real(kind=CUSTOM_REAL) :: siz_str,siz_dip
- integer :: hypo_cell_str,hypo_cell_dip
- real(kind=CUSTOM_REAL) :: hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ integer :: relz_num,sub_relz_num
+ integer :: num_cell_str,num_cell_dip
+ real(kind=CUSTOM_REAL) :: siz_str,siz_dip
+ integer :: hypo_cell_str,hypo_cell_dip
+ real(kind=CUSTOM_REAL) :: hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
- integer, dimension(:), allocatable :: inp_nx,inp_nz
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
- dyn_fc,swcd,cohes,tim_forcedRup
+ integer, dimension(:), allocatable :: inp_nx,inp_nz
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
+ dyn_fc,swcd,cohes,tim_forcedRup
- !integer, dimension(:), allocatable :: iLoc
+ !integer, dimension(:), allocatable :: iLoc
- integer :: iLoc
+ integer :: iLoc
- NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
+ NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
+ NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
- read(IIN_BIN) bc%nspec,bc%nglob
- if (bc%nspec==0) return
- allocate( bc%ibulk1(bc%nglob) )
- allocate( bc%ibulk2(bc%nglob) )
- allocate( ibool1(NGLLSQUARE,bc%nspec) )
- allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
- allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
+ read(IIN_BIN) bc%nspec,bc%nglob
+ if (bc%nspec==0) return
+ allocate( bc%ibulk1(bc%nglob) )
+ allocate( bc%ibulk2(bc%nglob) )
+ allocate( ibool1(NGLLSQUARE,bc%nspec) )
+ allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
+ allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
- allocate(bc%coord(3,(bc%nglob)))
- read(IIN_BIN) ibool1
- read(IIN_BIN) jacobian2Dw
- read(IIN_BIN) normal
- read(IIN_BIN) bc%ibulk1
- read(IIN_BIN) bc%ibulk2
- read(IIN_BIN) bc%coord(1,:)
- read(IIN_BIN) bc%coord(2,:)
- read(IIN_BIN) bc%coord(3,:)
- bc%dt = dt
+ allocate(bc%coord(3,(bc%nglob)))
+ read(IIN_BIN) ibool1
+ read(IIN_BIN) jacobian2Dw
+ read(IIN_BIN) normal
+ read(IIN_BIN) bc%ibulk1
+ read(IIN_BIN) bc%ibulk2
+ read(IIN_BIN) bc%coord(1,:)
+ read(IIN_BIN) bc%coord(2,:)
+ read(IIN_BIN) bc%coord(3,:)
+ bc%dt = dt
- allocate( bc%B(bc%nglob) )
- bc%B = 0e0_CUSTOM_REAL
- allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
- nx = 0e0_CUSTOM_REAL
- ny = 0e0_CUSTOM_REAL
- nz = 0e0_CUSTOM_REAL
- do e=1,bc%nspec
- do ij = 1,NGLLSQUARE
- k = ibool1(ij,e)
- nx(k) = nx(k) + normal(1,ij,e)
- ny(k) = ny(k) + normal(2,ij,e)
- nz(k) = nz(k) + normal(3,ij,e)
- bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
- 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
+end subroutine init_one_fault_TPV16
- !---------------------------------------------------------------------
+!---------------------------------------------------------------------
- subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+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), 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 / 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
- 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 (Ageing Law)
+ ! Set friction parameters and initialize friction variables
+ allocate( bc%rs )
+ allocate( bc%rs%V0(bc%nglob) )
+ allocate( bc%rs%f0(bc%nglob) )
+ allocate( bc%rs%a(bc%nglob) )
+ allocate( bc%rs%b(bc%nglob) )
+ allocate( bc%rs%L(bc%nglob) )
+ allocate( bc%rs%V_init(bc%nglob) )
+ allocate( bc%rs%theta_init(bc%nglob) )
+ ! WARNING: if V_HEALING is negative we turn off healing
- V0 =1.e-6_CUSTOM_REAL
- f0 =0.6_CUSTOM_REAL
- a =0.0080_CUSTOM_REAL
- b =0.0120_CUSTOM_REAL
- L =0.0135_CUSTOM_REAL
- V_init =1.e-12_CUSTOM_REAL
- theta_init =1.084207680000000e+09_CUSTOM_REAL
+ V0 =1.e-6_CUSTOM_REAL
+ f0 =0.6_CUSTOM_REAL
+ a =0.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=RS)
+ bc%rs%V0 = V0
+ bc%rs%f0 = f0
+ bc%rs%a = a
+ bc%rs%b = b
+ bc%rs%L = L
+ bc%rs%V_init = V_init
+ bc%rs%theta_init = theta_init
+ call init_2d_distribution(bc%rs%V0,bc%coord,IIN_PAR,nV0)
+ call init_2d_distribution(bc%rs%f0,bc%coord,IIN_PAR,nf0)
+ call init_2d_distribution(bc%rs%a,bc%coord,IIN_PAR,na)
+ call init_2d_distribution(bc%rs%b,bc%coord,IIN_PAR,nb)
+ call init_2d_distribution(bc%rs%L,bc%coord,IIN_PAR,nL)
+ call init_2d_distribution(bc%rs%V_init,bc%coord,IIN_PAR,nV_init)
+ call init_2d_distribution(bc%rs%theta_init,bc%coord,IIN_PAR,ntheta_init)
- bc%rs%theta = bc%rs%theta_init
- allocate(bc%MU(bc%nglob)) !Not necessary?
- !bc%MU = rs_mu(bc%rs)
- endif
+ bc%rs%theta = bc%rs%theta_init
+ allocate(bc%MU(bc%nglob)) !Not necessary?
+ bc%MU = rs_mu(bc%rs,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
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc,bc%nglob)
+ allocate( bc%asp )
+ allocate( bc%asp%AspTx(bc%nglob) )
+ allocate( bc%asp%Fload(bc%nglob) )
- end subroutine init_one_fault
+ AspTx =0.e0_CUSTOM_REAL
+ Fload =0.e0_CUSTOM_REAL
+ nAsp =0
+ nFload =0
- !---------------------------------------------------------------------
- subroutine compute_R(R,nglob,nx,ny,nz)
+ 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)
- integer :: nglob
- real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
- real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
+ endif
- 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
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc,bc%nglob)
- 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
+end subroutine init_one_fault
- 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
+!---------------------------------------------------------------------
+subroutine compute_R(R,nglob,nx,ny,nz)
- end subroutine compute_R
+ integer :: nglob
+ real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
+ real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
- !---------------------------------------------------------------------
- ! 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), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
- real(kind=CUSTOM_REAL), intent(inout) :: a(:)
- real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
- integer, intent(in) :: iin,n
+ ! 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
- 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
+ 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
- NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
+ 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
- if (n==0) return
+end subroutine compute_R
- 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
+!---------------------------------------------------------------------
+! adds a value to a fault parameter inside an area with prescribed shape
+subroutine init_2d_distribution(a,coord,iin,n)
- 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
+ real(kind=CUSTOM_REAL), intent(inout) :: a(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: iin,n
- where (b /= 0e0_CUSTOM_REAL) a = b
- enddo
+ 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
- end subroutine init_2d_distribution
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
- !---------------------------------------------------------------------
- elemental function heaviside(x)
+ if (n==0) return
- real(kind=CUSTOM_REAL), intent(in) :: x
- real(kind=CUSTOM_REAL) :: heaviside
+ 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
- if (x>=0e0_CUSTOM_REAL) then
- heaviside = 1e0_CUSTOM_REAL
- else
- heaviside = 0e0_CUSTOM_REAL
- endif
+ 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
- end function heaviside
+ where (b /= 0e0_CUSTOM_REAL) a = b
+ enddo
- !=====================================================================
- ! adds boundary term Bt into Force array for each fault.
- !
- subroutine bc_dynflt_set3d_all(F,Vel,Dis)
+end subroutine init_2d_distribution
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
- real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
+!---------------------------------------------------------------------
+elemental function heaviside(x)
- integer :: iflt
+ real(kind=CUSTOM_REAL), intent(in) :: x
+ real(kind=CUSTOM_REAL) :: heaviside
- 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 (x>=0e0_CUSTOM_REAL) then
+ heaviside = 1e0_CUSTOM_REAL
+ else
+ heaviside = 0e0_CUSTOM_REAL
+ endif
- end subroutine bc_dynflt_set3d_all
+end function heaviside
- !---------------------------------------------------------------------
- subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
+!=====================================================================
+! adds boundary term Bt into Force array for each fault.
+!
+subroutine bc_dynflt_set3d_all(F,Vel,Dis)
- use specfem_par, only:it,NSTEP
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
- 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
+ integer :: 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
+ 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
- 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,:))
+end subroutine bc_dynflt_set3d_all
- ! 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
+!---------------------------------------------------------------------
+subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
- ! rotate to fault frame (tangent,normal)
- ! component 3 is normal to the fault
- dD = rotate(bc,dD,1)
- dV = rotate(bc,dV,1)
- dA = rotate(bc,dA,1)
+ 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
- ! 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,:) )
+ 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
- !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
+ half_dt = 0.5e0_CUSTOM_REAL*bc%dt
+ theta_old = bc%swf%theta
+ Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- ! add initial stress
- T = T + bc%T0
+ ! get predicted values
+ dD = get_jump(bc,D) ! dD_predictor
+ dV = get_jump(bc,V) ! dV_predictor
+ dA = get_weighted_jump(bc,MxA) ! dA_free
- ! Solve for normal stress (negative is compressive)
- ! Opening implies free stress
- if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+ ! rotate to fault frame (tangent,normal)
+ ! component 3 is normal to the fault
+ dD = rotate(bc,dD,1)
+ dV = rotate(bc,dV,1)
+ dA = rotate(bc,dA,1)
- ! Update slip 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)
+ ! 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,:) )
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ !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
- ! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
+ ! add initial stress
+ T = T + bc%T0
- ! 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 normal stress (negative is compressive)
+ ! Opening implies free stress
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
- tnew = min(tnorm,strength)
- T(1,:) = tnew * t1
- T(2,:) = tnew * t2
+ ! 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)
- ! Save total tractions
- bc%T = T
+ ! Update friction coeficient
+ bc%MU = swf_mu_TPV16(bc%swf,it*bc%dt,bc%nglob)
- ! Subtract initial stress
- T = T - bc%T0
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
- ! 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 strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
- ! Update slip and slip rate, in fault frame
- bc%D = dD
- bc%V = dV + half_dt*dA
+ ! 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
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
+ tnew = min(tnorm,strength)
+ T(1,:) = tnew * t1
+ T(2,:) = tnew * t2
- ! 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,:)
+ ! Save total tractions
+ bc%T = T
- 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,:)
+ ! 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)
- !-- intermediate storage of outputs --
- Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
- Vnorm_old, Vnorm, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
+ ! Update 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)
- !-- 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)
+ ! 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,:)
- end subroutine BC_DYNFLT_set3d_TPV16
+ 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,:)
- !---------------------------------------------------------------------
- subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
- use specfem_par, only:it,NSTEP
+ !-- 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)
- 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
+ !-- 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)
- 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
+end subroutine BC_DYNFLT_set3d_TPV16
- 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,:))
+!---------------------------------------------------------------------
+subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
- ! 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
+ use specfem_par, only:it,NSTEP
- ! 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)
+ 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
- ! 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,:) )
+ 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
- !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
+ 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
- ! add initial stress
- T = T + bc%T0
+ 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,:))
- ! Solve for normal stress (negative is compressive)
- ! Opening implies free stress
- if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+ ! 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
- 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)
+ ! rotate to fault frame (tangent,normal)
+ ! component 3 is normal to the fault
+ dD = rotate(bc,dD,1)
+ dV = rotate(bc,dV,1)
+ dA = rotate(bc,dA,1)
- ! Update friction coeficient
- bc%MU = swf_mu(bc%swf)
+ ! 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
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ ! add initial stress
+ T = T + bc%T0
- ! Update strength
- strength = -bc%MU * 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)
- ! 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
+ 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)
- tnew = min(tnorm,strength)
- T(1,:) = tnew * t1
- T(2,:) = tnew * t2
+ ! Update friction coeficient
+ bc%MU = swf_mu(bc%swf)
- else ! Update rate and state friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
- Vf = sqrt(bc%V(1,:)**2 + bc%V(2,:)**2)
- call rs_update_state(Vf,bc%dt,bc%rs)
- tau_ratio = T(2,:)/T(1,:)
- tStick = sqrt(T(1,:)**2 + T(2,:)**2)
- ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
- call NRsearch(Vf1,tau1,bc%rs%f0,bc%rs%V0,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr,bc%nglob)
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
- ! 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
+ ! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
- 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
+ ! 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
- ! Updating state variable: 2nd loop
- Vfavg = 0.5e0_CUSTOM_REAL**abs(Vf + Vf1)
- call rs_update_state(Vfavg,bc%dt,bc%rs)
+ tnew = min(tnorm,strength)
+ T(1,:) = tnew * t1
+ T(2,:) = tnew * t2
- ! 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)
+ 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)
- where (Vf2 < 0._CUSTOM_REAL .or. tau2 < 0._CUSTOM_REAL .or. ierr == 1)
- tau2 = tStick
- Vf2 = 0._CUSTOM_REAL
- endwhere
+ ! 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
- ! 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
+ tStick = sqrt( (T(1,:) - bc%T0(1,:) + TxExt)**2 + T(2,:)**2)
+ ! tStick = sqrt(T(1,:)**2 + T(2,:)**2)
- where(T(2,:) > -VERYSMALLVAL)
- T(2,:) = tau2*abs(tau_ratio)/sqrt(ONE+tau_ratio**2) - bc%T0(2,:)
- elsewhere
- T(2,:) = -tau2*abs(tau_ratio)/sqrt(ONE+tau_ratio**2) - bc%T0(2,:)
- endwhere
+ 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)
- 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
+ ! 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
- endif
+ 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)
- ! Save total tractions
- bc%T = T
+ ! 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)
- ! Subtract initial stress
- T = T - bc%T0
+ where (Vf2 < 0._CUSTOM_REAL .or. tau2 < 0._CUSTOM_REAL .or. ierr == 1)
+ tau2 = tStick
+ Vf2 = 0._CUSTOM_REAL
+ endwhere
- ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
- dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
- dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
- dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
+ ! 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
- ! Update slip and slip rate, in fault frame
- bc%D = dD
- bc%V = dV + half_dt*dA
+ 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
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
+ 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
- ! 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,:)
+ endif
- 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,:)
+ ! Save total tractions
+ bc%T = T
- !-- 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)
+ ! Subtract initial stress
+ T = T - bc%T0
- !-- outputs --
- ! write dataT every NTOUT time step or at the end of simulation
- if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
- ! write dataXZ every NSNAP time step
- if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
- if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
+ ! Update 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)
- end subroutine BC_DYNFLT_set3d
+ ! Update slip and slip rate, in fault frame
+ bc%D = dD
+ bc%V = dV + half_dt*dA
- !===============================================================
- function get_jump (bc,v) result(dv)
+ ! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
- real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
+ ! 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,:)
- ! 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)
+ 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,:)
- end function get_jump
+ !-- 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)
- !---------------------------------------------------------------------
- function get_weighted_jump (bc,f) result(da)
+ !-- 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)
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
+end subroutine BC_DYNFLT_set3d
- real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
+!===============================================================
+function get_jump (bc,v) result(dv)
- ! 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)
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
+ real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
- end function get_weighted_jump
+ ! 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)
- !----------------------------------------------------------------------
- function rotate(bc,v,fb) result(vr)
+end function get_jump
- 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)
+!---------------------------------------------------------------------
+function get_weighted_jump (bc,f) result(da)
- ! Percy, tangential direction Vt, equation 7 of Pablo's notes in agreement with SPECFEM3D
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
- ! 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
+ real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
- ! 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
+ ! 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)
- endif
+end function get_weighted_jump
- end function rotate
+!----------------------------------------------------------------------
+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)
- !=====================================================================
- subroutine swf_update_state(dold,dnew,vold,f)
+ ! Percy, tangential direction Vt, equation 7 of Pablo's notes in agreement with SPECFEM3D
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
- type(swf_type), intent(inout) :: f
+ ! 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
- real(kind=CUSTOM_REAL) :: vnorm
- integer :: k,npoin
+ ! 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
- f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
+ 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
+end function rotate
- !=====================================================================
- subroutine rs_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
+!=====================================================================
+subroutine swf_update_state(dold,dnew,vold,f)
- real(kind=CUSTOM_REAL) :: vDtL(size(V))
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
+ type(swf_type), intent(inout) :: f
- vDtL = V * dt * f%L
+ real(kind=CUSTOM_REAL) :: vnorm
+ integer :: k,npoin
- 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
+ f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
- end subroutine rs_update_state
+ 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
- !=====================================================================
- ! Friction coefficient
- function swf_mu_TPV16(f,time,nglob) result(mu)
+!=====================================================================
+subroutine rs_update_state(V,dt,f)
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- real(kind=CUSTOM_REAL) :: time
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ type(rs_type), intent(inout) :: f
+ real(kind=CUSTOM_REAL), intent(in) :: dt
- integer :: i,nglob
+ real(kind=CUSTOM_REAL) :: vDtL(size(V))
- !-- 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
+ vDtL = V * dt * f%L
- end function swf_mu_TPV16
+ 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
- !=====================================================================
- ! Friction coefficient
- function swf_mu(f) result(mu)
+end subroutine rs_update_state
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+!=====================================================================
+! Friction coefficient
+function swf_mu_TPV16(f,time,nglob) result(mu)
- !-- linear slip weakening:
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- mu = max( mu, f%mud)
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+ real(kind=CUSTOM_REAL) :: time
- end function swf_mu
+ integer :: i,nglob
- !=====================================================================
- ! Friction coefficient
- function rs_mu(f,V) result(mu)
+ !-- 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
- type(rs_type), intent(in) :: f
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
- real(kind=CUSTOM_REAL) :: mu(size(V,2))
+end function swf_mu_TPV16
- !-- 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))
+!=====================================================================
+! Friction coefficient
+function swf_mu(f) result(mu)
- end function rs_mu
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- !===============================================================
- ! OUTPUTS
- subroutine init_dataT(DataT,coord,nglob,NT,iflt)
- ! NT = total number of time steps
+ !-- linear slip weakening:
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ mu = max( mu, f%mud)
- integer, intent(in) :: nglob,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
- type (dataT_type), intent(out) :: DataT
+end function swf_mu
- real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
- integer :: i, iglob , IIN, ier, jflt, np, k
- character(len=70) :: tmpname
+!=====================================================================
+! Friction coefficient
+function rs_mu(f,V) result(mu)
- ! 1. read fault output coordinates from user file,
- ! 2. define iglob: the fault global index of the node nearest to user
- ! requested coordinate
+ 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))
- 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)
+ !-- 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))
- if (DataT%npoin == 0) return
+end function rs_mu
- allocate(DataT%iglob(DataT%npoin))
- allocate(DataT%name(DataT%npoin))
+!===============================================================
+! OUTPUTS
+subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+ ! NT = total number of time steps
- 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)
+ integer, intent(in) :: nglob,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
+ type (dataT_type), intent(out) :: DataT
- 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
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
- ! 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
+ ! 1. read fault output coordinates from user file,
+ ! 2. define iglob: the fault global index of the node nearest to user
+ ! requested coordinate
- 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)
- end subroutine init_dataT
+ if (DataT%npoin == 0) return
- !---------------------------------------------------------------
- subroutine store_dataT(dataT,d,v,t,itime)
+ allocate(DataT%iglob(DataT%npoin))
+ allocate(DataT%name(DataT%npoin))
- type(dataT_type), intent(inout) :: dataT
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
- integer, intent(in) :: itime
+ 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)
- integer :: i,k
+ 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 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
+ ! 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
- end subroutine store_dataT
+ close(IIN)
- !-----------------------------------------------------------------
- subroutine write_dataT_all(nt)
+end subroutine init_dataT
- integer, intent(in) :: nt
+!---------------------------------------------------------------
+subroutine store_dataT(dataT,d,v,t,itime)
- integer :: i
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ integer, intent(in) :: itime
- if (.not.allocated(faults)) return
- do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
- enddo
+ integer :: i,k
- end subroutine write_dataT_all
+ 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
- !------------------------------------------------------------------------
- subroutine SCEC_write_dataT(dataT,DT,NT)
+end subroutine store_dataT
- type(dataT_type), intent(in) :: dataT
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT
+!-----------------------------------------------------------------
+subroutine write_dataT_all(nt)
- integer :: i,k,IOUT
- character :: NTchar*5
- integer :: today(3), now(3)
+ integer, intent(in) :: nt
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+ integer :: i
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+ if (.not.allocated(faults)) return
+ do i = 1,size(faults)
+ call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
+ enddo
- write(NTchar,1) NT
- NTchar = adjustl(NTchar)
+end subroutine write_dataT_all
-1 format(I5)
- do i=1,dataT%npoin
+!------------------------------------------------------------------------
+subroutine SCEC_write_dataT(dataT,DT,NT)
- 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)
+ 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)
-1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+ 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
- enddo
+ write(NTchar,1) NT
+ NTchar = adjustl(NTchar)
- end subroutine SCEC_write_dataT
+1 format(I5)
+ do i=1,dataT%npoin
- !-------------------------------------------------------------------------------------------------
- subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
+ 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)
- 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)
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+ enddo
- write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
+end subroutine SCEC_write_dataT
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+!-------------------------------------------------------------------------------------------------
+subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
- 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
+ type(dataXZ_type), intent(in) :: dataXZ
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT,iflt
- close(IOUT)
+ integer :: i,IOUT
+ character(len=70) :: filename
+ integer*4 today(3), now(3)
-1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
- end subroutine SCEC_Write_RuptureTime
- !-------------------------------------------------------------------------------------------------
- subroutine init_dataXZ(DataXZ,bc,nglob)
+ write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
- type(dataXZ_type), intent(inout) :: DataXZ
- type(bc_dynflt_type) :: bc
- integer, intent(in) :: nglob
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
- 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))
+ 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
- !Percy , setting up initial rupture time null for all faults.
- DataXZ%tRUP = 0e0_CUSTOM_REAL
- DataXZ%tPZ = 0e0_CUSTOM_REAL
+ close(IOUT)
- end subroutine init_dataXZ
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
- !---------------------------------------------------------------
- subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
+end subroutine SCEC_Write_RuptureTime
- 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
+!-------------------------------------------------------------------------------------------------
+subroutine init_dataXZ(DataXZ,bc,nglob)
- integer :: i
+ type(dataXZ_type), intent(inout) :: DataXZ
+ type(bc_dynflt_type) :: bc
+ integer, intent(in) :: nglob
- ! "stg" : strength .
+ 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))
- dataXZ%stg = stg
+ !Percy , setting up initial rupture time null for all faults.
+ DataXZ%tRUP = 0e0_CUSTOM_REAL
+ DataXZ%tPZ = 0e0_CUSTOM_REAL
- do i = 1,size(stg)
- ! process zone time = first time when slip = dc (break down process).
- ! with linear time interpolation
- if (dataXZ%tPZ(i)==0e0_CUSTOM_REAL) then
- if (dold(i)<=dc(i) .and. dnew(i) >= dc(i)) then
- dataXZ%tPZ(i) = time-dt*(dnew(i)-dc(i))/(dnew(i)-dold(i))
- endif
- endif
- ! rupture time = first time when slip velocity = vc
- ! with linear time interpolation
- ! vc should be pre-defined as input data .
+end subroutine init_dataXZ
- if (dataXZ%tRUP(i)==0e0_CUSTOM_REAL) then
- if (vold(i)<=V_RUPT .and. vnew(i)>=V_RUPT) dataXZ%tRUP(i)= time-dt*(vnew(i)-V_RUPT)/(vnew(i)-vold(i))
- endif
- enddo
+!---------------------------------------------------------------
+subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
- ! To do : add stress criteria (firs time strength is reached).
+ 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
- ! note: the other arrays in dataXZ are pointers to arrays in bc
- ! they do not need to be updated here
+ integer :: i
- end subroutine store_dataXZ
+ ! "stg" : strength .
- !---------------------------------------------------------------
- subroutine write_dataXZ(dataXZ,itime,iflt)
+ dataXZ%stg = stg
- type(dataXZ_type), intent(in) :: dataXZ
- integer, intent(in) :: itime,iflt
+ 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 .
- character(len=70) :: filename
+ 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
- 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.
+ ! To do : add stress criteria (firs time strength is reached).
- ! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
+ ! note: the other arrays in dataXZ are pointers to arrays in bc
+ ! they do not need to be updated here
- open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
+end subroutine store_dataXZ
- 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)
+!---------------------------------------------------------------
+subroutine write_dataXZ(dataXZ,itime,iflt)
- end subroutine write_dataXZ
+ type(dataXZ_type), intent(in) :: dataXZ
+ integer, intent(in) :: itime,iflt
- !---------------------------------------------------------------
+ character(len=70) :: filename
- ! 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)
+ 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.
- !implicit none
- !include "constants.h"
+ ! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
- integer :: iglob,nglob
+ open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
- 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
+ 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
- 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
+! 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)
- 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
+ !implicit none
+ !include "constants.h"
- END DO
+ integer :: iglob,nglob
- !compute the updated Vf before exiting NR search
+ 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)
- enddo
- end subroutine NRsearch
+ 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