[cig-commits] r19348 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Tue Jan 10 21:52:51 PST 2012
Author: surendra
Date: 2012-01-10 21:52:51 -0800 (Tue, 10 Jan 2012)
New Revision: 19348
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Added nucleation procedure for TPV16/17
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-01-10 01:20:11 UTC (rev 19347)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-01-11 05:52:51 UTC (rev 19348)
@@ -59,7 +59,7 @@
private
integer :: kind
logical :: healing = .false.
- real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), theta=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), theta=>null(), T=>null()
end type swf_type
type bc_dynflt_type
@@ -91,6 +91,8 @@
logical, save :: SIMULATION_TYPE_DYN = .false.
+ logical, save :: TPV16 = .false.
+
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
@@ -165,7 +167,11 @@
allocate( faults(nbfaults) )
do iflt=1,nbfaults
read(IIN_PAR,nml=BEGIN_FAULT,end=100)
- call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
+ 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)
@@ -191,6 +197,170 @@
!---------------------------------------------------------------------
+ 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
+
+ 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 :: 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 :: iLoc
+
+ 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))
+
+ allocate(bc%coord(3,(bc%nglob)))
+ read(IIN_BIN) ibool1
+ read(IIN_BIN) jacobian2Dw
+ read(IIN_BIN) normal
+ read(IIN_BIN) bc%ibulk1
+ read(IIN_BIN) bc%ibulk2
+ read(IIN_BIN) bc%coord(1,:)
+ read(IIN_BIN) bc%coord(2,:)
+ read(IIN_BIN) bc%coord(3,:)
+ bc%dt = dt
+
+ allocate( bc%B(bc%nglob) )
+ bc%B = 0e0_CUSTOM_REAL
+ allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
+ nx = 0e0_CUSTOM_REAL
+ ny = 0e0_CUSTOM_REAL
+ nz = 0e0_CUSTOM_REAL
+ do e=1,bc%nspec
+ do ij = 1,NGLLSQUARE
+ k = ibool1(ij,e)
+ nx(k) = nx(k) + normal(1,ij,e)
+ ny(k) = ny(k) + normal(2,ij,e)
+ nz(k) = nz(k) + normal(3,ij,e)
+ bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
+ enddo
+ enddo
+ do k=1,bc%nglob
+ norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
+ nx(k) = nx(k) / norm
+ ny(k) = ny(k) / norm
+ nz(k) = nz(k) / norm
+ enddo
+
+ allocate( bc%R(3,3,bc%nglob) )
+ call compute_R(bc%R,bc%nglob,nx,ny,nz)
+
+! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
+ allocate(bc%invM1(bc%nglob))
+ allocate(bc%invM2(bc%nglob))
+ bc%invM1 = Minv(bc%ibulk1)
+ bc%invM2 = Minv(bc%ibulk2)
+
+! Fault impedance, Z in : Trac=T_Stick-Z*dV
+! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
+! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
+! NOTE: same Bi on both sides, see note above
+ allocate(bc%Z(bc%nglob))
+ bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
+
+ allocate(bc%T(3,bc%nglob))
+ allocate(bc%D(3,bc%nglob))
+ allocate(bc%V(3,bc%nglob))
+ bc%T = 0e0_CUSTOM_REAL
+ bc%D = 0e0_CUSTOM_REAL
+ bc%V = 0e0_CUSTOM_REAL
+
+
+ allocate(inp_nx(bc%nglob))
+ allocate(inp_nz(bc%nglob))
+ allocate(loc_str(bc%nglob))
+ allocate(loc_dip(bc%nglob))
+ allocate(sigma0(bc%nglob))
+ allocate(tau0_str(bc%nglob))
+ allocate(tau0_dip(bc%nglob))
+ allocate(Rstress_str(bc%nglob))
+ allocate(Rstress_dip(bc%nglob))
+ allocate(static_fc(bc%nglob))
+ allocate(dyn_fc(bc%nglob))
+ allocate(swcd(bc%nglob))
+ allocate(cohes(bc%nglob))
+ allocate(tim_forcedRup(bc%nglob))
+
+ open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
+ read(IIN_NUC,*) relz_num,sub_relz_num
+ read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
+ read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ do ipar=1,bc%nglob
+ read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
+Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
+ enddo
+ close(IIN_NUC)
+
+ allocate(bc%T0(3,bc%nglob))
+ allocate( bc%swf )
+ allocate( bc%swf%mus(bc%nglob) )
+ allocate( bc%swf%mud(bc%nglob) )
+ allocate( bc%swf%Dc(bc%nglob) )
+ allocate( bc%swf%theta(bc%nglob) )
+ allocate( bc%swf%T(bc%nglob) )
+
+ allocate(iLoc(bc%nglob))
+ do ipar=1,bc%nglob
+ iLoc = minloc( (loc_str(ipar)-bc%coord(1,:))**2 + (-loc_dip(ipar)-bc%coord(3,:))**2 ) !iloc_dip is negative of Z-coord
+print*,minval(iLoc)
+ bc%T0(1,iLoc(1)) = sigma0(ipar)
+ bc%T0(2,iLoc(1)) = tau0_str(ipar)
+ bc%T0(3,iLoc(1)) = tau0_dip(ipar)
+
+ bc%swf%mus(iLoc(1)) = static_fc(ipar)
+ bc%swf%mud(iLoc(1)) = dyn_fc(ipar)
+ bc%swf%Dc(iLoc(1)) = swcd(ipar)
+ bc%swf%T(iLoc(1)) = tim_forcedRup(ipar)
+
+enddo
+
+ ! 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)
+
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc,bc%nglob)
+
+
+end subroutine init_one_fault_TPV16
+
+!---------------------------------------------------------------------
+
subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
type(bc_dynflt_type), intent(inout) :: bc
@@ -329,10 +499,6 @@
call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
call init_dataXZ(bc%dataXZ,bc,bc%nglob)
-!---- Initial Fault Parameters -------------
- call write_init_fault_par(bc,iflt)
-
-
end subroutine init_one_fault
!---------------------------------------------------------------------
@@ -403,20 +569,14 @@
read(iin,DIST2D)
select case(shape)
case ('circle')
- b = heaviside(r-sqrt((coord(1,:)-xc)**2 + &
- (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2)) * val
+ 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)) * val
+ 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 ('cilinder')
- b = heaviside(r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2)) * &
- heaviside((lz/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) * &
@@ -462,12 +622,139 @@
if (.not. allocated(faults)) return
do iflt=1,size(faults)
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
+ if(.not.TPV16) then
+ if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
+ else
+ if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d_TPV16(faults(iflt),F,Vel,Dis,iflt)
+ endif
enddo
end subroutine bc_dynflt_set3d_all
!---------------------------------------------------------------------
+ subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
+
+ use specfem_par, only:it,NSTEP
+
+ real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
+ type(bc_dynflt_type), intent(inout) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
+ integer,intent(in) :: iflt
+
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
+ real(kind=CUSTOM_REAL) :: half_dt
+! integer :: k
+
+ half_dt = 0.5e0_CUSTOM_REAL*bc%dt
+ theta_old = bc%swf%theta
+ Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+
+! get predicted values
+ dD = get_jump(bc,D) ! dD_predictor
+ dV = get_jump(bc,V) ! dV_predictor
+ dA = get_weighted_jump(bc,MxA) ! dA_free
+
+! rotate to fault frame (tangent,normal)
+! component 3 is normal to the fault
+ dD = rotate(bc,dD,1)
+ dV = rotate(bc,dV,1)
+ dA = rotate(bc,dA,1)
+
+
+! T_stick
+ T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
+ T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
+ T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
+
+!Warning : dirty particular free surface condition z = 0.
+! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
+! do k=1,bc%nglob
+! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
+! end do
+
+! add initial stress
+ T = T + bc%T0
+
+! Solve for normal stress (negative is compressive)
+ ! Opening implies free stress
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+
+! 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)
+
+! combined with time-weakening for nucleation
+! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+
+! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
+
+! Solve for shear stress
+ tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnorm = max(tnorm,1e0_CUSTOM_REAL)
+ t1 = T(1,:)/tnorm
+ t2 = T(2,:)/tnorm
+
+ tnew = min(tnorm,strength)
+ T(1,:) = tnew * t1
+ T(2,:) = tnew * t2
+
+! Save total tractions
+ bc%T = T
+
+! Subtract initial stress
+ T = T - bc%T0
+
+! Update slip acceleration da=da_free-T/(0.5*dt*Z)
+ dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
+ dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
+ dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
+
+! Update slip and slip rate, in fault frame
+ bc%D = dD
+ bc%V = dV + half_dt*dA
+
+! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
+
+! Add boundary term B*T to M*a
+ MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
+ MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
+ MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
+
+ MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
+ MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
+ MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
+
+
+!-- intermediate storage of outputs --
+ Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
+ Vnorm_old, Vnorm, it*bc%dt,bc%dt)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
+
+
+!-- outputs --
+! write dataT every NTOUT time step or at the end of simulation
+ if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
+! write dataXZ every NSNAP time step
+ if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
+ if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
+
+ end subroutine BC_DYNFLT_set3d_TPV16
+
+!===============================================================
+
+
+!---------------------------------------------------------------------
subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
use specfem_par, only:it,NSTEP
@@ -662,7 +949,29 @@
endif
end subroutine swf_update_state
+!=====================================================================
+! Friction coefficient
+ function swf_mu_TPV16(f,time,nglob) result(mu)
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+ real(kind=CUSTOM_REAL) :: time
+
+ integer :: i,nglob
+
+ !-- linear slip weakening:
+
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ do i=1,nglob
+ if ( (f%theta(i) >= f%Dc(i)) .or. (time >= f%T(i)) ) then
+ mu(i) = f%mud(i)
+ endif
+ enddo
+
+ end function swf_mu_TPV16
+
+
+
!=====================================================================
! Friction coefficient
function swf_mu(f) result(mu)
@@ -802,6 +1111,13 @@
integer :: i,k,IOUT
character :: NTchar*5
+
+ integer*4 today(3), now(3)
+
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+
+
IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
write(NTchar,1) NT
@@ -811,16 +1127,14 @@
do i=1,dataT%npoin
open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "# problem=TPV15"
- write(IOUT,*) "# author=Galvez, Ampuero, Nissen-Meyer"
- write(IOUT,*) "# date=2011/xx/xx"
- write(IOUT,*) "# code=SPECFEM3D_FAULT "
+ write(IOUT,*) "# problem=TPV16"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ write(IOUT,1000) "# date=",today(2), today(1), today(3), now
+ write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=100 m (*4 GLL nodes)"
+ write(IOUT,*) "# element_size=75 m (*3 GLL nodes)"
write(IOUT,*) "# time_step=",DT
- write(IOUT,*) "# num_time_steps=",NT
write(IOUT,*) "# location=",trim(dataT%name(i))
- write(IOUT,*) "# Time series in 8 column of E15.7"
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)"
@@ -839,6 +1153,12 @@
dataT%t3(k,i)/1.0e6_CUSTOM_REAL
enddo
close(IOUT)
+
+
+ 1000 format ( 'Date ', i2.2, '/', i2.2, '/', i4.4, '; time ',
+ & i2.2, ':', i2.2, ':', i2.2 )
+
+
enddo
end subroutine SCEC_write_dataT
@@ -854,19 +1174,23 @@
integer :: i,IOUT
character(len=70) :: filename
+ integer*4 today(3), now(3)
+
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+
+
write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
open(IOUT,file=trim(filename),status='replace')
- write(IOUT,*) "# problem=TPV5"
- write(IOUT,*) "# author=Galvez, Ampuero, Tarje"
- write(IOUT,*) "# date=2011/xx/xx"
- write(IOUT,*) "# code=SPECFEM3D_FAULT"
+ write(IOUT,*) "# problem=TPV16"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ write(IOUT,1000) "# date=",today(2), today(1), today(3), now
+ write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=100 m (*4 GLL nodes)"
- write(IOUT,*) "# time_step=",DT
- write(IOUT,*) "# num_time_steps=",NT
+ write(IOUT,*) "# element_size=75 m (*3 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)"
@@ -877,6 +1201,10 @@
close(IOUT)
+ 1000 format ( 'Date ', i2.2, '/', i2.2, '/', i4.4, '; time ',
+ & i2.2, ':', i2.2, ':', i2.2 )
+
+
end subroutine SCEC_Write_RuptureTime
!-------------------------------------------------------------------------------------------------
@@ -983,39 +1311,5 @@
end subroutine write_dataXZ
-!---------------------------------------------------------------
- subroutine write_init_fault_par(bc,iflt)
-
- type(bc_dynflt_type), intent(inout) :: bc
- integer, intent(in) :: iflt
-
- character(len=70) :: filename
-
-
- write(filename,"('OUTPUT_FILES/Fault_init_par','_F',I0,'.bin')") iflt
-! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
-! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
-! compilers.
-
-! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
-
- open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
-
- write(IOUT) bc%coord(1,:)
- write(IOUT) bc%coord(2,:)
- write(IOUT) bc%coord(3,:)
- write(IOUT) bc%swf%mus
- write(IOUT) bc%swf%mud
- write(IOUT) bc%swf%Dc
- write(IOUT) bc%T0(1,:)
- write(IOUT) bc%T0(2,:)
- write(IOUT) bc%T0(3,:)
-
- close(IOUT)
-
- end subroutine write_init_fault_par
-
-!---
-
end module fault_solver
More information about the CIG-COMMITS
mailing list