c     PROGRAM VC_SER.F
c
c	VIRTUAL_SOUTHERN_CALIFORNIA, serial version
c
c	Modified October 3, 2002
c
c	*************************************************************
c
c     General Description:
c
c	This is the Monte Carlo code that generates simulated,realistic
c	earthquakes on an arbitrary fault surface mesh.  The topology
c	of the fault mesh is defined by the user.  The stress Green's
c	functions are then computed in codes VC_Stress_Green.f
c	and VC_SG_Compress.f and then used as input to
c	the Virtual_California family of codes.  The stress Green's
c	functions are then used, together with a user-defined friction
c	model, as input to a initializer code VC_INIT_SER.f
c	The initializer code is run several times until the user sees
c	all initial transient effects have ceased.  Finally, this
c	code, VC_SER.f is used to generate simulated
c	earthquakes.
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, 2. Application
c		to Southern California, J. Geophys. Res., 93, 6255
c		(1988)
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	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 data and static stress
c				Green's function factors, updates
c				plate tectonic time, and contains
c				mcs loop
c
c	Subroutine TSUM()   --	Computes time factors for elastic
c				or viscoelastic stress Green's
c				functions at time t.  Must be
c				called prior to calling STRESS().
c
c	Subroutine STRESS() --	Computes elastic or viscoelastic
c				stress for the fault segments at time
c				t
c
c	Subroutine STREL()  --	Computes only the elastic stress
c				on segments.  Called during mcs
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.  Current program state variables read in from file
c		storing current state information from previous
c		code runs.  Variables are read in as single precision
c		variables, then converted to double precision
c		variables.  The double precision variables are:
c
c		TBIAS	- An unimportant dummy variable
c		hplate	- Elastic plate thickness in km
c		vplatx	- East-component of plate velocity in cm/yr
c		vplaty	- North-component of plate velocity in cm/yr
c		tauast	- Maxwell viscoelastic asthenosphere relaxation
c				time in yrs.
c		taufar	- Far field background shear stress in bars
c		bulkm	- Elastic bulk modulus in bars
c		amu	- Elastic shear modulus in bars
c		tmin	- Time at which present run begins
c		tstep	- Plate tectonic time step ("long time step")
c		ptrade	- Small parameter that quantifies the additive
c				nonlinear term (no longer used)
c		pdecy   - Small parameter that quantifies the additive
c				nonlinear term for each segment.
c		
c		
c		Note:  The following variables are single precision
c		versions of the above:
c
c		TBIS, hpl, vplx, vply, taua, tauf, bulkms, amuu,
c		tminn, tstepp, ptrad, pdec
c
c
c     This code takes the stress coefficients generated by 
c     strcoeff and generates synthetic earthquakes from them 
c     using a constant stress drop condition on the patches
c
c     This version of the code uses the garden variety CA
c     algorithm, but with precursory slip
c     to smooth the stress distribution.  A segment does
c     not heal until nflag sweeps have occurred during an event.
c
c     There is no aging or logarithmic friction dependence in
c     this code.
c
c	Uses fault-variable aseismic slip factor
c
c	This code can only be used to continue simulations after
c		VC_INIT_SER.F has been run
c
c	Best implementation of Karner-Marone aseismic slip to date
c
c	Second try at removing positive strdif's in creeping section
c		This time we use the normal stress change in dtauxy
c		as well as the shear stress change.
c
c
c	This code allows one to adjust the facrat() parameters.
c
c       
c       Modified to allow timing of code runs using implicit
c       function etime
c
c
c     ******************************************************
c
      implicit real*8(a-h,o-z)
      INCLUDE 'EQPARM_4.FOR'
      character*1 bell
      real*4 hpl,vplx,vply,taua,tauf,amuu,tminn,tstepp,ptrad
      real*4 rhofcc,dslipp,tslipp,taub,TBIS,dfr,timi
      real*4 cfr,slpv,delts,bulkms, slpdf, facsd
      real*4 pdec
      real*4 rand
      real*4 dslpT
      real*4 pdcy
      real*4 tsec_1, tsec_2, tsec_3, tsec_4, tdfsec, jsecond
      character resp3*1, rptime*1, rpts*1, respc*1
      character rsprs*1, respfb*1, respas*1
      character rpdec*1
      character fnstr*20, fnfric*20, fnmod*20
      character fnout*20, fndum*20, fnin*20
      dimension rhofcc(NF)
      dimension dslipp(NF,NTSTP), tslipp(NF,NTSTP)
      dimension dfr(NF), timpol(100)
      dimension cfric(NF), dfric(NF)
      dimension cfrtim(NF), cfrdif(NF), cfrred(NF)
      dimension signum(NF)
      dimension cfr(NF),slpv(NF),delts(NF)
      dimension txybac(NF), taub(NF)
      dimension depth(NF), rhogd(NF)
      dimension isltot(NTSTP,NF), displt(NF), dsltot(NF)
      dimension slpdef(NF)
      dimension Sxyvec(NF), Snvec(NF), Slpvel(NF)
      dimension tauxy(NF), taun(NF)
      dimension tstart(NF), dslip(NF,NTSTP), tslip(NF,NTSTP)
      dimension tfact(10,NF), islip(NF), iflag(NF)
      dimension CFF(NF)
      dimension timint(NF), timi(NF)
      dimension ifltt(NF), vcof(NF,NF), vcovec(NFF)
      dimension veldum(NF),vncof(NF,NF),vnvec(NFF)
      dimension slpdfd(NF), slpdf(NF)
      dimension Sxy(NFF),Sn(NFF),Sxym(NMF),Snm(NMF)
      dimension dsl(NF),gamma(NF)
      dimension sdrpL(NF), sdrpS(NF), strdrp(NF), sdrop(NF)
      dimension disnmL(NF), disnmS(NF)
      dimension pdecay(NF), pdec(NF)
      dimension dtauxy(NF)
      dimension cfrict(NF), cfricr(NF), Tlast(NF)
      dimension Strshr(NF), Strnml(NF)
      dimension Factor(NF)
      dimension aslrun(NF), dslrun(NF), frunsl(NF)
      dimension dslpT(NF,NTSTP)
      dimension facsdp(NF), facsd(NF)
      dimension slpdav(NF), slpav(NF)
      dimension strdfl(NF), indstr(NF), delCFF(NF)
      dimension STrate(NF), facrat(NF), ASrate(NF), Tr(NF)

      common/str/Sxy,Sn,Sxym,Snm,Sxyvec,Snvec
      external sgn, jsecond
c
c     ******************************************************
c
c
      print *,' '
      print *,' Enter the name of the existing (input) file:'
      read(*,'(a20)') fnin
      print *,' Enter a name for the output file:'
      read(*,'(a20)') fnout
c
c     Aseismic slip occurs
c
      respas = 'y'

      open(unit=30,name=fnout,status='new')
      close(30)
c
C     Sets up the run:
c
      open(unit=30,name=fnin,status='old')
c
c     Read in all the restart variables
c
      read(30,'(a20)') fnmod
      read(30,'(a20)') fnstr
c
c     fnmod: file with fault network topology
c     fnstr: file with fault stress coefficients
c
      read(30,*) nfault,hpl,vplx,vply,taua,tauf,amuu,
     .tminn,tstepp,ntime,TBIS,bulkms,ntm,ncycle,pdcy
      read(30,*) ptrad
      read(30,'(a1)') respfb

      TBIAS=dble(TBIS)
      hplate=dble(hpl)
      vplatx=dble(vplx)
      vplaty=dble(vply)
      tauast=dble(taua)
      taufar=dble(tauf)
      bulkm=dble(bulkms)
      amu=dble(amuu)
      tmin=dble(tminn)
      tstep=dble(tstepp)
      ptrade = dble(ptrad)
      pdecy = dble(pdcy)
c	
c     Current program state variables read in from file
c	storing current state information from previous
c	code runs.  Variables are read in as single precision
c	variables, then converted to double precision
c	variables.  The double precision variables are:
c
c	TBIAS	- An unimportant dummy variable
c	hplate	- Elastic plate thickness in km
c	vplatx	- East-component of plate velocity in cm/yr
c	vplaty	- North-component of plate velocity in cm/yr
c	tauast	- Maxwell viscoelastic asthenosphere relaxation
c			time in yrs.
c	taufar	- Far field background shear stress in bars
c	bulkm	- Elastic bulk modulus in bars
c	amu	- Elastic shear modulus in bars
c	tmin	- Time at which present run begins
c	tstep	- Plate tectonic time step ("long time step")
c	ptrade	- Small parameter that quantifies the additive
c			nonlinear term (no longer used)
c	pdecy   - Small parameter that quantifies the additive
c				nonlinear term for each segment.
c		
c     Note:  The following variables are single precision
c	versions of the above:
c
c	TBIS, hpl, vplx, vply, taua, tauf, bulkms, amuu,
c       tminn, tstepp, ptrad, pdec
c
      read(30,*) (timi(n), n=1,nfault)
      read(30,*) (taub(n), n=1,nfault)
      read(30,*) (delts(n), n=1,nfault)
      read(30,*) (cfr(n), n=1,nfault)
      read(30,*) (dfr(n), n=1,nfault)
      read(30,*) (slpdf(n), n=1,nfault)
      read(30,*) (slpv(n), n=1,nfault)
      read(30,*) (rhofcc(n), n=1,nfault)
      read(30,*) (islip(n), n=1,nfault)
c
c     converts to double
c
      do i=1,nfault
         timint(i) = dble(timi(i))
         txybac(i) = dble(taub(i))
         displt(i) = dble(delts(i))
         cfrdif(i) = dble(cfr(i))
         dfric(i)  = dble(dfr(i))
c     Computed on the fly:
         cfric(i)  = dfric(i) + cfrdif(i)

         slpdef(i) = dble(slpdf(i))
         slpvel(i) = dble(slpv(i))
         rhogd(i)  = dble(rhofcc(i))
         
      end do
      do i=1,nfault
         read(30,*) (dslipp(i,k), k=1,islip(i))
         do j=1,islip(i)
            dslip(i,j)=dble(dslipp(i,j))
         end do
      end do
c
      do i=1,nfault
         read(30,*) (tslipp(i,k), k=1,islip(i))
         do j=1,islip(i)
            tslip(i,j)=dble(tslipp(i,j)) + dble(TBIS)
         end do
      end do
c
c     State variables for each segment: as above, read in as single
c        precision and converted to double.
c     timint(i): recurrence interval (time) of segment i (corrected)
c     txybac(i): background stress on segment i
c     displt(i): total slip of segment i over all time steps
c     dfric(i): dynamics coeficient of friction segment i
c     cfric(i): static friction coefficient of segment i
c     cfrdif(i): diff between static & dynamic friction coeff.
c     slpdef(i): slip deficit segment i: (slip - Vt)
c     slpvel(i): slip velocity of segment i
c     rhogd(i): pressure in bars at the midpoint of segment i
c     dslip(i,j): amount of slip of segment i when it slipped for
c	 the jth time
c     tslip(i,j): time at which segment i slipped for the jth time
c
c     Note: The following variables are single precision 
c        versions of the above:
c
c     timi(i), taub(i), delts(i), cfr(i), dfr(i), slpdf(i), slpv(i), 
c     rhofcc(i), dslipp(i,j), tslipp(i,j)
c
      print *, ' '
c
      nttime=ntime-ntm
      ncycle=ncycle + 1
      ntm=int(sngl(40.d0*tauast)/tstep) + 1

      if(ntm .gt. ntime) ntm=ntime
c
      nthi=ntime + 1001

c      if(nthi .lt. NTSTP) then
c         read(30,*) ((isltot(i,j), i=1,nthi), j=1,NF)
c      else
         read(30,*) ((isltot(i,j), i=1,ntime), j=1,NF)
c      end if

      read(30,*) (slpav(n), n=1,nfault)

      read(30,*) (facsd(n), n=1,nfault)

      do i=1,nfault
         slpdav(i) = dble(slpav(i))
         facrat(i) = dble(facsd(i))
         facsdp(i) = facrat(i)*timint(i)
      end do

c
c     More state variables:
c     slpdav(i): average slip deficit
c     facrat(i): renormalized fraction of stable slip
c     facsdp(i): un-renormalized fraction of stable slip
c
c     Single precision correspondents:
c     slpav(i), facsd(i)
c
      close(30)

      print *, ' '
      print *, ' Finished reading input data --'
      print *, ' Hit <CR> to continue.'
      read(*,'(a1)') resp3
c
c     Redefine matrices
c
c     ntm is the number of time steps at the end we want
c     to use to reinitialize the calculation when we restart.
c
c     nti is the time step we want to start with
c
      nti=ntime - ntm
c
      do i=1,nfault
         isli=isltot(nti,i)
         islf=isltot(ntime,i)
         Tlast(i) = tslip(i,islf)
         islip(i)=islf-isli

         do j=1,islip(i)
            jsp=isli + j
            dslip(i,j)=dslip(i,jsp)
            tslip(i,j)=tslip(i,jsp)
         end do
      end do
c
      do i=1,nfault
         jis=islip(i) + 1
         do j=jis,NTSTP
            dslip(i,j)=0.d0
            tslip(i,j)=0.d0
         end do
      end do
c
      do i=1,ntm
         ijnt=nti + i
         do j=1,nfault
            isltot(i,j)=isltot(ijnt,j) - isltot(nti,j)
         end do
         isltot(i,NF)=isltot(ijnt,NF)
      end do
c
c     Read the fault coefficients
c
      open(unit=10,name=fnstr,status='old')
      read(10,'(a20)') fndum
      read(10,*) nfault,mend,hplate,vplatx,vplaty,bulkm
      imax=nfault*nfault
      imaxm=imax*mend
      read(10,*) (timpol(n), n=1,mend)
      read(10,*) (depth(n), n=1,nfault)
      read(10,*) (Sxy(n), n=1,imax)
      read(10,*) (Sn(n), n=1,imax)
      read(10,*) (Sxym(n), n=1,imaxm)
      read(10,*) (Snm(n), n=1,imaxm)
      read(10,*) (gamma(n), n=1,nfault)
      read(10,*) (veldum(n), n=1,nfault)
      read(10,*) (veldum(n), n=1,nfault)
      read(10,*) (veldum(n), n=1,nfault)
      read(10,*) (vcovec(n), n=1,imax)
      read(10,*) (vnvec(n), n=1,imax)
      close(10)
c
c     Fault coefficients:
c     timpol(j): Time constant for term j out of the 8 terms in the
c	viscoelastic sum of decaying exponentials for each fault segment.  
c       Units of yrs.
c     depth(i): depth of midpoint of segment i
c     Sxy(k): Shear stress elastic Green's function for	influence 
c       of segment i (=k/NF) on segment j =(k#NF). Units are in 
c       inverse cm (cm^-1). Will be multiplied by amu prior to using 
c       in Virtual California.
c	Note that the N*N coefficients are stored by columns as usual, 
c       with first column being source fault i, all N target faults j.
c     Sn(k): Normal stress elastic Green's function for influence of 
c       segment i (=k/NF) on segment j =(k#NF). Units are in inverse 
c       cm (cm^-1). Will be multiplied by amu prior to using in 
c       Virtual California.
c	Note that the N*N coefficients are stored by columns as usual, 
c       with first column being source fault i, all N target faults j.
c     Sxym(k): Shear stress viscoelastic Green's function for influence 
c       of segment i (=k/NF) on segment j (=k#NF). Units are in inverse 
c       cm (cm^-1).  Will be multiplied by amu prior to using in 
c       Virtual California.
c       Note that the N*N*mend coefficients are stored by columns as usual, 
c       with first column being source fault i, all N target faults j, all
c	mend images.
c     Snm(k): Normal stress viscoelastic Green function for influence of 
c       segment i on segment j. Units are in inverse cm (cm^-1).  Will be 
c       multiplied by amu prior to using in Virtual California.
c	Note that the N*N*mend coefficients are stored by columns as usual, 
c       with first column being source fault i, all N target faults j, all
c	mend images.
c     gamma(NF): areas of segments in km^2
c     Sxyvec(NF): represents the total shear stress change per year, in 
c       inverse yr, (yr^-1) (i.e., normalized by amu in bars), due to
c	the loading velocity on all the faults.
c     Snvec(NF): Represents the total normal stress change per year, in 
c       inverse yr, (yr^-1) (i.e., normalized by amu in bars), due to
c	the loading velocity on all the faults.
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     vcovec(NFF): Used in computing svec(NF).  Represents the total 
c       influence of a unit of slip on fault segment j on the shear stress 
c       change on fault segment i.
c     vnvec(NFF) - Used in computing snvec(NF).  Represents the total 
c       influence of a unit of slip on fault segment j on the normal stress 
c       change on fault segment i.
c
c
c     Multiply by amu. Also, convert to km (multiply by 1.e-5 km/cm)
c
      do i=1,imax
         Sxy(i)=amu*Sxy(i)*1.d-5
         Sn(i)=amu*Sn(i)*1.d-5
         vcovec(i)=amu*vcovec(i)*1.d-5
         vnvec(i)=amu*vnvec(i)*1.d-5
      end do
      do i=1,imaxm

        Sxym(i)=amu*Sxym(i)*1.d-5
        Snm(i)=amu*Snm(i)*1.d-5
c
c	***************************
c	No viscoelasticity for now
c	***************************
c
         Sxym(i)=0.d0
         Snm(i)=0.d0

      end do
      do i=1,nfault
         Sxyvec(i) = amu*Sxyvec(i)*1.d-5
         Snvec(i) = amu*Snvec(i)*1.d-5
      end do
      indv=0
      do i=1,nfault
         do j=1,nfault
            indv=indv + 1
            vcof(j,i) = vcovec(indv)
            vncof(j,i) = vnvec(indv)
         end do
      end do
c
      do i=1,nfault
         ind=(i-1)*nfault + i
         disnmL(i) = slpvel(i)*timint(i)
         disnmS(i) =  facsdp(i)*disnmL(i)
         sdrpL(i) = Sxy(ind)*disnmL(i)
         sdrpS(i) = facsdp(i)*sdrpL(i)
c
c	Depth is depth to midpoint of fault.  rhogd is pressure at
c	that depth in bars
c
         rhogd(i)  =  depth(i) * 100.d0
      end do
c
      print *, ' '
      print *, ' Finished reading stress interaction coefficients --'
      print *, ' Hit <CR> to continue.'
      read(*,'(a1)') resp3
c
c     Enter values for time steps
c
      tstpo = tstep

      print *,' '
      print *,' Current time step is (years):', tstep
      print *,' Change time step? (y/n)'
      read *, rpts
      if(rpts .eq. 'y') then
         print *,' Enter time step (years)'
         read *, tstep
      end if

      NTMAX=ntm
      NTMAX=NTSTP-NTMAX

      print *,' '
      print *,' First time value will be', tmin
      print *,' '
      print *,' Maximum number of time steps is',NTMAX
      print *,' How many time steps?'
      read *, ntime
c
c     Start into time loop
c
      ktest = ntime / 10
      t = tmin - tstep
      tpres = tmin
c
      print *,' '
      print *,' Current Time (years):  ',tpres
      print *,' '
      print *,' '
      print *,' ******************************'
      print *,' '
      print *,' BEGIN TIMED PART OF RUN...'
      print *,' '
      print *,' ******************************'
      print *,' '
      print *,' '
c
      KKOUNT=0
c
c     nflag is number of mcs for a segment to heal
c
      nflag = 50
c
c     Naftsl is the afterslip, or healing, parameter...it determines
c		how many sweeps after an unstable slip the healing is
c		enforced.
c
      Naftsl = 5

      do ir = 1,nfault
         signum(ir)=sgn(slpvel(ir))
         dslrun(ir) = 0.d0
         aslrun(ir) = 0.d0
c
c	dslrun(i): total of unstable slip of segment i in this run
c	aslrun(i): total increment of stable slip of segment i in this run
c
      end do

      Nsegmx = 0

      do ir=1,nfault
         cfrict(ir) = dfric(ir) + cfrdif(ir)
         pdecay(ir) = ptrade
      end do
c
c	Begin timing the fortran program
c
      tsec_1 = jsecond()

      print *, ' Hit <CR> to continue. '
      read(*,'(a1)') resp3

      do it=1,ntime
c
c     Increment time
c
         t = t + tstep
         
         KKOUNT=KKOUNT + 1
c
c     Updates slip deficits
c
         do ir=1,nfault
            slpdef(ir)= slpdef(ir) - slpvel(ir) * tstep
         end do
c
c     Calculate slpdfd by subtracting the recent displacement from
c          slpdef
c
         do i=1,nfault
            slpdfd(i) = slpdef(i)
            do j=1,islip(i)
               tbart = (t - tslip(i,j))/tauast
c
c     See analogous statments in TSUM
c
               if(tbart .le. 10.) then
                  slpdfd(i) = slpdfd(i) - dslip(i,j)
               end if
            end do
         end do
c
        call TSUM(t,tauast,dslip,tslip,tfact,slpvel,slpdef,
     &  	timpol,slpdfd,nfault,mend,islip)
c
c     TSUM returns (updates only) tfact
c
        call STRESS(t,tauxy,taun,rhogd,tfact,slpvel,
     &  	nfault,mend,slpdef)
c
c     Calculate effective shear and normal stress, including
c        background stress (txybac) and stabilizing term (sinhfc)
c
        do ir = 1,nfault
           Strnml(ir) = - taun(ir)
           Strshr(ir) = tauxy(ir)* signum(ir) + txybac(ir)
        end do
c
c	Now calculate stress difference
c
	do ir = 1,nfault
           sinhfc = signum(ir) * dabs(sdrpL(ir))*dsinh(pdecay(ir)*
     &		(slpdef(ir) - slpdav(ir)) )
           CFF(ir) = Strshr(ir) - cfrict(ir)*Strnml(ir) - sinhfc
	end do

c
c	If first time through for real dynamics,
c		renormalize txybac so that one site is at failure
c

	if (it .eq. 1 .and. ncycle .eq. 1) then
c
           strmax = -1.d10
           do ir = 1,nfault
              if(CFF(ir) .gt. strmax) strmax = CFF(ir)
           end do

           do ir = 1,nfault
              txybac(ir) = txybac(ir) - strmax
              CFF(ir) = CFF(ir) - strmax
           end do

           print *,' INITIAL CFFs (bars) --'
           print *, (CFF(n), n=1,nfault)
           
           print *, ' '
           print *, ' Hit <CR> to continue. '
           read(*,'(a1)') resp3

        end if
c
c     Allow precursory slip to occur.  Do not turn failure flag on
c		(iflag) since this precursory slip is stable.
c
        do ir =1,nfault	
c
c	Note that strdif =< 0.
c
           strttt = signum(ir) * sdrpL(ir)
           
           dsl(ir) = 0.d0
    	
           if (it .gt. 1) then
c
c	With the following formulation of dsl, dsl ~ 0 when stress is
c		low, therefore stress increases rapidly
c		on the segment.  But when stress is near zero,
c		Factor ~ 1 and dsl increases maximally rapidly
c		[it is at most facsdp(ir) * V * tstep / timint(ir) ]

c
c	***************	THIS IS THE ASEISMIC SLIP PART ************
c
c	If facsdp .ne. 1, we relieve some stress via aseismic
c		 slip at each time step
c	We don't do this for the creeping part...creep is considered
c		to be small earthquakes
c
              tstfac = dabs(1.d0 - facsdp(ir))
              if (tstfac .gt. .00001d0) then

                 ind=(ir-1)*nfault + ir


                 delCFF(ir) = CFF(ir) - (signum(ir) * sdrpL(ir) )

                 if(delCFF(ir) .gt. 0.d0) then
                    dsl(ir)  =  - delCFF(ir)*signum(ir) * facrat(ir)
     $                   * tstep / Sxy(ind)
                 end if

              end if
c
c	*******************************************************
c
              
           end if

           slpdef(ir) = slpdef(ir) + dsl(ir)
           displt(ir) = displt(ir) + dsl(ir)
           slpdfd(ir) = slpdfd(ir) + dsl(ir)
           aslrun(ir) = aslrun(ir) + dsl(ir)
           dsltot(ir) = 0.d0
           dslpT(ir,it) = sngl(dsl(ir))
c
c	dsl(i): increment of slip - here, stable - of segment i
c	displt(i): total slip of segment i over all time steps since begining
c	aslrun(i): total increment of stable slip of segment i in this run
c	dsltot(i): total increment of unstable slip of segment i
c
        end do
c
c     Having updated the slip with some stable slip,
c           we now need to update the current values
c	    of stress
c
        call TSUM(t,tauast,dslip,tslip,tfact,slpvel,slpdef,
     &  	timpol,slpdfd,nfault,mend,islip)
c
        call STRESS(t,tauxy,taun,rhogd,tfact,slpvel,
     &  	nfault,mend,slpdef)
c
c	iflag(i) keeps track of how many total times segment i has failed
c
        do ir=1,nfault
           iflag(ir)=0
           delCFF(ir) = 0.d0
        end do

        respc='c'
c
c	This is the failure loop at time t
c
        Nsweep = 0
	
        do while(respc .eq. 'c')

           Nsweep = Nsweep + 1
c
c	ifltt keeps track of which faults have just failed
c
           do ir=1,nfault
              ifltt(ir)=0
           end do
c
c	Check to see which faults have failed.
c       Also, the quantity jp is the stage variable; at any stage of the
c       calculation, it tells you how many faults have triggered
c
c
c	Figure out if any sites are at failure
c
           respc = 's'
           nflagn = 0
	
           do ir = 1,nfault
c
c	Define stress difference - computes Coulomb failure function (CFF) 
c         for each fault - at this stage of the while loop...
c	  some stresses may have changed due to slip on
c		that fault or other faults.
c
              sinhfc = signum(ir) * dabs(sdrpL(ir))* dsinh(pdecay(ir)*
     &             (slpdef(ir) - slpdav(ir)) )
              Strnml(ir) = - taun(ir)
              Strshr(ir) = tauxy(ir)* signum(ir) + txybac(ir)
              CFF(ir) = Strshr(ir) - cfrict(ir)*Strnml(ir) - sinhfc	
c
c	Increment iflag if strdif .GT. 0.d0
c
              if(CFF(ir) .gt. 0.d0) then
                 
                 nflagn = 1
                 iflag(ir) = iflag(ir) + 1
c
c	Also, increment iflag if site *ir* ever failed during
c		this time step
c
              else
                 if(iflag(ir) .gt. 0) iflag(ir) = iflag(ir) + 1
              end if
           end do	
c
c	Initialize number of failed sites to zero ( jp = 0 )
c
           jp = 0
c
c	There must be at least one new rupture in this or the previous
c		time step.  Check to see if this condition is true.
c		If inflag = 1, there is at least one new rupture.
c
           inflag = 0
           do ir =1,nfault
c     if (iflag(ir) .eq. 1 .or. iflag(ir) .eq. 2) inflag = 1
              if (iflag(ir) .ge. 1 .and. iflag(ir) .le. Naftsl) 
     &             inflag = 1
           end do
c
c	Only if inflag = 1 or 2, i.e., a new failure existed during
c		last two time steps do we proceed
c
c	NOTE:  THIS REQUIREMENT IS PROBABLY RESPONSIBLE FOR
c		THE FACT THAT THE CREEPING SECTION DOES
c		NOT ALWAYS FULLY QUIESCE FOR LARGE TIME STEPS
c
c	ANOTHER NOTE: WE HAVE NOW REPLACED THIS inflag STATEMENT
c		WITH THE Naftsl, OR *AFTERSLIP* PARAMETER
c
           if (inflag .eq. 1) then
c
c	Check for maximum number of sweeps of a given site
c
              do ir = 1,nfault
                 if(iflag(ir) .gt. 0 .and. iflag(ir) .le. nflag) then
                    respc = 'c'
                    jp = jp + 1
                    ifltt(jp) = ir
c
c     ifltt(jp): contains identifier (#) of jpth segment to fail
c
                 end if
              end do
           end if
c
c	Calculate slip increment dsl
c
c
c	New Failure Algorithm:  This is a CA-jump algorithm
c
c	A fault stays in the loop for nflag times before it heals
c
           do ir = 1,nfault
              dsl(ir) = 0.d0
           end do

           if (jp .gt. 0) then
	
              do ijr=1,jp
                 ikr = ifltt(ijr)
                 ind=(ikr-1)*nfault + ikr
c
c	Allow for triggered sites going above the failure
c		threshold momentarily.  Note that CFF(*)
c		is always a postive number, whereas sdrpL(*)
c		is negative for positive slip (right lateral)
c		and positive for negative slip (left lateral)
c		signum(*) = +1 for right lateral, -1 for
c		left lateral
c
                 strdrp(ikr) = - CFF(ikr) * signum(ikr) + sdrpL(ikr)
c
c	If segment has previously slipped during this event, do not
c		allow healing for nflag time steps..i.e., the stress
c		drop for subsequent events is stress increment in
c		the same sense as the stress buildup
c
                 if(iflag(ikr) .gt. 1) then
c
c	If a site has failed, we want to keep the stress at or below
c		the value of sdrpL.  So we want to compute change in
c		CFF once the site has failed...this will be delCFF(*).
c
                    delCFF(ikr) = CFF(ikr) - (signum(ikr) * sdrpL(ikr))
                    strdrp(ikr) = 0.d0
    	
                    if(delCFF(ikr) .gt. 0.d0) then
                       strdrp(ikr) = - delCFF(ikr)*signum(ikr)
                    end if   	
    	
                 end if
c
c	End of part where we determine static stress drop.  Now compute
c		nominal slip
c
                 dsl(ikr)  =  strdrp(ikr)/ Sxy(ind)
c
c	dsl(i): increment of slip - here, unstable - of segment i
c	strdrp(i): variable static stress drop of segment i
c
              end do
           end if
c
c	Increment slip
c
c
c	Introduce Noise.  +/- 5% Random Jumps.  This lets system
c		explore phase space.  Change to +/- 25% jumps
c
           do ir = 1,nfault
              dsl(ir)  =  dsl(ir)
c     &		* (1.d0 + .003d0*  (2.d0*(.5d0 - dble(rand())) ) )
     &             * (1.d0 + .050d0*  (2.d0*(.5d0 - dble(rand())) ) )
c     &		* (1.d0 + .010d0*  (2.d0*(.5d0 - dble(rand())) ) )
c     &		* (1.d0 + .25d0*  (2.d0*(.5d0 - dble(rand())) ) )
           end do
c
c	Adjust elastic stresses
c
           do ik = 1,nfault
              dsltot(ik) = dsltot(ik) + dsl(ik)
              displt(ik) = displt(ik) + dsl(ik)
              dslrun(ik) = dslrun(ik) + dsl(ik)
              slpdef(ik) = slpdef(ik) + dsl(ik)
c
c	dsltot(i): total increment of unstable slip of segment i
c	displt(i): total slip of segment i over all time steps since begining
c	dslrun(i): total of unstable slip of segment i in this run
c
           end do
c
c	Adjust tauxy and taun due to elastic slip
c
           call STREL(tauxy,taun,dsl,nfault)
c
c	End do for "do while"
c
c	Stop if too many sweeps
c
           if(Nsweep .ge. nflag) respc = 's'
c
        end do
c
c	Store total slip
c
        do ir = 1,nfault
           dslpT(ir,it) =  sngl(dsltot(ir)) + dslpT(ir,it)
        end do
c
        do ik =1,nfault

           denrun = aslrun(ik) + dslrun(ik)
           if (dabs(denrun) .gt. 0.d0) then
              frunsl(ik) = dslrun(ik) /denrun
           end if
        end do

        nf_hlf = nfault/2

        if(KKOUNT .eq. 25) then
           print *, ' '
           print *,' TOTAL DISPLACEMENT DURING THIS RUN (cm) --'
           print *, (dslrun(n), n=1,nfault)
           print *,' '
           print *,' FRACTION OF SLIP THAT IS UNSTABLE --'
           print *, (frunsl(n), n=1,nfault)

           KKOUNT=0
        end if
c
        do ir = 1,nfault
           if(iflag(ir) .gt. 0) then
              islip(ir)=islip(ir) + 1
              is=islip(ir)
              dslip(ir,is)=dsltot(ir)
              tslip(ir,is)=t
              Tlast(ir) = t
c
c     islip(i):   # of times segment i slipped
c     dslip(i,j): total unstable slip of segment i when it slipped for 
c        the jth time
c     tslip(i,j): time at which segment i slipped for the jth time
c     Tlast(i): time at which segment i slipped for the last time
c
           end if
        end do
c
c     Store the time-status vectors
c
        itmntm=it + ntm
        do jj=1,NF
           isltot(itmntm,jj)=islip(jj)
        end do
        isltot(itmntm,NF)=ns
c
c	Calculate number of segments that slipped
c
	Nsegsl = 0
	
	do ir = 1,nfault
           if (iflag(ir) .gt. 0) Nsegsl = Nsegsl + 1
	end do
	
        if (it .gt. 10) then
           if(Nsegsl .gt. Nsegmx) then
              Nsegmx = Nsegsl
              itseg = it
           end if
        end if

        nzero = 0
        do ijn=1,nfault
           Tlslp = t - Tlast(ijn)
           if (Tlslp .gt. 10000.d0) nzero = nzero + 1
        end do
c
c     Searches maximum, minimum values of CFF
c
        CFFmax = -1.d20
        CFFmin =  1.d20
        CFFmen =  0.d0
        CFFstd =  0.d0
        jinds = 0
        do ir = 1,nfault
           if(CFF(ir) .gt. CFFmax) CFFmax = CFF(ir)
           if(CFF(ir) .lt. CFFmin) CFFmin = CFF(ir)
           CFFmen = CFFmen + CFF(ir)
           if (CFF(ir) .gt. 0.d0) then
              jinds = jinds + 1
              indstr(jinds) = ir
           end if
        end do

        CFFmen = CFFmen / dfloat(nfault)

        do ir = 1,nfault
           CFFstd = CFFstd + (CFF(ir)- CFFmen)*
     $          (CFF(ir) - CFFmen)
        end do

        CFFstd = dsqrt(CFFstd/dfloat(nfault))
	
	print *, ' '
        print *,'We have completed time step: ', it
        print *,'  ==>  Total Time Steps Requested: ', ntime
        print *,'  ==>  Current Time: ', t
        print *,'  ==>  Number of Segments that Slipped: ', Nsegsl
        print *,'  ==>  Maximum Number of Segments: ', Nsegmx
        print *,'           on Time Step: ', itseg
        print *,'  ==>  Number of Sweeps: ', Nsweep
	print *,'  ==>  MAX, MIN CFF: ', CFFmax, CFFmin
	print *,'  ==>  MEAN CFF: ', CFFmen
        print *,'  ==>  Number Locked in Last 10000 years: ',nzero
	if (jinds .gt. 0) then
           print *,'  ==>  Segments with strdif > 0 : ',
     &		(indstr(ijk), ijk = 1,jinds)
	end if
c
c	*******************************************************
c
c
c     Write out output file every .1 * ntime steps
c
	ktime = (it / ktest) * ktest
	
	if (ktime .eq. it) then
           tmin_1 = tmin + dfloat(it)*tstep
           itime = itmntm

           ntlo=itime + 1
           nthi=ntlo + 1000
           if(nthi .lt. NTSTP) then
              do i=ntlo,nthi
                 do j=1,NF
                    isltot(i,j)=isltot(itime,j)
                 end do
              end do
           end if
c
c     Write out slip histories:
c        converts to single...
c
           hpl=sngl(hplate)
           vplx=sngl(vplatx)
           vply=sngl(vplaty)
           TBIS=sngl(TBIAS)
           taua=sngl(tauast)
           tauf=sngl(taufar)
           amuu=sngl(amu)
           tminn=sngl(tmin_1)
           tstepp=sngl(tstep)
           bulkms=sngl(bulkm)
           ptrad = sngl(ptrade)
           pdcy = sngl(pdecy)
c
           do i=1,nfault
              pdec(i)=sngl(pdecay(i))
              facsd(i) = sngl(facrat(i))
              timi(i)=sngl(timint(i))
              taub(i)=sngl(txybac(i))
              rhofcc(i)=sngl(rhogd(i))
              cfr(i)=sngl(cfrdif(i))
              dfr(i)=sngl(dfric(i))
              slpdf(i)=sngl(slpdef(i))
              slpv(i)=sngl(slpvel(i))
              delts(i)=sngl(displt(i))
              slpav(i) = sngl(slpdav(i))
              do j=1,islip(i)
                 dslipp(i,j)=sngl(dslip(i,j))
                 tslipp(i,j)=sngl(tslip(i,j))-sngl(TBIAS)
              end do
           end do
c
           open(unit=30,name=fnout,status='old')
c
c        ...and writes out
c
           write(30,'(a20)') fnmod
           write(30,'(a20)') fnstr
           write(30,*) nfault,hpl,vplx,vply,taua,tauf,amuu,
     .          tminn,tstepp,itime,TBIS,bulkms,ntm,ncycle,pdcy
           write(30,*) ptrad
           write(30,'(a1)') respfb
           write(30,*) (timi(n), n=1,nfault)
           write(30,*) (taub(n), n=1,nfault)
           write(30,*) (delts(n), n=1,nfault)
           write(30,*) (cfr(n), n=1,nfault)
           write(30,*) (dfr(n), n=1,nfault)
           write(30,*) (slpdf(n), n=1,nfault)
           write(30,*) (slpv(n), n=1,nfault)
           write(30,*) (rhofcc(n), n=1,nfault)
           write(30,*) (islip(n), n=1,nfault)
c
           do i=1,nfault
              write(30,*) (dslipp(i,k), k=1,islip(i))
           end do
c
           do i=1,nfault
              write(30,*) (tslipp(i,k), k=1,islip(i))
           end do
c
c           if(nthi .lt. NTSTP) then
c              write(30,*) ((isltot(i,j), i=1,nthi), j=1,NF)
c           else
              write(30,*) ((isltot(i,j), i=1,itime), j=1,NF)
c           end if
           write(30,*) (slpav(n), n=1,nfault)
           write(30,*) (facsd(n), n=1,nfault)
           write(30,*) ((dslpT(i,j), i=1,nfault), j=1,it)
	   close(30)
c
	end if
c
c     "End do" for the time loop
c
      end do
c
c     Finish timing the fortran code
c
      tsec_2 = jsecond()
      tdfsec = tsec_2 - tsec_1
c
      print *, ' '
      print *, ' Duration of time-stepping loop in cpu secs: ', tdfsec
c
      stop
      end
c
      subroutine TSUM(t,tauast,dslip,tslip,tfact,slpvel,slpdef,
     &	timpol,slpdfd,nfault,mend,islip)
c
c     INTERFACE:
c     t: current time
c     tauast: Maxwell viscoelastic asthenosphere relaxation time in yrs.
c     dslip(i,j): amount by which segment i slipped for the jth time 
c     tslip(i,j): time at which segment i slipped for the jth time
c     tfact(i,j): total time factor for the ith term of the Green's function 
c                 expansion for segment j 
c     slpvel(i): slip velocity of segment i
c     slpdef(i): slip deficit of segment i
c     timpol(i): ...for the ith term of the expansion
c     slpdfd(i): slip deficit of segment i, corrected for recent displacements
c     nfault: number of segments
c     mend: # of terms in the Green's function expansion
c     islip(i): running total of number of times segment i has slipped
c
c     RETURN VALUE:
c     tfact(i,j): see above
c
c     DESCRIPTION:
c     Computes  the multiplicative summed time factors for the entire 
c     slip history of a given fault for elastic or viscoelastic stress 
c     Green's functions at time t.
c     Must be called prior to calling STRESS()
c
      implicit real*8(a-h,o-z)
      INCLUDE 'EQPARM_4.FOR'
      dimension tfact(10,NF), dslip(NF,NTSTP),
     .tslip(NF,NTSTP), islip(NF), slpdef(NF)
      dimension Sxy(NFF),Sn(NFF),Sxym(NMF),
     .Snm(NMF), Sxyvec(NF), Snvec(NF)
      dimension slpvel(NF), timpol(100)
      dimension slpdfd(NF)
      common/str/Sxy,Sn,Sxym,Snm,Sxyvec,Snvec
c
c     Enter fault loop
c
      do if=1,nfault
c
c     For each fault (if), compute all the normalized slip times
c     and then for each term of the Green's function expansion (m), 
c     calculate and sum the time factors to get one time factor for each m.
c
         do m=1,mend
            tfact(m,if)=0.d0
         end do
c
         jend=islip(if)
c
         if(jend .gt. 0)then
c
            do j=1,jend
               tbart=(t-tslip(if,j))/(tauast)
c          
               if(tbart .le. 10.d0) then
c
                  do m=1,mend
                     tbar=(t-tslip(if,j))/(timpol(m)*tauast)
                     tfa=1.d0 - dexp(-tbar)
                     tfact(m,if)=tfact(m,if) + tfa*dslip(if,j)
                  end do
               end if
            end do
         end if
      end do
c
      do if=1,nfault
         do m=1,mend
            tfact(m,if)=tfact(m,if) + slpdfd(if)
         end do
      end do
c
      return
      end      
c

      subroutine STRESS(t,tauxy,taun,rhogd,tfact,slpvel,
     &	nfault,mend,slpdef)
c
c     INTERFACE:
c     t: current time
c     tauxy(i): shear stress on segment i
c     taun(i): normal stress on segment i
c     rhogd(i): pressure in bars at midpoint of segment i
c     tfact(i,j): total time factor for the ith term of the Green's function 
c                 expansion for segment j 
c     slpvel(i): slip velocity of segment i
c     nfault: number of segments
c     mend: # of terms in the Green's function expansion
c     slpdef(i): slip deficit of segment i
c
c     RETURN VALUE:
c     tauxy(i): see above
c     taun(i): see above
c
c     DESCRIPTION:
c
c     This subroutine computes the stresses at the time t
c
      implicit real*8(a-h,o-z)
      INCLUDE 'EQPARM_4.FOR'
      dimension tauxy(NF),taun(NF),tfact(10,NF)
      dimension Sxy(NFF),Sn(NFF),Sxym(NMF),
     .Snm(NMF),Sxyvec(NF),Snvec(NF)
      dimension slpdef(NF), rhogd(NF)
      dimension slpvel(NF)
      common/str/Sxy,Sn,Sxym,Snm,Sxyvec,Snvec
      external sgn
c
c     Compute stresses
c
      do if=1,nfault
         tauxy(if)=0.d0
         taun(if)=-rhogd(if)
      end do
c
c     Couples segment ic with all other segments ir
c
      do ic=1,nfault
         ind=(ic-1)*nfault
         do ir=1,nfault
            id=ind + ir
            tauxy(ir)=tauxy(ir) + slpdef(ic)*Sxy(id)
            taun(ir)=taun(ir) +   slpdef(ic)*Sn(id)
         end do
      end do
c
      do ic=1,nfault
         do ir=1,nfault
            indm=(ic-1)*mend*nfault + (ir-1)*mend
            do m=1,mend
               idm=indm + m
               tfat=tfact(m,ic)
               tauxy(ir)=tauxy(ir) + tfat*Sxym(idm)
               taun(ir)=taun(ir) + tfat*Snm(idm)
            end do
         end do
      end do

      return
      end
c
      subroutine STREL(tauxy,taun,dsl,nfault)
c
c     INTERFACE:
c     tauxy(i): shear stress on segment i
c     taun(i): normal stress on segment i
c     dsl(i): precursory slip of segment i
c     nfault: number of segments
c
c     RETURN VALUE:
c     tauxy(i): see above
c     taun(i): see above
c
c     DESCRIPTION:
c
c     This subroutine adjusts the elastic stresses
c
      implicit real*8(a-h,o-z)
      INCLUDE 'EQPARM_4.FOR'
      dimension tauxy(NF), dsl(NF), taun(NF)
      dimension ifltt(NF)
      dimension slpdef(NF), slpdav(NF), pdecay(NF),signum(NF)
      dimension slpdf2(NF), sdrpL(NF)
      dimension Sxy(NFF),Sn(NFF),Sxym(NMF),
     &	Snm(NMF),Sxyvec(NF),Snvec(NF)
      common/str/Sxy,Sn,Sxym,Snm,Sxyvec,Snvec
c
c     Compute changes in elastic component of stress
c
      do ic=1,nfault
         ind=(ic-1)*nfault
         do ir=1,nfault
            id=ind + ir
            tauxy(ir)=tauxy(ir) + dsl(ic)*Sxy(id)
            taun(ir)=taun(ir) + dsl(ic)*Sn(id)
         end do
      end do
c
      return
      end
c
      function sgn(x)
c
c     Returns the sign of x (0 if x=0)
c
      implicit real*8(a-h,o-z)
      sgn=-1.d0
      if(x .eq. 0.d0) sgn=0.d0
      if(x .gt. 0.d0) sgn=1.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