c PROGRAM VC_STRESS_GREEN.F c c Computation of Stress Green's functions (influence functions) c for interaction of fault segments. Present code c allows possibility of both elastic and viscoelastic c Greens functions according to the metod of c Rundle (1988a) c c Modified July 2, 2002 c c ************************************************************* c c General Description: c c This is the code that computes the stress green's functions c according to the methods described in Rundle (1988a,b) for c a prescribed group of vertical strike slip fault segments c embedded in an elastic half space, or in an elastic layer c overlying an viscoelastic half space. The stress greens c function specifies the stress change on fault j due to c a unit of slip on fault i. Relative geometry of the fault c segments, as well as elastic or viscoelastic properties of c the medium, are important. The output file of this code is c used as input file to VC_SG_COMPRESS.F, which compresses c the viscoelastic part into a sum of 8 decaying exponentials c (collocation technique) for each fault segment's influence c functions. The output file of VC_SG_COMPRESS.F is then c used as input to the first Virtual California code, c VC_INIT_SER.f. The output of that code is then used as c input to VC_SER.f c c VC_STRESS_GREEN.f is therefore a Boundary Element Method c (BEM). c c References: c c The basic ideas behind these codes are described in the c following series of papers: c c Rundle, JB, A physical model for earthquakes, 1., Fluctuations c and interactions, J. Geophys. Res., 93, 6237 (1988a) c c Rundle, JB, A physical model for earthquakes, 2. Application c to Southern California, J. Geophys. Res., 93, 6255 c (1988b) c c Rundle, JB, Linear pattern dynamics in nonlinear threshold c systems, Phys. Rev. E, 61, 2418 (2001) c c Rundle, PB, JB Rundle, KF Tiampo, JSS Martins, S McGinnis, and c W Klein, Phys. Rev. Lett., 87, 148501 (2001) c c Rundle, JB, PB Rundle, W Klein, JSS Martins, KF Tiampo, c A Donnellan and LH Kellogg, GEM plate boundary c simulations for the plate boundary observatory: c Understanding the physics of earthquakes c on complex fault systems, Pure. Appl. Geophys., c in press (2002) c c c ************************************************************* c c Virtual California Code Physics: c c Virtual California uses topologically realistic networks of c independent fault segments that are mediated by elastic c interactions. (Note that earlier versions had viscoelastic c interactions as well, but such are not included in the present c code). VC is a "backslip" model, inasmuch as the plate c tectonic stress increases is produced by means of applying a c negative ("backslip") velocity to each segment whose magnitude c is that of the long-term rate of slip on the segment. Since c "positive slip" reduces the stress on a fault segment, c "negative slip" due to the backslip increases the stress. c On each time step, all faults are checked to determine whether c the shear stress has reached the failure threshold. Once c at least one segment reaches the threshold, the "long time c steps" stop, and "short (failure) time steps" (a.k.a. Monte c Carlo Sweeps, or mcs) begins. An mcs begins with a check of c of each site to determine whether it has failed, followed by c a parallel updating of each segment. An update of a segment c consists of increasing the sudden seismic slip on each segment c so that the stress of the segment, considered in isolation, c drops to a residual value, plus or minus a random overshoot/ c undershoot. The elastic stress on all segments is then c recalculated, and another mcs is carried out. This iterative c process repeats until all segments are below the failure c threshold, at which time the mcs time steps cease and the c long plate tectonic time steps begin again. Note that VC c also includes a stress-dependent "precursory slip", or c "stress leakage" of the type that has been observed in c laboratory experiments by TE Tullis (1996) and S Karner c and C Marone (2001). The physics of this process is that c as the stress on a segment increases, a small amount of stable c sliding occurs that is proportional to the level of the stress c above the residual. Lab experiments and field data suggest c that the parameter called "alpha" (in Rundle et al. 2001) c is of the order of a few percent. This parameter is called c facsdp() or facrat() or facsd() in this code. c c In addition, as described in Rundle (1988), the fault system c topology + elasticity + failure law may lead to an unstable, c "runaway" dynamics due to the presence of positive c eigenvalues, so we introduce a nonlinear, "self-adapting" or c "self-correcting" term. The VC initializer code estimates c the magnitude of these terms for each segment using a simple c but unrealistic dynamics. They are then readjusted by c running the VC code with real dynamics until all startup c transients (unrealistically large stresses) disappear. The c time scale for the transients to disappear depends on the c complexity of the problem. Physically, one can imagine c that the system evolves on a rough energy landscape, and c on average sits at the bottom of a free energy well. Large c earthquakes may tend to displace the system from the well c allowing non-stationary evolution to occur before the system c settles into a new well. The configuration of the well is c determined by the magnitude of the positive eigenvalues c discussed above. The eigenvalues for each segment are c proportional to the derivative of the shear stress with c respect to the slip on each fault segment. c c c ************************************************************* c c Program Modules: c c This code uses several modules and sub-programs, whose c variables will be described below. These modules are: c c Main code -- Reads in fault topology data, c rotates fault segment i into c a standard coordinate system c by an angle theta(), calls c subroutines to compute stress change c on fault segment j due to a unit c slip on fault segment i,then rotates c stresses back into the original c coordinate system. This routine c calls both CALT() and CALTRM() c for the actual stress computations. c c Subroutine CALT() -- Master subroutine that computes the c elastic stress influence coefficients. c For each fault segment i, calculates c the stress on all segments j at c a grid of 36 points uniformly c distributed over a rectangular grid c on fault segment j. Then averages c the stresses over the 36 points c to obtain an average stress on fault c segment j due to slip on fault segment c i. CALT() calls Subroutine STRELA() c to do the most complex parts of the c computation. c c Subroutine CALTRM() -- Master subroutine that computes the c viscoelastic stress influcence c coefficients due to the Maxwell c viscoelastic half space. For this c computation, we use the method of c images as discussed in Rundle (1988b). C The stress is computed on the target c fault segment j due to slip on the c master fault segment i. Averages c of stress are taken over 6 points on c the target fault. c c Subroutine STRELA() -- Computes the elastic strains on the c target fault segment j due to unit c of slip on master fault segment i. c Method used is that of Chinnery (BSSA, c 51, 355, 1961 as corrected), or c Maruyama (Bull. EQ Res. Inst. Tokyo, c 42, 289, 1964), or Mansinha & Smylie c (BSSA, 61, 1433, 1971). c c Subroutine STRVIS() -- Computes the strains at location x_j c due to slip at a point source at c location x_i. These strains are used c to compute the viscoelastic strains c due to mend image sources. Point c sources are used instead of rectangular c sources since only long wavelength, c or far-field effects are important c when computing the viscoelastic effects. c c Function sgn() -- Signum function (x < 0, or x > 0) c c Function jsecond() -- Computes cpu time for purposes of c program timing. c c ************************************************************* c c Dimensions from INCLUDE statement: c c NF -- MAX number of faults -- should be same as in EQPARM.FOR c file c NTSTP -- Max number of time steps c NFT -- Same as NTSTP c NF1 -- Same as NF c NF2 -- Same as NF c NFF -- NF * NF c NMF -- > NF * NF * 8 c NPT -- > NF * 36 c c c ************************************************************* c c Major Program Variables: c c 1. INPUT FILE. The primary input is a file with the c coordinates of the fault segments at which the c stress Green's functions are to be computed. All c fault segments are vertical strike slip faults. c Note that the coordinate system is a East = x, North c = y rectangular coordinate system. Units are in km c from the point 31 deg N Lat, 121 deg W Long. c All computations are carried out in double precision c in FORTRAN 77. The input variables are: c c numflt - Number of fault segments c mend - Number of viscoelastic images c vplatx - x-coordinate of far field plate tectonic c velocity in cm/yr c vplaty - y-coordinate of far field plate tectonic c velocity in cm/yr c c dbn(NF) - Depth in km to bottom of fault segment c dtn(NF) - Depth in km to top of fault segment c xfwn(NF)- x-coordinate of west end of fault segment c yfwn(NF)- y-coordinate of west end of fault segment c xfen(NF)- x-coordinate of east end of fault segment c yfen(NF)- y-coordinate of east end of fault segment c slpvel(NF) - Slip velocity in cm/yr. Convention is c "+" (>0) for right lateral strike slip, "-" c (<0) for left lateral strike slip. c c ************************************************************* c c 2. OUTPUT FILE. The primary output (other than program c diagnostic printouts to screen) is a file that c contains the shear stress and normal stress Green's c functions, or influence coefficients. Both elastic c and viscoelastic Green's functions are stored. c Convention is "+" is right lateral shear stress, c "-" is left lateral shear stress. c c The primary outputs variables are: c c fname5 - Character string name of input file c numflt - Number of fault segments c mend - Number of viscoelastic images c hplate - Elastic plate thickness in km c vplatx - x-coordinate of far field plate tectonic c velocity in cm/yr c vplaty - y-coordinate of far field plate tectonic c velocity in cm/yr c bulkm - Elastic bulk modulus in bars c rhofcc(NF) - Pressure at center of fault segment in c bars c c Ssh(NFF)- Shear stress elastic Green's function for c influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j. c c Snorm(NFF)- Normal stress elastic Green's function for c influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j. c c Sshm(NMF)- Shear stress viscoelastic Green's function c for influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N*mend coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j, all c mend images. c c Snormm(NMF)- Normal stress viscoelastic Green function c for influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N*mend coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j, all c mend images. c c area(NF) - areas of segments in km^2 c c svec(NF) - represents the total shear stress change c per year, in inverse yr, (yr^-1) c (i.e., normalized by amu in bars), due to c the loading velocity on all the faults. c c snvec(NF) - Represents the total normal stress change c per year, in inverse yr, (yr^-1) c (i.e., normalized by amu in bars), due to c the loading velocity on all the faults. c c slpvel(NF) - Slip velocity in cm/yr. Convention is c "+" (>0) for right lateral strike slip, "-" c (<0) for left lateral strike slip. c c vcovec(NFF) - Used in computing svec(NF). Represents c the total influence of a unit of slip on c fault segment j on the shear stress change c on fault segment i. c c vnvec(NFF) - Used in computing snvec(NF). Represents c the total influence of a unit of slip on c fault segment j on the normal stress change c on fault segment i. c c ************************************************************* c ************************************************************* c implicit real*8(a-h,o-z) c integer recrd character fm7*23,fm8*15,fm9*11 character fma80*7 character resp4*1 character fmt1*8,fmt2*15,fmt3*6,fmt4*14,fma*7 character fname*20, fname5*20 real*4 tsec_1, tsec_2, tsec_3, tsec_4, tdfsec, jsecond INCLUDE 'EQPARM_4.FOR' dimension x(NPT),y(NPT),z(NPT),alen(NF) dimension x1(NPT), y1(NPT), z1(NPT) dimension dbn(NF),dtn(NF),xfwn(NF),yfwn(NF),slpvel(NF), & xfen(NF),yfen(NF),unitx(NF),unity(NF) dimension xxfe(NF), yyfe(NF), xxfw(NF), yyfw(NF) dimension thetan(NF),xmidn(NF),ymidn(NF),area(NF) dimension Sxy(NF),Sn(NF),Sxym(20,NF),Snm(20,NF) dimension unx(NF),uny(NF),vcovec(NFF) dimension vcof(NF,NF),svec(NF), vnvec(NFF) dimension vncof(NF,NF),snvec(NF) dimension Ssh(NFF), Snorm(NFF), Sshm(NMF), Snormm(NMF) dimension rhofcc(NF) external sgn, jsecond c c ************************************************************* c c Open file faultdata c c ************************************************************* c data fmt1,fmt2,fmt3/'(2i10)','(7g20.10,i10)','(a1)'/ data fmt4/'(2g20.10,a1)'/ data fma80,fma /'(a80)','(a20)'/ data fm8,fm9 /'(a80,2f20.10)','(6f20.10)'/ data fm7 /'(2f20.10,i10,2f20.10)'/ c c Get the data file name c print *,' What data file name?' read(*,fma) fname open(unit=20,name=fname,status='old', .access='direct',form='formatted',recl=160) fname5=fname c c Open the output file c print *,' What output file name?' read(*,fma) fname open(unit=30,name=fname,status='new') c c Query for printing stress coefficients c resp4 = 'n' c c Define pi c pi=3.1415926535d0 c c Get the number of faults and number of terms in the Green function c expansion c read(20,fmt=*) numflt,mend c c print *, ' What value for Elastic Plate Thickness (km)? ' c read *, hplate c print *, ' ' print *, ' Viscoelastic physics is turned off in this code' print *, ' for now.' c c However, we must assume some value for elastic plate thickness c for code computations. Hence we take hplate = 30 km c even though it is not used in the final simulator. c hplate = 30. c c ...since only elastic deformations are allowed for the moment c bulkm = 1.66666 alam=bulkm - 2.d0/3.d0 c c Set up loop over the faults c do i=1,numflt recrd=i+1 read(20,fmt=*) dbn(i),dtn(i),xfwn(i),yfwn(i), . xfen(i),yfen(i),slpvel(i) c c Change sign of y so that +=south c yfen(i)=-yfen(i) yfwn(i)=-yfwn(i) rhofcc(i)=2.646d0*(dbn(i) + dtn(i))/2.d0 c end do c c Get the plate velocities c recrd=numflt + 2 read(20,fmt=*) vplatx,vplaty c c Change sign of vplate y to coincide with south facing y-axis c vplaty=-vplaty c c Indexes for Ssh and Sshm c indn=0 indnm=0 do i=1,numflt c c Define variables c xfe=xfen(i) yfe=yfen(i) xfw=xfwn(i) yfw=yfwn(i) dt=dtn(i) db=dbn(i) alen(i)=dsqrt((xfe-xfw)**2 + (yfe-yfw)**2) area(i)=alen(i)*(db-dt) c c Compute the fault midpoints c xmidn(i)=.5d0*(xfw+xfe) ymidn(i)=.5d0*(yfw+yfe) c c Compute the rotation angle c c Test to see if fault is E-W c testy=yfw-yfe if(dabs(testy) .le. 1.d-10) then xmax=dmax1(xfw,xfe) xmin=dmin1(xfw,xfe) xfw=xmax xfe=xmin thetan(i)=0.d0 end if c delx=xfe-xmidn(i) dely=yfe-ymidn(i) if(delx .le. 1.d-10) delx=(1.d-10)*sgn(delx) thetan(i)=datan(dely/delx) c c Compute fault normal unit vector c theta=thetan(i) unitx(i)=-dsin(theta) unity(i)=dcos(theta) end do npoint=numflt * 36 c c Rotate observation points and compute stresses c write(30,'(a20)') fname5 write(30,*) numflt,mend,hplate,vplatx,vplaty, . bulkm c c ************************************************************* c c Begin timing the code c c Begin timing the fortran program c tsec_1 = jsecond() print *, ' ' print *, ' tsec_1: ',tsec_1 c c ************************************************* c c In the following loop, i is the column index, c and generally, j is the row index c c ************************************************* c do i=1,numflt st=dsin(thetan(i)) ct=dcos(thetan(i)) xmid=xmidn(i) ymid=ymidn(i) c do j=1,numflt unx(j)= unitx(j)*ct + unity(j)*st uny(j)=-unitx(j)*st + unity(j)*ct xxfe(j)= (xfen(j)-xmid)*ct + (yfen(j)-ymid)*st yyfe(j)=-(xfen(j)-xmid)*st + (yfen(j)-ymid)*ct xxfw(j)= (xfwn(j)-xmid)*ct + (yfwn(j)-ymid)*st yyfw(j)=-(xfwn(j)-xmid)*st + (yfwn(j)-ymid)*ct end do c c Calculate the observation points for elastic c sources c n=0 do nr=1,numflt zstart=(dbn(nr)-dtn(nr))/12.d0 + dtn(nr) delzs=(dbn(nr)-dtn(nr))/6.d0 do nj=1,6 xstart=(xxfe(nr)-xxfw(nr))/12.d0 + xxfw(nr) ystart=(yyfe(nr)-yyfw(nr))/12.d0 + yyfw(nr) delxs=(xxfe(nr)-xxfw(nr))/6.d0 delys=(yyfe(nr)-yyfw(nr))/6.d0 do ni=1,6 n=n + 1 x1(n)= xstart y1(n)= ystart z1(n)= zstart ystart=ystart + delys xstart=xstart + delxs end do zstart=zstart + delzs end do end do c ns=(i-1)*36 do ni=1,36 nii=ns + ni y1(nii)=y1(nii) + .0001d0 end do c c Calculate the observation points for viscoelastic c m-sources c n=0 do nr=1,numflt zstart=(dbn(nr)-dtn(nr))*.25d0 + dtn(nr) delzs=(dbn(nr)-dtn(nr))*.5d0 do nj=1,2 xstart=(xxfe(nr)-xxfw(nr))/6.d0 + xxfw(nr) ystart=(yyfe(nr)-yyfw(nr))/6.d0 + yyfw(nr) delxs=(xxfe(nr)-xxfw(nr))/3.d0 delys=(yyfe(nr)-yyfw(nr))/3.d0 do ni=1,3 n=n + 1 x(n)= xstart y(n)= ystart z(n)= zstart ystart=ystart + delys xstart=xstart + delxs end do zstart=zstart + delzs end do end do c ns=(i-1)*6 do ni=1,6 nii=ns + ni y(nii)=y(nii) + .0001d0 end do c sl=.5d0*dsqrt((xfwn(i)-xfen(i))**2+(yfwn(i)-yfen(i))**2) xfe=xfen(i) yfe=yfen(i) xfw=xfwn(i) yfw=yfwn(i) dt=dtn(i) db=dbn(i) c c Calculate displacements c amu=1.d0 c c Change sign so that right lateral slip is positive c uds=-1.d0 c if(i .eq. 1) then print *,' ' print *,' ******************************************' print *,' ' print *,' REMINDER: Right lateral slip is positive,' print *,' and right lateral shear stress on the fault' print *,' is positive!' print *,' ' print *,' ******************************************' print *, ' ' end if c xc = 0.d0 yc = 0.d0 zc = ( dtn(i) + dbn(i) ) / 2.d0 call CALT(Sxy,Sn,x1,y1,z1,xc,yc,zc,uds,sl,db,dt, & amu,alam,npoint,numflt,i,unx,uny,dbn,dtn,xxfe,xxfw,yyfe, & yyfw, area) c call CALTRM(Sxym,Snm,x,y,z,hplate,uds,area, & sl,db,dt,amu,alam,mend,numflt,i,unx,uny) c do n=1,numflt indn = indn + 1 Ssh(indn)=Sxy(n) Snorm(indn)=Sn(n) end do c do n=1,numflt do m=1,mend indnm=indnm + 1 Sshm(indnm)=Sxym(m,n) Snormm(indnm)=Snm(m,n) end do end do c do j=1,numflt sum=0.d0 sumn=0.d0 do m=1,mend sum=sum + Sxym(m,j) sumn=sumn + Snm(m,j) end do vcof(i,j)=sum + Sxy(j) vncof(i,j)=sumn + Sn(j) end do c c Print out the coefficients c c This is the end do that goes with the loop over i c print *, ' ' print *, 'JUST FINISHED FAULT: ', i end do c c Finish timing the fortran code c tsec_2 = jsecond() tdfsec = tsec_2 - tsec_1 print *, ' ' print *, ' CPU Time for main part of code: ', tdfsec c c Symmetrize the shear stress matrices c do ir=2,numflt do ic=1,ir-1 vcl=vcof(ir,ic)*area(ic) vcu=vcof(ic,ir)*area(ir) vcof(ir,ic)=.5d0*(vcu + vcl)/area(ic) vcof(ic,ir)=.5d0*(vcu + vcl)/area(ir) indl=(ic-1)*numflt + ir indu=(ir-1)*numflt + ic sxyl=Ssh(indl)*area(ir) sxyu=Ssh(indu)*area(ic) Ssh(indl)=.5d0*(sxyl + sxyu)/area(ir) Ssh(indu)=.5d0*(sxyl + sxyu)/area(ic) indlm=(indl-1)*mend indum=(indu-1)*mend do m=1,mend indlmm=indlm + m indumm=indum + m shml=Sshm(indlmm)*area(ir) shum=Sshm(indumm)*area(ic) Sshm(indlmm)=.5d0*(shml + shum)/area(ir) Sshm(indumm)=.5d0*(shml + shum)/area(ic) end do end do end do c c Write the output into file 30 c indmx=numflt*numflt indnmx=numflt*numflt*mend write(30,*) (rhofcc(n), n=1,numflt) write(30,*) (Ssh(n), n=1,indmx) write(30,*) (Snorm(n), n=1,indmx) write(30,*) (Sshm(n), n=1,indnmx) write(30,*) (Snormm(n), n=1,indnmx) c vmag=dsqrt(vplatx**2 + vplaty**2) vpx=vplatx/vmag vpy=vplaty/vmag c c Write the areas out c write(30,*) (area(n), n=1,numflt) c do n=1,numflt svec(n)=0.d0 snvec(n)=0.d0 do i=1,numflt svec(n)=svec(n) - slpvel(i)*vcof(i,n) snvec(n)=snvec(n) - slpvel(i)*vncof(i,n) end do end do c indvx=0 do n=1,numflt do j=1,numflt indvx=indvx + 1 vcovec(indvx)=vcof(n,j) vnvec(indvx)=vncof(n,j) end do end do write(30,*) (svec(n), n=1,numflt) write(30,*) (snvec(n), n=1,numflt) write(30,*) (slpvel(n), n=1,numflt) write(30,*) (vcovec(n), n=1,indvx) write(30,*) (vnvec(n), n=1,indvx) c close(unit=30,disp='keep') close(unit=20,disp='keep') stop end c subroutine CALT(Sxy,Sn,xpt,ypt,zpt,xc,yc,zc,uds, & sl,db,dt,amu,alam,npoint,nfault,ifault,unitx,unity, & dbn,dtn,xxfe,xxfw,yyfe,yyfw,area) c c ************************************************************* c c INTERFACE: c Sxy(i)--shear stress change on the i-th fault due to c slip on the ifault-th fault c Sn(i)--normal stress on the i-th fault due to slip on the c ifault-th fault c xpt(36), ypt(36), zpt(36)-- Center points of all 36 c sub-segments at which we are computing strains c xc, yc, zc -- Center points of ifault-th fault c uds -- Slip on ifault-th fault (= 1 unit of right lateral c strike slip, or -1.0) c sl--the fault semilength (km) c db -- depth to bottom of ifault-th fault in km c dt -- depth to top of ifault-th fault in km c amu -- shear modulus (but set to 1.0 here) c alam -- second lame parameter (set to 1.0 also) c npoint -- total number of subsegments, = nfault*36 c nfault--the number of faults (up to NF) c ifault--the i-th fault which we are presently on and which is c slipping by a unit amount c unitx--an array which contains the x-component (along strike) c of the j-th fault in the coordinate system of the i-th c fault (the i-th fault is the one considered here) c unity--an array which contains the y-component (perp. c to strike) of the j-th fault in the coordinate c system of the i-th fault c dbn(i) -- depth to bottom of i-th fault segment c dtn(i) -- depth to bottom of i-th fault segment c xxfe(i) -- x-coordinate of east end of i-th fault segment c xxfw(i) -- x-coordinate of west end of i-th fault segment c yyfe(i) -- y-coordinate of east end of i-th fault segment c yyfw(i) -- y-coordinate of west end of i-th fault segment c area(i)--areas of all the fault patches c c RETURN VALUE: c Sxy, Sn (see above) c c DESCRIPTION: c Master subroutine that computes the elastic stress influence c coefficients. c For each fault segment i, calculates the stress on all segments c j at a grid of 36 points uniformly distributed over a rectangular c grid on fault segment j. Then averages the stresses over the c 36 points to obtain an average stress on fault segment j due to c slip on fault segment i. c CALT() calls Subroutine STRELA() to do the most complex parts of c the computation. c c AUXILIARY VARIABLES: c exx(36)--normal strain along strike at each of 36 c sub-segments c eyy(36)--normal strain perpendicular to strike at each c of 36 sub-segments c ezz(36)--normal strain in vertical direction at each of c 36 sub-segments c exy(36)--shear strain along fault at each of 36 sub- c segments c c ************************************************************* c implicit real*8(a-h,o-z) INCLUDE 'EQPARM_4.FOR' c dimension xpt(NPT),ypt(NPT),zpt(NPT) dimension xj(36), yj(36), zj(36) dimension xr(36), yr(36), zr(36) dimension exx(36),eyy(36),exy(36),ezz(36) dimension exxsub(36),eyysub(36) dimension exysub(36),ezzsub(36) dimension unitx(NF),unity(NF) dimension Sxy(NF),Sn(NF) dimension dbn(NF),dtn(NF) dimension area(NF) dimension xxfe(NF),xxfw(NF),yyfe(NF),yyfw(NF) c c ******************************************************* c c Zero out stress variables c do i=1,nfault Sxy(i)=0.d0 Sn(i)=0.d0 end do c npoint=nfault*36 xi=xc yi=yc zi=zc c c We are on fault *ifault*, and we want to compute its c influence on all the other faults c c c *************** c c Loop over faults c do jfault = 1,nfault jps=(jfault-1)*36 do j=1,36 xj(j)=xpt(jps + j) yj(j)=ypt(jps + j) zj(j)=zpt(jps + j) end do do k=1,36 exx(k)=0.d0 eyy(k)=0.d0 ezz(k)=0.d0 exy(k)=0.d0 end do c c Now we are going to loop over all of the target points c on fault = jfault. c c The fault we are on is ifault c c ktargt is loop for target subfault c do ktargt = 1,36 xr(ktargt) = xj(ktargt) - xi yr(ktargt) = yj(ktargt) - yi & - .0001d0 zr(ktargt) = zj(ktargt) end do do k=1,36 exxsub(k)=0.d0 eyysub(k)=0.d0 ezzsub(k)=0.d0 exysub(k)=0.d0 end do c c Define fault parameters for the ksourc-th subfault c uds = -1.d0 nsbseg = 36 call STRELA(db,dt,sl,xr,yr,zr,uds, & exxsub,eyysub,ezzsub,exysub, & amu,alam,nsbseg) do ktargt = 1,36 exx(ktargt) = exx(ktargt) + exxsub(ktargt) eyy(ktargt) = eyy(ktargt) + eyysub(ktargt) ezz(ktargt) = ezz(ktargt) + ezzsub(ktargt) exy(ktargt) = exy(ktargt) + exysub(ktargt) end do c c Find average of all 36 points c sumx = 0.d0 sumy = 0.d0 sumz = 0.d0 sumxy = 0.d0 do k = 1,36 sumx = sumx + exx(k) sumy = sumy + eyy(k) sumz = sumz + ezz(k) sumxy = sumxy + exy(k) end do sumx = sumx / 36.d0 sumy = sumy / 36.d0 sumz = sumz / 36.d0 sumxy = sumxy / 36.d0 delt = sumx + sumy + sumz c uxsq=unitx(jfault)*unitx(jfault) uysq=unity(jfault)*unity(jfault) uxuy=unitx(jfault)*unity(jfault) c Sxy(jfault)=2.d0*amu*(uxuy*(sumx-sumy) + & sumxy*(uysq - uxsq)) Sxy(jfault)=-Sxy(jfault) Sn(jfault)=(alam*delt + 2.d0*amu*(uxsq*sumx + uysq*sumy & + 2.d0*uxuy*sumxy)) c c End loop over faults c end do c return end subroutine STRELA(db,dt,sl,xpt,ypt,zpt,uds,exx,eyy, . ezz,exy,amu,alam,nsbseg) c c *************************************************** c c INTERFACE: c db--the depth to the fault bottom (km) c dt--the depth to the fault top (km) c sl--the fault semilength (km) c xpt(i),ypt(i),zpt(i)--coordinates of the observation c point in chinnery's coordinates (km) c uds--fault slip (m) c exx(i)--normal strain in direction parallel to fault c eyy(i)--normal strain in direction perpendicular to fault c ezz(i)--normal strain in vertical direction c exy(i)--shear strain in horizontal plane c amu -- shear modulus (but set to 1.0 here) c alam -- second lame parameter (set to 1.0 also) c nsbseg--the number of target fault sub-segments c c RETURN VALUE: c exx, eyy, ezz, exy (see above) c c DESCRIPTION: c This subroutine evaluates the strains arising from c slip on a vertical strike slip fault in c a homogeneous, elastic half space. Uses the methods c as described in Rundle (1988a). c c ****************************************************** c implicit real*8(a-h,o-z) INCLUDE 'EQPARM_4.FOR' dimension exx(36),eyy(36),exy(36),ezz(36), & xpt(36),ypt(36),zpt(36) c c ******************************************************** c c Note that we have to multiply x,u1,and u by -1 to be c consistent with a coordinate system in which the x axis is c reversed from Chinnery, and right lateral slip is positive c fourpi=12.5663706140d0 u=(uds)/fourpi al=sl ddb=db ddt=dt dalam=alam damu=amu alph=(dalam+damu)/(dalam+2.d0*damu) beta=(dalam-damu)/(dalam+damu) cee=damu*dalam/((dalam+damu)**2) alph1=1.d0-alph beta1=1.d0-beta cee1=1.d0+cee c do jk=1,nsbseg x=xpt(jk) y=ypt(jk) z=zpt(jk) c exx(jk)=0.0d0 eyy(jk)=0.0d0 ezz(jk)=0.0d0 exy(jk)=0.0d0 c c Define parameters c du1dx=0.d0 du1dy=0.d0 du2dx=0.d0 du2dy=0.d0 du3dz=0.d0 do j=1,4 if(j .eq. 1) then x1=al x3=ddb sign=1.d0 else if(j .eq. 2) then x1=-al x3=ddb sign=-1.d0 else if(j .eq. 3) then x1=al x3=ddt sign=-1.d0 else if(j .eq. 4) then x1=-al x3=ddt sign=1.d0 end if c p=x3+z q=x3-z t=x1-x s1=dsqrt(t*t+y*y+q*q) s2=dsqrt(t*t+y*y+p*p) s1sq=t*t+y*y+q*q s2sq=t*t+y*y+p*p s1q=s1+q s2p=s2+p s2psq=s2p*s2p s1qsq=s1q*s1q s2pl=s2+s2p s1ql=s1+s1q pqmsq=p*p-q*q s1m=1.d0/s1 s2m=1.d0/s2 c c ************************************************************* c c Define quantities and then calculate strain c c ************************************************************* c vxy=1.d0/(s1*(s1q)) . + ((cee1)*s2+(beta1)*p+q)/(s2*s2psq) . - ((p*p-q*q)*(s2pl))/(2.d0*s2*s2sq*s2psq) c dend4=s2*s2psq d4x=t/dend4 d4y=-y/dend4 c dend5=s1*s1sq*s1qsq anum5=s1+s1q d5=anum5/dend5 d5x=t*d5 d5y=-y*d5 c dend6=s2*s2sq*s2psq anum6=s2pl d6=anum6/dend6 d6x=t*d6 d6y=-y*d6 c d7=(s2+s2pl)/(s2*s2sq*s2p*s2psq) d7x=t*d7 d7y=-y*d7 c d8=(2.d0*s2+3.d0*s2p)/(s2sq*s2sq*s2*s2psq) d8x=t*d8 d8y=-y*d8 c dat1=y*y*s1sq+q*q*t*t dat1=dat1*s1 dat2=y*y*s2sq+p*p*t*t dat2=dat2*s2 ddat1x=-y*q*(y*y+q*q)/dat1 ddat2x=-y*p*(y*y+p*p)/dat2 ddat1y=-q*t*(y*y+s1sq)/dat1 ddat2y=-p*t*(y*y+s2sq)/dat2 c den1=s2sq*s2psq den2=s2sq*den1 term1=cee1/den1 term2=cee1*s2 + beta1*p + q term3=pqmsq/den2 term4=.5d0*pqmsq*s2pl term5=alph1/(s1*s1q) term6=(alph1-alph*cee)/(s2*s2p) term7=alph*((beta+cee)*p - q) term8=.5d0*alph*pqmsq term9=d5x - t*term1 . + term2*d7x + t*term3 . -term4*d8x term0=d5y + y*term1 + term2*d7y - y*term3 . - term4*d8y c ds1dz=-q*s1m ds2dz=p*s2m ds1mdz=-.5d0*ds1dz*s1m*s1m*s1m ds2mdz=-.5d0*ds2dz*s2m*s2m*s2m d1z=s2m/s2p dd1zdz=-(2.d0*s2*ds2dz + s2 + p*ds2dz)*(d1z*d1z) sm3=s2m*s2m*s2m dsm3dz=-3.d0*ds2dz*sm3*s2m*s2m c du1dx=du1dx + sign*( -y*alph*vxy . + y*t*alph*(term9) + ddat1x + ddat2x) c du1dy=du1dy + sign*( t*alph*vxy + y*t*alph* . (term0) + ddat1y + ddat2y) c du2dx=du2dx + sign*( -t*(term5 + term6) . - term7*d4x - term8*d6x -y*y*alph*term9 ) c du2dy=du2dy + sign*( y*(term5 + term6) . - term7*d4y - term8*d6y -2.d0*y*alph*vxy . - y*y*alph*term0) c du3dz=du3dz + sign*(y*(alph*(ds1mdz-ds2mdz) + . d1z*(ds2dz*beta1 + 2.d0*alph - 2.d0*alph1) + . dd1zdz*(s2*beta1 + 2.d0*alph*p + 2.d0*q*alph1) + . alph*(p+q)*sm3 + .5d0*alph*pqmsq*dsm3dz )) end do c exx(jk)=du1dx*u eyy(jk)=du2dy*u ezz(jk)=du3dz*u exy(jk)=.5d0*(du1dy + du2dx)*u c end do return end c c c subroutine CALTRM(Sxym,Snm,x,y,z,h,uds,area, & sl,db,dt,amu,alam,mend,nfault,ifault,unitx,unity) c c ************************************************************* c c INTERFACE: c Sxym -- m-dependent terms for viscoelastic contribution c to shear stress Green's function. Units of inverse c slip. c Snm -- m-dependent terms for viscoelastic contribution c to normal stress Green's function. Units of inverse c slip c x(i),y(i),z(i) -- coordinates of points on target faults c at which viscoelastic stress is to be computed (km). c Here we use only 6 points on each target fault segment c rather than the 36 used for the elastic stress, since c viscoelastic contribution is long wavelength and c does not vary strongly in space. c h--the layer thickness (km) c uds -- unit of slip on fault, set to -1.0 assuming right c lateral strike slip c area(n) -- area of c sl -- fault segment semi-length of source (ifault-th) fault c db -- depth to bottom of source (ifault-th) fault c dt -- depth to top of source (ifault-th) fault c amu -- shear modulus (but set to 1.0 here) c alam -- second lame parameter (set to 1.0 also) c mend--the total number of terms kept in the c summation (number of images) c nfault--the number of faults (up to NF) c ifault--the i-th fault which we are presently on c unitx--an array which contains the x-component (along strike) c of the j-th fault in the coordinate system of the i-th c fault (the i-th fault is the one considered here) c unity--an array which contains the y-component (perp. to strike) c of the j-th fault in the coordinate system of the i-th c fault c c RETURN VALUE: c Sxym, Snm (see above) c c DESCRIPTION: c This subroutine calculates the m-dependent c terms for the Rundle-Jackson green function c See JB Rundle and DD Jackson (Geophys. J. Roy. Astr. c Soc., 49, 575, 1977). c c AUXILIARY VARIABLES: c exx(i)--normal strain in direction along strike c eyy(i)--normal strain in direction perpendicular to strike c exy(i)--shear strain c c ************************************************** c implicit real*8(a-h,o-z) INCLUDE 'EQPARM_4.FOR' dimension xpt(NPT),ypt(NPT),zpt(NPT),zmh(NPT),zph(NPT) dimension xs(6),ys(6),zs(6),x(NPT),y(NPT),z(NPT) dimension exxm(NF,20),eyym(NF,20),exym(NF,20) dimension exmm(NPT),eymm(NPT),eshmm(NPT) dimension expp(NPT),eypp(NPT),eshpp(NPT) dimension unitx(NF),unity(NF),area(NF) dimension Sxym(20,NF),Snm(20,NF) c c ******************************************************* c npntv=nfault*6 ips=(ifault-1)*6 do i=1,6 xs(i)=x(ips + i) ys(i)=y(ips + i) - .05d0*h zs(i)=z(ips + i) end do c if(npntv .gt. NPT) pause 'number of observation . points too big' c do m=1,mend do i=1,nfault exxm(i,m)=0.d0 eyym(i,m)=0.d0 exym(i,m)=0.d0 end do end do c do np=1,6 depsor=zs(np) c do i=1,npntv xpt(i)=x(i) - xs(np) ypt(i)=y(i) - ys(np) zpt(i)=z(i) end do c do m=1,mend do ipoint=1,npntv zmh(ipoint)=2.d0*dfloat(m)*h - zpt(ipoint) zph(ipoint)=2.d0*dfloat(m)*h + zpt(ipoint) end do c c Call the half space fault calculations for varying m c call STRVIS(xpt,ypt,zmh,uds,depsor,exmm,eymm,eshmm, & npntv) c call STRVIS(xpt,ypt,zph,uds,depsor,expp,eypp,eshpp, & npntv) c c Write strains into strain matrices c do i=1,nfault do ni=1,6 k=(i-1)*6 + ni exxm(i,m)=exxm(i,m) + (exmm(k) + expp(k))*.5d0 eyym(i,m)=eyym(i,m) + (eymm(k) + eypp(k))*.5d0 exym(i,m)=exym(i,m) + (eshmm(k) + eshpp(k))*.5d0 end do end do c end do c c c Now do case with faults getting deeper with time c do m=1,mend dm=2.d0*dfloat(m)*h - depsor dp=2.d0*dfloat(m)*h + depsor c c Call the half space fault calculations for varying m c call STRVIS(xpt,ypt,zpt,uds,dm,exmm,eymm,eshmm,npntv) c call STRVIS(xpt,ypt,zpt,uds,dp,expp,eypp,eshpp,npntv) c c Write strains into strain matrices c do i=1,nfault do ni=1,6 k=(i-1)*6 + ni exxm(i,m)=exxm(i,m) + (exmm(k) + expp(k))*.5d0 eyym(i,m)=eyym(i,m) + (eymm(k) + eypp(k))*.5d0 exym(i,m)=exym(i,m) + (eshmm(k) + eshpp(k))*.5d0 end do end do c end do end do c do i=1,nfault c uxsq=unitx(i)*unitx(i) uysq=unity(i)*unity(i) uxuy=unitx(i)*unity(i) c c IMPORTANT: Right lateral shear stress is positive, c so we need to multiply Sxy and Sxym by -1, so that right c lateral slip brings about a decrease in right lateral c shear stress on the fault. c c Do m-th order terms c do m=1,mend sumxm= exxm(i,m) sumym= eyym(i,m) sumxym= exym(i,m) c c Define stress for m-th order term c deltm=sumxm + sumym Sxym(m,i)=2.d0*amu*(uxuy*(sumxm-sumym) + . sumxym*(uysq-uxsq)) Sxym(m,i)=-Sxym(m,i)*area(ifault) Snm(m,i)=(alam*deltm + 2.d0*amu*(uxsq*sumxm + uysq*sumym . + 2.d0*uxuy*sumxym)) Snm(m,i)=Snm(m,i)*area(ifault) c end do end do c do m=1,mend do i=1,nfault Sxym(m,i)=Sxym(m,i)/36.d0 Snm(m,i)=Snm(m,i)/36.d0 end do end do c return end c subroutine STRVIS(xpt,ypt,zpt,uds,depsor,exx, . eyy,exy,npntv) c c ************************************************************* c c INTERFACE: c xpt(i),ypt(i),zpt(i)--coordinates of the observation c point in chinnery's coordinates (km) c uds--fault slip, set to -1.0 here for right lateral c strike slip c depsor--depth of the source below the surface c exx(i)--normal strain in direction parallel to fault c eyy(i)--normal strain in direction perpendicular to fault c exy(i)--shear strain in horizontal plane c npntv--the number of observation points at which strain c is computed c c RETURN VALUE: c exx, eyy, exy (see above) c c DESCRIPTION: c This subroutine calculates the strains at the c observation points due to an viscoelastic image point source c at depth depsor at the origin c c ****************************************************** c implicit real*8(a-h,o-z) INCLUDE 'EQPARM_4.FOR' dimension exx(NPT),eyy(NPT),exy(NPT), . xpt(NPT),ypt(NPT),zpt(NPT) c c ******************************************************** c c Note that we have to multiply x,u1,and u by -1 to be c consistent with a coordinate system in which the x axis is c reversed from Chinnery, and right lateral slip is positive c atepi=25.13274123d0 u=(uds)/atepi c c=depsor do jk=1,npntv x=xpt(jk) y=ypt(jk) z=zpt(jk) c exx(jk)=0.0d0 eyy(jk)=0.0d0 exy(jk)=0.0d0 c c Define parameters c p=c+z s1=dsqrt(x*x+y*y+(z+c)*(z+c)) s2=x*x + y*y + (z+c)*(z+c) s3=s2*s1 s4=s2*s2 s5=s4*s1 s6=s2*s2*s2 s7=s6*s1 s8=s4*s4 c r1=dsqrt(x*x + y*y + (z-c)*(z-c)) r2=x*x + y*y + (z-c)*(z-c) r3=r1*r2 r4=r2*r2 r5=r4*r1 r6=r2*r4 c rr1=dsqrt(x*x + y*y) rr2=x*x + y*y rr3=rr2*rr1 rr4=rr2*rr2 rr5=rr4*rr1 c drdx=x/r1 drdy=y/r1 dsdx=x/s1 dsdy=y/s1 c drrdx=x/rr1 drrdy=y/rr1 c dxrdx=2.d0*x/rr2 - 2.d0*x*x*drrdx/rr3 dxrdy=-2.d0*x*x*drrdy/rr3 dyrdx=-2.d0*y*y*drrdx/rr3 dyrdy=2.d0*y/rr2 - 2.d0*y*y*drrdy/rr3 c c ************************************************************* c c Define quantities and then calculate strain c c ************************************************************* c A1=(4.d0*((p-2.d0*c)**2)/3.d0) A2=(2.d0*((p-2.d0*c)**4)/3.d0) c A=2.d0*r1/3.d0 - A1/r1 + A2/r3 c Ap1=(10.d0*p*p/3.d0 + 8.d0*p*c - 8.d0*c*c) Ap2=(16.d0*p*p*c*c - 2.d0*p*p*p*p/3.d0 - 16.d0*p*p*p*c) Ap3=(8.d0*p*p*p*p*p*c - 8.d0*p*p*p*p*c*c) c Ap=-4.d0*p + 4.d0*s1/3.d0 + Ap1/s1 + Ap2/s3 . + Ap3/s5 c B1=12.d0*(p-2.d0*c)*(p-2.d0*c) B2=12.d0*((p-2.d0*c)**4) B3=4.d0*((p-2.d0*c)**6) c B=4.d0*r1 - B1/r1 + B2/r3 - B3/r5 c Bp1=40.d0*c*c - 36.d0*p*p - 40.d0*p*c Bp2=(26.d0*p*p + 120.d0*p*c - 120.d0*c*c)*p*p Bp3=(120.d0*c*c - 8.d0*p*p - 120.d0*p*c)*p*p*p*p Bp4=40.d0*(p - c)*p*p*p*p*p*p*c c Bp11=Bp1/s1 Bp22=Bp2/s3 Bp33=Bp3/s5 Bp44=Bp4/s7 Bp55=16.d0*p Bp=16.d0*p + 2.d0*s1 + Bp1/s1 + Bp2/s3 . + Bp3/s5 + Bp4/s7 c dAdr=2.d0/3.d0 + A1/r2 - 3.d0*A2/r4 c dApds=4.d0/3.d0 - Ap1/s2 - 3.d0*Ap2/s4 - 5.d0*Ap3/s6 c dBdr=4.d0 + B1/r2 - 3.d0*B2/r4 + 5.d0*B3/r6 c dBpds=2.d0 - Bp1/s2 - 3.d0*Bp2/s4 - 5.d0*Bp3/s6 . -7.d0*Bp4/s8 c du1dx=-4.d0*y*drrdx*(A+Ap + x*x*(B+Bp)/rr2)/rr5 . + (y/rr4)*((dAdr + x*x*dBdr/rr2)*drdx + (dApds + . x*x*dBpds/rr2)*dsdx) + (y/rr4)*(B + Bp)*dxrdx c term11=(1.d0/rr4 - 4.d0*y*drrdy/rr5) term12=(A+Ap + x*x*(B+Bp)/rr2) t1=A+Ap t2=x*x*(B+Bp)/rr2 t3=x t4=B+Bp t5=rr2 term1=term11*term12 term2=(y/rr4)*((dAdr+x*x*dBdr/rr2)*drdy + . (dApds + x*x*dBpds/rr2)*dsdy) term3= + (y/rr4)*(B + Bp)*dxrdy du1dy= term1 + term2 + term3 c du2dx=(1.d0/rr4 - 4.d0*x*drrdx/rr5)*(A+Ap + y*y* . (B+Bp)/rr2) + (x/rr4)*((dAdr+y*y*dBdr/rr2)*drdx + . (dApds + y*y*dBpds/rr2)*dsdx) + (x/rr4)*(B + Bp)* . dyrdx c du2dy=-4.d0*x*drrdy*(A+Ap + y*y*(B+Bp)/rr2)/rr5 . + (x/rr4)*((dAdr + y*y*dBdr/rr2)*drdy + (dApds + . y*y*dBpds/rr2)*dsdy) + (x/rr4)*(B + Bp)*dyrdy exx(jk)=du1dx*u eyy(jk)=du2dy*u exy(jk)=.5d0*(du1dy + du2dx)*u c end do return end c c function sgn(x) c c This function computes the sign of the supplied c quantity x c implicit real*8(a-h,o-z) sgn=1.d0 if(x .lt. 0.) sgn=-1.d0 if(x .eq. 0.) sgn=0.d0 return end c real function jsecond() c c This function does the timing of the fortran code c real junk, tarray(2) junk = etime(tarray) jsecond = tarray(1) return end