[cig-commits] r19927 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Sun Apr 8 13:29:40 PDT 2012
Author: surendra
Date: 2012-04-08 13:29:40 -0700 (Sun, 08 Apr 2012)
New Revision: 19927
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Added Boxcar distribution. Removed NRsearch and switched to traditional Newton-Raphson solver. Reworked fault outputs to include rate and state parameters. Cleaned up RS solver
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-04-06 02:56:24 UTC (rev 19926)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-04-08 20:29:40 UTC (rev 19927)
@@ -41,7 +41,7 @@
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
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1,v1,t1,d2,v2,t2,t3,theta
character(len=70), dimension(:), pointer :: name
end type dataT_type
@@ -230,8 +230,8 @@
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) :: Fload
- integer :: nFload
+ real(kind=CUSTOM_REAL) :: AspTx,Fload
+ integer :: nAsp,nFload
real(kind=CUSTOM_REAL) :: C,T
integer :: nC,nForcedRup
@@ -254,6 +254,16 @@
real(kind=CUSTOM_REAL) :: minX
+
+ real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
+ real(kind=CUSTOM_REAL) :: x,z
+ logical :: c1,c2,c3,c4
+ real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
+ integer iglob
+
+ double precision, dimension(:), allocatable :: mu1,mu2,mu3
+
+
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc,C,T,nC,nForcedRup
NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init,C,T,nC,nForcedRup
@@ -315,9 +325,6 @@
! NOTE: same Bi on both sides, see note above
allocate(bc%Z(bc%nglob))
bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
- ! WARNING: In non-split nodes at fault edges M is assembled across the fault.
- ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
- ! In a symmetric mesh (M1=M2) Z will be twice its intended value
allocate(bc%T(3,bc%nglob))
allocate(bc%D(3,bc%nglob))
@@ -499,11 +506,58 @@
call init_2d_distribution(bc%rsf%C ,bc%coord,IIN_PAR,nC)
call init_2d_distribution(bc%rsf%T ,bc%coord,IIN_PAR,nForcedRup)
+
+ W1=15000._CUSTOM_REAL
+ W2=7500._CUSTOM_REAL
+ w=3000._CUSTOM_REAL
+ hypo_z = -7500._CUSTOM_REAL
+ do iglob=1,bc%nglob
+ x=bc%coord(1,iglob)
+ z=bc%coord(3,iglob)
+ c1=abs(x)<W1+w
+ c2=abs(x)>W1
+ c3=abs(z-hypo_z)<W2+w
+ c4=abs(z-hypo_z)>W2
+ if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
+
+ if (c1 .and. c2) then
+ b11 = w/(abs(x)-W1-w)
+ b12 = w/(abs(x)-W1)
+ B1 = 0.5 * (ONE + tanh(b11 + b12))
+ elseif(abs(x)<=W1) then
+ B1 = 1._CUSTOM_REAL
+ else
+ B1 = 0._CUSTOM_REAL
+ endif
+
+ if (c3 .and. c4) then
+ b21 = w/(abs(z-hypo_z)-W2-w)
+ b22 = w/(abs(z-hypo_z)-W2)
+ B2 = 0.5 * (ONE + tanh(b21 + b22))
+ elseif(abs(z-hypo_z)<=W2) then
+ B2 = 1._CUSTOM_REAL
+ else
+ B2 = 0._CUSTOM_REAL
+ endif
+
+
+ bc%rsf%a(iglob) = 0.008 + 0.008 * (ONE - B1*B2)
+ elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
+ bc%rsf%a(iglob) = 0.008
+ else
+ bc%rsf%a(iglob) = 0.016
+ endif
+
+ bc%rsf%theta(iglob) = L/V0 * exp( (bc%rsf%a(iglob) * log(2.0*sinh(-S1/S3/bc%rsf%a(iglob))) - f0 - bc%rsf%a(iglob)*log(V_init/V0) )/b )
+
+
+ enddo
+
+
+
allocate(bc%MU(bc%nglob))
- bc%MU = rsf_mu(bc%rsf,bc%rsf%V_init) ! Should we switch to regularized form in rsf_mu?
- bc%T0(1,:) = bc%MU * bc%T0(3,:)
- bc%T0(2,:) = 0._CUSTOM_REAL
+ bc%V(1,:) = bc%rsf%V_init
allocate( bc%asp )
allocate( bc%asp%Fload(bc%nglob) )
@@ -516,6 +570,7 @@
endif
+
call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
call init_dataXZ(bc%dataXZ,bc,bc%nglob)
@@ -560,7 +615,6 @@
!---------------------------------------------------------------------
! adds a value to a fault parameter inside an area with prescribed shape
subroutine init_2d_distribution(a,coord,iin,n)
-!JPA refactor: background value shuld be an argument
real(kind=CUSTOM_REAL), intent(inout) :: a(:)
real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
@@ -595,7 +649,7 @@
case ('circle-exp')
r1 = sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 )
where(r1<r)
- b =exp (r1**2/(r1**2 - r**2) )
+ b =exp(r1**2/(r1**2 - r**2) )
elsewhere
b =0._CUSTOM_REAL
endwhere
@@ -642,9 +696,6 @@
!=====================================================================
! adds boundary term Bt into Force array for each fault.
!
-! NOTE: On non-split nodes at fault edges, dD=dV=dA=0
-! and the net contribution of B*T is =0
-!
subroutine bc_dynflt_set3d_all(F,Vel,Dis)
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
@@ -756,29 +807,14 @@
TxExt = DTau0 * bc%asp%Fload * GLoad
T(1,:) = T(1,:) + TxExt
-!JPA the solver below can be refactored into a loop with two passes
Vf = Vnorm_old
theta_old = bc%rsf%theta
call rsf_update_state(Vf,bc%dt,bc%rsf)
- ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
tStick = sqrt(T(1,:)**2 + T(2,:)**2)
do iNode=1,bc%nglob
- psi = log( bc%rsf%theta(iNode) * bc%rsf%V_init(iNode) / bc%rsf%L(iNode) )
- call NRsearch(Vf1(iNode),tau1(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),-T(3,iNode),psi,bc%Z(iNode),ta(iNode),tStick(iNode),ierr)
-
- ! If Vf1 or tau1 is negative, set Vf1 = 0 and tau1 = TauStick
- if (Vf1(iNode) < 0._CUSTOM_REAL .or. tau1(iNode) < 0._CUSTOM_REAL .or. ierr == 1) then
- tau1(iNode) = tStick(iNode)
- Vf1(iNode) = 0._CUSTOM_REAL
- endif
-
- if (ierr == 1) then
- print*, 'NR search 1st loop failed!'
- print*, 'iteration = ', it
- !call exit_MPI(myrank,'NR search 1st loop failed!')
- stop
- endif
+ Vf1(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
+ tau1(iNode) = tStick(iNode) - bc%Z(iNode)*Vf1(iNode)
enddo
! Updating state variable: 2nd loop
@@ -788,20 +824,8 @@
! NR search 2nd loop
do iNode=1,bc%nglob
- psi = log( bc%rsf%theta(iNode) * bc%rsf%V_init(iNode) / bc%rsf%L(iNode) )
- call NRsearch(Vf2(iNode),tau2(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),-T(3,iNode),psi,bc%Z(iNode),ta(iNode),tStick(iNode),ierr)
-
- if (Vf2(iNode) < 0._CUSTOM_REAL .or. tau2(iNode) < 0._CUSTOM_REAL .or. ierr == 1) then
- tau2(iNode) = tStick(iNode)
- Vf2(iNode) = 0._CUSTOM_REAL
- endif
-
- if (ierr == 1) then
- print*, 'NR search 2nd loop failed!'
- print*, 'iteration = ', it
- !call exit_MPI(myrank,'NR search 2nd loop failed!')
- stop
- endif
+ Vf2(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
+ tau2(iNode) = tStick(iNode) - bc%Z(iNode)*Vf2(iNode)
enddo
tStick = max(tStick,1e0_CUSTOM_REAL)
@@ -849,7 +873,7 @@
call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
Vnorm_old, Vnorm, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
!-- outputs --
! write dataT every NTOUT time step or at the end of simulation
@@ -887,9 +911,6 @@
da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
- ! NOTE: In non-split nodes at fault edges M and f are assembled across the fault.
- ! Hence, f1=f2, invM1=invM2=1/(M1+M2) instead of invMi=1/Mi, and da=0.
-
end function get_weighted_jump
!----------------------------------------------------------------------
@@ -955,17 +976,82 @@
function rsf_mu(f,V) result(mu)
type(rsf_type), intent(in) :: f
- !real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
- !real(kind=CUSTOM_REAL) :: mu(size(V,2))
real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
real(kind=CUSTOM_REAL) :: mu(size(V))
+ mu = f%a * asinh( V/2.0/f%V0 * exp((f%f0 + f%b*log(f%theta*f%V0/f%L))/f%a ) ) ! Regularized
- !mu = f%f0 + (f%a * log(sqrt(V(1,:)**2 + V(2,:)**2)/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
- mu = f%f0 + (f%a * log(V/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
-
end function rsf_mu
!---------------------------------------------------------------------
+subroutine funcd(x,fn,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+
+ real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
+ double precision :: arg,fn,df,x
+
+ arg = exp((f0+dble(b)*log(V0*theta/L))/a)/2._CUSTOM_REAL/V0
+ fn = tStick - Z*x - a*Seff*asinh(x*arg)
+ df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
+end subroutine funcd
+
+!---------------------------------------------------------------------
+function rtsafe(funcd,x1,x2,xacc,tStick,Seff,Z,f0,V0,a,b,L,theta)
+
+ integer, parameter :: MAXIT=200
+ real(kind=CUSTOM_REAL) :: x1,x2,xacc
+ EXTERNAL funcd
+ integer :: j
+ !real(kind=CUSTOM_REAL) :: df,dx,dxold,f,fh,fl,temp,xh,xl
+ double precision :: df,dx,dxold,f,fh,fl,temp,xh,xl,rtsafe
+ real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
+
+ call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ if( (fl>0 .and. fh>0) .or. (fl<0 .and. fh<0) ) stop 'root must be bracketed in rtsafe'
+ if(fl==0.) then
+ rtsafe=x2
+ return
+ elseif(fh==0.) then
+ rtsafe=x2
+ return
+ elseif(fl<0) then
+ xl=x1
+ xh=x2
+ else
+ xh=x1
+ xl=x2
+ endif
+
+ rtsafe=0.5d0*(x1+x2)
+ dxold=abs(x2-x1)
+ dx=dxold
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ do j=1,MAXIT
+ if( ((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f)>0 .or. abs(2.*f)>abs(dxold*df) ) then
+ dxold=dx
+ dx=0.5d0*(xh-xl)
+ rtsafe=xl+dx
+ if(xl==rtsafe) return
+ else
+ dxold=dx
+ dx=f/df
+ temp=rtsafe
+ rtsafe=rtsafe-dx
+ if(temp==rtsafe) return
+ endif
+ if(abs(dx)<xacc) return
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ if(f<0.) then
+ xl=rtsafe
+ else
+ xh=rtsafe
+ endif
+ enddo
+ stop 'rtsafe exceeding maximum iterations'
+ return
+
+end function rtsafe
+
+!---------------------------------------------------------------------
subroutine rsf_update_state(V,dt,f)
real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
@@ -991,91 +1077,6 @@
end subroutine rsf_update_state
-!---------------------------------------------------------------
-! Newton-Raphson search for rate-and-state solver
-!
-! NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror)
-! Newton-Raphson search used to find slip rates
-!
-! When Newton-Raphson search fails (happens near the free surface points),
-! use bisection method which will not fail to find the root if the specified
-! interval contains a root.
-! NOTE: When an interval contains more than one root, the bisection method
-! can find one of them.
-!
-!
-! INPUT VARIABLES fo,Vo,cca,ccb,Seff,tau,FaultVFree,FaultZ,ierror
-!
-! OUTPUT VARIABLES Vfout,tauout
-!
-! Yoshihiro Kaneko: ykaneko at gps.caltech.edu
-subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror)
-
- real(kind=CUSTOM_REAL) :: cA,eps,delta,fa,help,help1,help2,Vfprime
- real(kind=CUSTOM_REAL) :: Vfout,tauout,fo,Vo,cca,ccb,Seff,tau,psi,FaultZ,TauStick
- real(kind=CUSTOM_REAL) :: f_bisec,f_bisec_Vflow,Vfbisec,Vfhigh,Vflow
- integer :: j,k
- integer :: ierror
-
-
- cA = cca*Seff
- eps = 0.001*cA
- k=0
- delta = HUGEVAL
-
- DO WHILE ( abs(delta) > eps )
- fa = tau/(Seff*cca)
- help = -(fo+ccb*psi)/cca
- help1 = exp(help+fa)
- help2 = exp(help-fa)
- Vfout = Vo*(help1 - help2)
- Vfprime = (Vo/cA)*(help1 + help2)
- delta = (TauStick - HALF*FaultZ*Vfout - tau)/(ONE + HALF*FaultZ*Vfprime)
- tau = tau + delta
- k = k + 1
- if ( abs(delta) > 1.e10 .OR. k > 100 ) then
- ! try bisection if NR search fails
- Vflow = 0._CUSTOM_REAL ! lower bound of Vf for bisection
- Vfhigh = 5._CUSTOM_REAL ! upper bound of Vf for bisection
- do j = 1,200
- Vfbisec = (Vflow + Vfhigh)/TWO
- f_bisec = TauStick - HALF*FaultZ*Vfbisec - cA*asinh(Vfbisec/(TWO*Vo)*exp(-help))
- if ( abs(f_bisec) < eps ) then
- Vfout = Vfbisec
- tauout = TauStick - HALF*FaultZ*Vfout
- return
- endif
- f_bisec_Vflow = TauStick - HALF*FaultZ*Vflow - cA*asinh(Vflow/(TWO*Vo)*exp(-help))
- if ( f_bisec*f_bisec_Vflow < 0._CUSTOM_REAL ) then
- Vfhigh = Vfbisec
- else
- Vflow = Vfbisec
- endif
- enddo
-
- ierror = 1
- !compute updated Vf before exiting NR search
- fa = tau/(Seff*cca)
- help = -(fo+ccb*psi)/cca
- help1 = exp(help+fa)
- help2 = exp(help-fa)
- Vfout = Vo*(help1 - help2)
- tauout = tau
- return
- endif
-
- END DO
-
- !compute the updated Vf before exiting NR search
- fa = tau/(Seff*cca)
- help = -(fo+ccb*psi)/cca
- help1 = exp(help+fa)
- help2 = exp(help-fa)
- Vfout = Vo*(help1 - help2)
- tauout = tau
-
-end subroutine NRsearch
-
!===============================================================
! OUTPUTS
subroutine init_dataT(DataT,coord,nglob,NT,iflt)
@@ -1139,6 +1140,7 @@
allocate(DataT%v2(NT,DataT%npoin))
allocate(DataT%t2(NT,DataT%npoin))
allocate(DataT%t3(NT,DataT%npoin))
+ allocate(DataT%theta(NT,DataT%npoin))
DataT%d1 = 0e0_CUSTOM_REAL
DataT%v1 = 0e0_CUSTOM_REAL
DataT%t1 = 0e0_CUSTOM_REAL
@@ -1146,16 +1148,18 @@
DataT%v2 = 0e0_CUSTOM_REAL
DataT%t2 = 0e0_CUSTOM_REAL
DataT%t3 = 0e0_CUSTOM_REAL
+ DataT%theta = 0e0_CUSTOM_REAL
close(IIN)
end subroutine init_dataT
!---------------------------------------------------------------
-subroutine store_dataT(dataT,d,v,t,itime)
+subroutine store_dataT(dataT,d,v,t,theta,itime)
type(dataT_type), intent(inout) :: dataT
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
integer, intent(in) :: itime
integer :: i,k
@@ -1169,6 +1173,7 @@
dataT%t1(itime,i) = t(1,k)
dataT%t2(itime,i) = t(2,k)
dataT%t3(itime,i) = t(3,k)
+ dataT%theta(itime,i) = theta(k)
enddo
end subroutine store_dataT
@@ -1210,12 +1215,12 @@
do i=1,dataT%npoin
open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "# problem=TPV16"
+ write(IOUT,*) "# problem=TPV102"
write(IOUT,*) "# author=Surendra Nadh Somala"
write(IOUT,1000) today(2), today(1), today(3), now
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=75 m (*5 GLL nodes)"
+ write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
write(IOUT,*) "# time_step=",DT
write(IOUT,*) "# location=",trim(dataT%name(i))
write(IOUT,*) "# Column #1 = Time (s)"
@@ -1226,15 +1231,26 @@
write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
write(IOUT,*) "# Column #8 = normal stress (MPa)"
+ if(Rate_AND_State) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
write(IOUT,*) "#"
write(IOUT,*) "# The line below lists the names of the data fields:"
- 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, &
+ if(.not. Rate_AND_State) then
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(8(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
-dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
dataT%t3(k,i)/1.0e6_CUSTOM_REAL
- enddo
+ enddo
+ else
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
+ enddo
+ endif
close(IOUT)
@@ -1265,12 +1281,12 @@
IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
open(IOUT,file=trim(filename),status='replace')
- write(IOUT,*) "# problem=TPV16"
+ write(IOUT,*) "# problem=TPV102"
write(IOUT,*) "# author=Surendra Nadh Somala"
write(IOUT,1000) today(2), today(1), today(3), now
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=75 m (*5 GLL nodes)"
+ write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
write(IOUT,*) "# Column #1 = horizontal coordinate, distance along strike (m)"
write(IOUT,*) "# Column #2 = vertical coordinate, distance down-dip (m)"
write(IOUT,*) "# Column #3 = rupture time (s)"
More information about the CIG-COMMITS
mailing list