c     PROGRAM VC_INIT_SER.F
c	
c	***	***	***	***	***	***	***
c
c	THIS IS THE CODE THAT ANNEALS THE SYSTEM INTO THE PROPER
c		WELL ON THE DYNAMICAL LANDSCAPE.  
c
c	Each segment can only fail ONCE on a time step to prevent
c		sympathetic slip that may tend to diverge with
c		time.
c
c	This code has a simplified dynamics in which the slip
c		displacement adapts.  There is no unstable slip.
c
c	See comment header file of VD_SER.f for details of
c		program variables and their usage.
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, slpav
      real*4 rand
      real*4 dslpT
      real*4 pdcy
      character resp3*1, rptime*1, rpts*1, respc*1
      character rpdec*1, rsncyc*1
      character rsprs*1, respfb*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), iflagp(NF), iflagn(NF)
      dimension strdif(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)
      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 dslsm(NF), dslslp(NF), dslint(NF), strdf0(NF)
      dimension slpstr(NF), sdrpLT(NF)
      dimension slpdav(NF), slpav(NF), pfactr(NF)
      dimension txybc(NF), slpdf0(NF), strdfl(NF)
      dimension delstr(NF)

      common/str/Sxy,Sn,Sxym,Snm,Sxyvec,Snvec
      external sgn
c
c     ******************************************************
c
c     Query for a restart
c
      print *,' '
      print *,' Do you want to restart'
      print *,' an existing calculation? (y/n)'
      read(*,'(a1)') rsprs
      if(rsprs .eq. 'y') then
         print *,' '
         print *,' Enter name of the existing (input) file:'
         read(*,'(a20)') fnin
         print *,' Enter name for the output file:'
         read(*,'(a20)') fnout
         
         open(unit=30,name=fnout,status='new')
         close(30)
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
         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
         do i=1,nfault
            if (islip(i) .gt. 0) then
               read(30,*) (dslipp(i,k), k=1,islip(i))
               do j=1,islip(i)
                  dslip(i,j)=dble(dslipp(i,j))
               end do
            end if
         end do
c
         do i=1,nfault
            if (islip(i) .gt. 0) then
               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 if
         end do
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
            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))
            cfric(i) = dfric(i) + cfrdif(i)
            slpdef(i)=dble(slpdf(i))
            slpvel(i)=dble(slpv(i))
            rhogd(i)=dble(rhofcc(i))
            slpdav(i) = dble(slpav(i))
            facsdp(i) = dble(facsd(i))
            pdecay(i) = pdecy
         end do
         close(30)

	 print *, ' '
	 print *, ' Finished reading in input data.'
         print *, ' Hit Carriage Return 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)
c
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, 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
c
c	***************************
c	No viscoelasticity for now
c	***************************
c
c        Sxym(i)=amu*Sxym(i)*1.d-5
c        Snm(i)=amu*Snm(i)*1.d-5
c
            Sxym(i)=0.d0
            Snm(i)=0.d0
         end do
c
         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 Greens function'
	 print *, ' and fault data.'
c
         print *, ' Hit Carriage Return to Continue'
         read(*,'(a1)') resp3
c
c     This is the "else" that goes with rsprs='y': first run
      else
c
c     Ask for coefficient file name
c
         print *,'  Enter the name of the file containing the'
         print *,'  stress Greens function coefficients:'
         read(*,'(a20)') fnstr
c
c     Open the coefficient file
c
         open(unit=10,name=fnstr,status='old')
         read(10,'(a20)') fnmod
         read(10,*) nfault,mend,hplate,vplatx,vplaty,bulkm
         read(10,*) (timpol(n), n=1,mend)
         read(10,*) (depth(n), n=1,nfault)
c
c     Description of variables: see above
c
c
c     Ask for name for the output file
c
         print *,' Enter a name for the output file:'
         read(*,'(a20)') fnout

         open(unit=30,name=fnout,status='new')
         close(30)
c
c     Ask for name for the friction file
c
         print *,' Enter the name of the friction file:'
         read(*,'(a20)') fnfric
c
c     No feedback in this code!!!
c
         open(unit=20,name=fnfric,status='old')
c
         read(20,'(2g20.10)') amu, pfacl,cfrcl
         do i=1,nfault
            read(20,*) timint(i), tstart(i), facsdp(i)
c
c     timint(i): recurrence interval (time) for segment i
c     tstart(i): slip rate of segment i
c     facsdp(i): fraction of aseismic slip of segment i
c
c     This next bit increases all the aseismic creeps
c		to ensure bigger events
c
            if(facsdp(i) .le. .5d0) facsdp(i) = 2.d0* facsdp(i)
c
         end do
         rewind 20

         close(20)
c
c     Read the fault coefficients
c
         imax=nfault*nfault
         imaxm=nfault*nfault*mend
         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,*) (Sxyvec(n), n=1,nfault)
         read(10,*) (Snvec(n), n=1,nfault)
         read(10,*) (Slpvel(n), n=1,nfault)
         read(10,*) (vcovec(n), n=1,imax)
         read(10,*) (vnvec(n), n=1,imax)
c
c     Description of variables: see above
c
	 print *, ' '
	 print *, ' Fault slip velocities in cm/yr are:'
         print *, (Slpvel(n), n=1,nfault)
c
c     Multiply by amu. Also, 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
c     Depth is depth to midpoint of fault.  rhogd is pressure at
c	that depth in bars
c
         do ir = 1,nfault
            rhogd(ir)  =  depth(ir) * 100.d0
         end do
c
	 print *, ' '
	 print *, ' Finished reading Greens function'
	 print *, ' and fault interaction data.'
         print *, ' '
         print *, ' Hit carriage return to continue'
         read(*,'(a1)') resp3
c
         print *,' '
         print *,' Change the NOMINAL TIME INTERVALS by a factor? (y/n)'
         read *, rptime
         if(rptime .eq. 'y') then
            print *, ' '
            print *, ' Enter factor:'
            read *, timfac
            do i=1,nfault
               timint(i) = timint(i)*timfac
               if (timint(i) .gt. 1000.d0) timint(i) = 1000.d0
            end do
         endif

         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)
         end do
c
c     Calculate cfric and dfric
c
         ivar = 583936739
c
c     Static friction coefficient will be 20% random
c
         do i=1,nfault
            dfric(i) = .1d0
            rndfac = (1.d0 + .2d0* (2.d0*(.5d0 - dble(rand()))))
            cfrdif(i) = rndfac * dabs(sdrpL(i)) /
     &           dabs(-rhogd(i) + Snvec(i)*slpvel(i)*timint(i))

            cfric(i) = dfric(i) + cfrdif(i)
	
            txybac(i) = cfric(i) * dabs(rhogd(i))
         end do
c
         vmag = dsqrt(vplatx*vplatx + vplaty*vplaty)
c     
         print *,' '
         tauast = 5.d0
         print *,' NO ASTHENOSPHERIC RELAXATION'
c
c	Last time to slip was 0.d0
c
         do i = 1,nfault
            Tlast(i) = 0.d0
         end do

c
c	This next code overwrites the fault friction file
c		It is disabled for the public access code.
c
c         open(unit=20,name=fnfric,status='old')
c         write(20,'(2g20.10)') amu,pfacl,cfrcl
c         do i=1,nfault
c            write(20,*) timint(i),tstart(i), facsdp(i)
c         end do
c         close(20)

         tmin = 0.d0
c
c     Calculate the slip deficits
c
c
c     Copy final values into matrices
c
         do i=1,nfault
            displt(i) = 0.d0
         end do
c
         do j=1,NTSTP
            do i=1,nfault
               dslip(i,j)=0.d0
               tslip(i,j)=0.d0
            end do
         end do
         do i=1,nfault
            islip(i)=0
         end do
c
         ncycle=1
         ntm=0
c
c     This is the "end if" that goes with rsprs='y':
c
      end if
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

      pdrtcr = 0.d0
      if (ncycle .eq. 1) then
c
c	SET pdecay FOR ANNEALING DURING RUN
c
         pdecy = .0001d0
         pdrate = 0.d0
         rpdec = 'n'
c
         do i=1,nfault
            pdecay(i) = pdecy
         end do
c
c	End if for pdecay-ncycle logic
c
      end if
c
      ptrade = pdecy
      print *, ' '
      print *, ' Set FINAL PARAMETERS for running the'
      print *, ' dynamical code? (y/n)'
      print *, ' '
      print *, ' (MAKE SURE THAT ALL CFFs, WHETHER + OR -'
      print *, ' ARE AT LEAST SAME ORDER OF MAGNITUDE AS'
      print *, ' THE LARGEST NORMALIZED STRESS DROP!!! '
      print *, ' ...which means that all DELSTR should be'
      print *, ' between roughly +/- 1)'
      print *, ' '
      read(*,'(a1)') rsncyc
      ptrade = pdecy

      do ir = 1,nfault
         pfactr(ir) = pdecay(ir)
      end do
c
c	Find max dabs(sdrpL)
c
      sdrmax = -  1.d20
      do ir = 1, nfault
         dsr = dabs(sdrpL(ir))
         if (dsr .gt. sdrpL(ir)) sdrmax = dsr
      end do
c
c	Print values of sdrpL
c
      print *, ' '
      print *, ' NOMINAL STRESS DROPS: ', (sdrpL(n), n = 1,nfault)
      print *, ' Hit carriage return to continue'
      read(*,'(a1)') rpdec

      print *,' '
      print *,' Present Time (years):  ',tpres
      print *,' '
      print *,' '
      print *,' ******************************'
      print *,' '
      print *,' STARTING RUN...'
      print *,' '
      print *,' ******************************'
      print *,' '
      print *,' '
c
      KKOUNT=0
      nflag = 1

      do ir = 1,nfault
         signum(ir)=sgn(slpvel(ir))
         dslrun(ir) = 0.d0
         aslrun(ir) = 0.d0
         dslsm(ir) = disnmS(ir) * tstep / timint(ir)
      end do

      dslmx = - 100000.
      do ir=1,nfault
         if (dabs(disnmL(ir)) .gt. dslmx) then
            dslmx = dabs(disnmL(ir))
         end if
      end do

      do ir = 1,nfault
         dslslp(ir) = disnmL(ir)
         factst = dabs(facsdp(ir) - 1.d0)
         if (factst .lt.  .0001d0) then
            dslslp(ir) = signum(ir) * dslmx
         end if
      end do
c
c	Set value for slip deficit to be 1/pdecay
c
c
c	Set values
c
      do ir = 1,nfault
         slpdf0(ir) = slpdav(ir)
      end do

      if (ncycle .eq. 1) then
         do ir = 1,nfault
            slpdf0(ir) = 0.d0
            slpdef(ir) = slpdf0(ir)
         end do
      end if

      Nsegmx = 0

      do ir = 1,nfault
         dsltot(ir) = 0.d0
         cfrict(ir) = dfric(ir) + cfrdif(ir)
         sdrpLT(ir) =  -  4.d0 * dabs(sdrpL(ir))
         slpdav(ir) = 0.d0
      end do
c
c	!! ***********   Enter Time Loop   *********** !!
c
      do it=1,ntime
c
c     Increment time
c
         t = t + tstep

         KKOUNT=KKOUNT + 1
         print *,' '
c
         do ir=1,nfault
            slpdef(ir)= slpdef(ir) + dsltot(ir)- slpvel(ir) * tstep
            strdfl(ir) = strdif(ir)
         end do
c
         call TSUM(t,tauast,dslip,tslip,tfact,slpvel,slpdef,
     &        timpol,slpdef,nfault,mend,islip)
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 (hypsin)
c
         do ir = 1,nfault
	
            hypsin = signum(ir) * dabs(sdrpL(ir))*
     &           dsinh(pdecay(ir)* (slpdef(ir)-slpdf0(ir)))

            Strnml(ir) = - taun(ir)
            Strshr(ir) = tauxy(ir)* signum(ir) - hypsin + txybac(ir)
         end do
c
c	Now calculate stress difference
c
         do ir = 1,nfault
            strdif(ir) = Strshr(ir) - cfrict(ir)*Strnml(ir)
         end do
c
c	If strdif > 0, readjust txybac
c
         do ir = 1,nfault
            if (ncycle .eq. 1 .and. it .eq. 1) then
               if (strdif(ir) .ge. 0.d0) then
                  txybac(ir) = txybac(ir) - strdif(ir) -
     &                 dble(rand())* 10.d0
               end if

               if (strdif(ir) .le. -10.d0) then
                  txybac(ir) = txybac(ir) - strdif(ir) -
     &                 dble(rand())*10.d0
               end if
            end if

            hypsin = signum(ir) * dabs(sdrpL(ir))*
     &           dsinh(pdecay(ir)* (slpdef(ir)-slpdf0(ir)))

            Strnml(ir) = - taun(ir)
            Strshr(ir) = tauxy(ir)* signum(ir) - hypsin + txybac(ir)
            strdif(ir) = Strshr(ir) - cfrict(ir)*Strnml(ir)
         end do
c
c	Allow precursory slip to occur.  Do not turn failure flag on
c		(iflagp) since this precursory slip is stable.
c
	 do ir =1,nfault	
            Sigfac =  dabs(strdif(ir)/sdrpL(ir))
c
c	Turn off precursory slip until system has settled down
c
            Factor(ir) = 0.d0
c
c	System has settled down.  Also, strdif < 0.
c
            if (strdif(ir) .lt. 0.d0) then
               if (Sigfac .ge. 0.d0 .and.
     &              Sigfac .le. 1.d0) then
                  Factor(ir) = 1.d0 - Sigfac
               end if
            end if

            dsl(ir)= dslsm(ir)* Factor(ir)

            if (dabs(dsl(ir)) .gt. dabs(dslsm(ir)) ) dsl(ir) = dslsm(ir)

            dsl(ir) = 0.d0

            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))
	 end do
c
c	dsl:    increment of slip - here stable
c	displt: total slip over all time steps since begining
c	aslrun: total increment of stable slip in this run
c	dsltot: total increment of unstable slip
C
c	dslrun: total of unstable slip in this run
c	strdrp: variable static stress drop
c
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,slpdef,nfault,mend,islip)
c
         call STRESS(t,tauxy,taun,rhogd,tfact,slpvel,
     &        nfault,mend,slpdef)
c
         do ir = 1,nfault
	
	    hypsin = signum(ir) * dabs(sdrpL(ir))*
     &           dsinh(pdecay(ir)* (slpdef(ir)-slpdf0(ir)))

            Strnml(ir) = - taun(ir)
            Strshr(ir) = tauxy(ir)* signum(ir) - hypsin + txybac(ir)
            strdif(ir) = Strshr(ir) - cfrict(ir)*Strnml(ir)
         end do
c
c     This is the failure loop at time t
c
         do ir=1,nfault
            iflagp(ir)=0
            iflagn(ir)=0
         end do
c
c	iflagp: keeps track of how many total times a site has failed
c
c     Failure loop
c
         do ir = 1,nfault
            slpstr(ir) = slpdef(ir)
         end do

	respc = 'c'

	do while (respc .eq. 'c')
           do ik = 1,nfault
              respc = 's'
              dslint(ik) = 0.d0
              if (strdif(ik) .ge. 0.d0) then
                 respc = 'c'
                 iflagp(ik) = 1
                 rndfac = (1.d0 + .2d0* (2.d0*(.5d0 - dble(rand()))))
                 dslint(ik) = dslslp(ik)* rndfac
              end if
c
c     Goal here is to find the average slip deficit that keeps
c		the system within 1 stress drop of zero
c
              if (strdif(ik) .lt. sdrpLT(ik)) then
                 respc = 'c'
                 iflagn(ik) = 1
                 rndfac = (1.d0 + .2d0* (2.d0*(.5d0 - dble(rand()))))
                 dslint(ik) = - dslslp(ik)*rndfac
              end if

              dsltot(ik) = dsltot(ik) + dslint(ik)
              displt(ik) = displt(ik) + dslint(ik)
              dslrun(ik) = dslrun(ik) + dslint(ik)
              slpstr(ik) = slpstr(ik) + dslint(ik)
           end do

           call TSUM(t,tauast,dslip,tslip,tfact,slpvel,slpstr,
     &  	timpol,slpstr,nfault,mend,islip)
c
           call STRESS(t,tauxy,taun,rhogd,tfact,slpvel,
     &  	nfault,mend,slpstr)
c
           do ir = 1,nfault
	
              hypsin = signum(ir) * dabs(sdrpL(ir))*
     &             dsinh(pdecay(ir)* (slpstr(ir)-slpdf0(ir)))

              Strnml(ir) = - taun(ir)
              Strshr(ir) = tauxy(ir)* signum(ir) - hypsin + txybac(ir)
              strdif(ir) = Strshr(ir) - cfrict(ir)*Strnml(ir)
           end do
c
c	End do for 'do while'
c
        end do

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

        nf_hlf = nfault/2
	strdsm = 0.d0
	cffmx = -1.d20
	cffmn = 1.d20

        do ir = 1,nfault
           islip(ir) = 0
           delstr(ir) = strdif(ir)/sdrmax
           strdsm = delstr(ir)

           if(strdif(ir) .gt. cffmx)  cffmx = strdif(ir)/sdrmax
           if(strdif(ir) .lt. cffmn)  cffmn = strdif(ir)/sdrmax

        end do

        strdsm = strdsm / dfloat(nfault)
c
c     Store the time-status vectors
c
        itmntm=it + ntm
        do jj=1,NF
           isltot(itmntm,jj)=0
        end do
c
c	Calculate number of segments that slipped
c
	Nsegsp = 0
	Nsegsn = 0
	
	do ir = 1,nfault
           if (iflagp(ir) .gt. 0) Nsegsp = Nsegsp + 1
           if (iflagn(ir) .gt. 0) Nsegsn = Nsegsn + 1
	end do

	Nsegsp = Nsegsp + Nsegsn
	
        if (it .gt. 10) then
           if(Nsegsp .gt. Nsegmx) then
              Nsegmx = Nsegsp
              itseg = it
           end if
        end if

        nzero = 0
        do ijn=1,nfault
           Tlslp = t - Tlast(ijn)
           if (Tlslp .gt. 5000.d0) nzero = nzero + 1
        end do
	
	print *, ' '
        print *,'We have completed time step: ', it
        print *,'  ==>  Total Time Steps Requested: ', ntime
        print *,'  ==>  Current Time: ', t
        print *,'  ==>  Number of Segments that Slipped: ', Nsegsp
        print *,'  ==>  Maximum Number of Segments: ', Nsegmx
        print *,'           on Time Step: ', itseg
	print *,'  ==>  Number of NEGATIVE slips: ', Nsegsn
c
        if(KKOUNT .eq. 25) then
           print *,' Fault Coulomb Failure Functions, CFFs (bars) --'
           print *, (strdif(n), n=1,nfault)
           print *, ' '
           print *,'  ==>  Average CFF 25 time steps ago:',strdsl
           print *,'  ==>  Average CFF Now:  ', strdsm
           print *,'  (If Total CFF is growing, repeat run with'
	   print *, '  reduced time step)'
           strdsl = strdsm
           print *, ' '
           print *, ' MAX Norm. CFF: ', cffmx
           print *, ' '
           print *, ' MIN Norm. CFF: ', cffmn
           
           KKOUNT=0
        end if
c
c	Compute average slip deficit
c
        do ir=1,nfault
           slpdav(ir) = slpdav(ir) * dfloat(it-1)
           slpdav(ir) = slpdav(ir) + slpstr(ir)
           slpdav(ir) = slpdav(ir) / dfloat(it)
           txybc(ir) = txybac(ir)
        end do
c
c	Write out output file every .1 * ntime steps
c
	ktime = (it / ktest) * ktest
	
	if (ktime .eq. it) then

           do ir = 1,nfault
	
              hypsin = signum(ir) * dabs(sdrpL(ir))*
     &             dsinh(pdecay(ir)* (slpstr(ir)-slpdf0(ir)))
              
              sinhfc = signum(ir) * dabs(sdrpL(ir))*
     &             dsinh(pdecay(ir)* (slpstr(ir)- slpdav(ir)) )
		
              txybc(ir) = txybac(ir) - hypsin + sinhfc

           end do

           tmin_1 = tmin + dfloat(it)*tstep
           itime = itmntm
c
c     Print out velocities
c
           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
           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
              facsd(i) = sngl(facsdp(i))
c
c	If this is the last time, store facrat = facsdp/timint here
c
              if (rsncyc .eq. 'y') facsd(i) = facsdp(i)/timint(i)
c
              timi(i)=sngl(timint(i))
              taub(i)=sngl(txybc(i))
              rhofcc(i)=sngl(rhogd(i))
              cfr(i)=sngl(cfrdif(i))
              dfr(i)=sngl(dfric(i))
              slpdf(i) = sngl(slpstr(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
c     set nycle = 1 for subsequent runs with real dynamics
c
           if (rsncyc .eq. 'y') ncycle = 1
c
           open(unit=30,name=fnout,status='old')
           
           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     Close files
c
      close(10)
      close(40)
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 *******************************************************
c
c  Definitions--
c
c             tfact is the total time factor for a given fault
c                   and a given slip
c
c ********************************************************
c
c     Enter fault loop
c
      do if=1,nfault
c
c     For each fault, compute all the normalized slip times
c     and then for each m, calculate and sum the time factors
c     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
            do j=1,jend
          
               tbart=(t-tslip(if,j))/(tauast)
          
               if(tbart .le. 10.d0) then

                  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
      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
      function sgn(x)
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