c
c     PROGRAM EQSYN_4_PSLIP.F (earthquake synthesizer)
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     ******************************************************
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
      real*4 pdec
      real*4 rand
      character resp3*1, rptime*1, rpts*1, respc*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,NFT), tslipp(NF,NFT)
      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,NF1), displt(NF), dsltot(NF)
      dimension slpdef(NF)
      dimension Sxyvec(NF), Snvec(NF), Slpvel(NF)
      dimension tauxy(NF), taun(NF)
      dimension tstart(NF), dslip(NF,NFT), tslip(NF,NFT)
      dimension tfact(10,NF), islip(NF), iflag(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), 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)

      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 *,' What name for the existing (input) file?'
      read(*,'(a20)') fnin
      print *,' What 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,facsdp
      read(30,'(f16.8)') 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)
      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
        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

        do i=1,nfault
        timint(i)=dble(timi(i))
        txybac(i)=dble(taub(i))
        
        dfric(i)=dble(dfr(i))
        cfrdif(i) = dble(cfr(i))
        
        cfric(i) = dfric(i) + cfrdif(i)
                
        rhogd(i)=dble(rhofcc(i))
        slpdef(i)=dble(slpdf(i))
        slpvel(i)=dble(slpv(i))
        displt(i)=dble(delts(i))
        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
      if(nthi .lt. NTSTP) then
      read(30,*) ((isltot(i,j), i=1,nthi), j=1,NF1)
      else
      read(30,*) ((isltot(i,j), i=1,ntime), j=1,NF1)
      end if
      read(30,*) (pdec(n), n=1,nfault)
      
        do i=1,nfault
        pdecay(i) = dble(pdec(i))
        end do
        
        close(30)
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
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,NFT
          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,NF1)=isltot(ijnt,NF1)
        end do
c
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     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        Sxym(i)=amu*Sxym(i)*1.d-5
c        Snm(i)=amu*Snm(i)*1.d-5

        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*disnmL(i)
        sdrpL(i)=Sxy(ind)*disnmL(i)
        sdrpS(i) = facsdp*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 *, (Slpvel(n), n=1,nfault)
c
      print *, ' Hit Carriage Return to Begin'
      read(*,'(a1)') resp3
c
c     This is the "else" that goes with rsprs='y':
      else
c
      
c
c     Ask for coefficient file name
c
      print *,'  What is the name of the file containing the'
      print *,' stress 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     Ask for name for the output file
c
      print *,' What 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 *,' What name for the friction file?'
      read(*,'(a20)') fnfric
      
      print *, ' '
      print *, ' Enter PRESLIP factor (0. < f < 1.)'
      
      read *, facsdp
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,'(3g20.10)') timint(i), tstart(i)
        end do
      rewind 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
      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
        end do
        do i=1,nfault
        
c        Sxyvec(i)=amu*Sxyvec(i)*1.d-5
c        Snvec(i)=amu*Snvec(i)*1.d-5

        Sxym(i)=0.d0
        Snm(i)=0.d0

        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 *, 'Hit carriage return to begin'
      read(*,'(a1)') resp3
c
c      print *,' '
c      print *,' The driving coefficients are (bars/yr)'
c      print *,' '
c      print *,' Fault, Shear, Normal Coefficients --'
c        do i=1,nfault
c        write(*,'(i5,2g20.10)') i,Sxyvec(i),Snvec(i)
c        end do


      print *,' '
      print *,' Change the nominal time interval 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
	end do
      endif

      
          do i=1,nfault
          ind=(i-1)*nfault + i
          disnmL(i) = slpvel(i)*timint(i)
          disnmS(i) = facsdp*disnmL(i)
          sdrpL(i)=Sxy(ind)*disnmL(i)
          sdrpS(i) = facsdp * 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(ivar)))))
        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      print *,' What value for asthenosphere relaxation'
c      print *,' time? (years)'
c      read *, tauast
c
c     Start everything off over the first 200 years
c
        do i=1,nfault
        dt = dble(rand(ivar))
        tstart(i) = dt * 200.d0
        end do

  
c
c	Last time to slip was 0.d0
c

	do i = 1,nfault
	Tlast(i) = 0.d0
	end do


      rewind 20
      write(20,'(2g20.10)') amu,pfacl,cfrcl
        do i=1,nfault
        write(20,'(3g20.10)') timint(i),tstart(i)
        end do
        
        tmin = 0.d0
        
c
c     Calculate the slip deficits
c

c     Copy final values into matrices
c
        do i=1,nfault
        id=(i-1)*nfault + i
        slpdef(i) =    slpvel(i)*tstart(i)
c	pdecay(i) = .0001d0
        pdecay(i) = dabs(0.001d0/(slpvel(i)*timint(i)))
c
c	Following is the old style version of pdecay
c


        displt(i) = 0.d0
        end do
c
      do j=1,NFT
        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

c
        print *,' '
        print *,' Present Time (years):  ',tpres
        print *,' '
        print *,' '
        print *,' ******************************'
        print *,' '
        print *,' PREPARE YOURSELF:'
        print *,' '
        print *,' THE ADVENTURE BEGINS...'
        print *,' '
        print *,' ******************************'
        print *,' '
        print *,' '
c
      KKOUNT=0
c      nflag = 10

      nflag = 50 

      
        do ir = 1,nfault
        signum(ir)=sgn(slpvel(ir))
        dslrun(ir) = 0.d0
        aslrun(ir) = 0.d0
        end do
        
        Nsegmx = 0

        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
	  cfrict(ir) = dfric(ir) + cfrdif(ir)
          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 analagous 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
        call STRESS(t,tauxy,taun,rhogd,tfact,slpvel,
     &  	nfault,mend,taufar,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) * dsinh(pdecay(ir)* slpdef(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
	 if (it .eq. 1) then
         
         print *,' INITIAL STRDIF --'
         print *, (strdif(n), n=1,nfault)

         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	
	 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)= disnmS(ir) * Factor(ir)* tstep / timint(ir)
         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
	 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
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,taufar,slpdef)
	 
c
c
c	dsl is increment of unstable slip
c	dsltot is total increment of unstable slip
c	dslrun is total of unstable slip in this run
c	aslrun is total increment of stable slip in this run
c	displt is total slip over all time steps since begining
c	strdrp is variable static stress drop
c	iflag keeps track of how many total times a site has failed
c
c
          do ir=1,nfault
          iflag(ir)=0
          dtauxy(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
c	Figure out if any sites are at failure
c

	  respc = 's'
	  nflagn = 0
	  
            do ir = 1,nfault
c
c	Define stress difference
c
	    hypsin = signum(ir) * dsinh(pdecay(ir)* slpdef(ir))
            Strnml(ir) = - taun(ir)
            Strshr(ir) = tauxy(ir)* signum(ir) - hypsin + txybac(ir)
            strdif(ir) = Strshr(ir) - cfrict(ir)*Strnml(ir)
    	
c
c	Increment iflag if strdif .GT. 0.d0
c
       	    if(strdif(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
     	  if (iflag(ir) .eq. 1 .or. iflag(ir) .eq. 2) inflag = 1
     	  end do

c
c	Only if inflag = 1, i.e., a new failure existed during last 
c		two time steps do we proceed 
c

	if (inflag .eq. 1) then

       	  do ir = 1,nfault
       	  if(iflag(ir) .gt. 0 .and. iflag(ir) .le. nflag) then
       	  respc = 'c'
       	  jp = jp + 1
       	  ifltt(jp) = ir
       	  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 strdif(*)
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) = - strdif(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 
    	  strtst = dtauxy(ikr) * signum(ikr)
    	  strdrp(ikr) = 0.d0
    	 
    	  if(strtst .gt. 0.d0) then
    	  strdrp(ikr) = - dtauxy(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)

	  end do
	  
	end if
 
c
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 + .050d0*  (2.d0*(.5d0 - dble(rand(ivar))) ) )
c     &		* (1.d0 + .010d0*  (2.d0*(.5d0 - dble(rand(ivar))) ) )
     &		* (1.d0 + .25d0*  (2.d0*(.5d0 - dble(rand(ivar))) ) )
            end do
            
c
c	Adjust elastic stresses
c

        call STREL(tauxy,taun,dsl,dtauxy,ifltt,jp,nfault)

     	    do ik = 1,nfault
      	    dsltot(ik) = dsltot(ik) + dsl(ik)
            displt(ik) = displt(ik) + dsl(ik)
            dslrun(ik) = dslrun(ik) + dsl(ik)
     	    end do
      
c
c	End do for "do while"
c
          end do
          
            do ik =1,nfault
            frunsl(ik) = dslrun(ik) / 
     &		(aslrun(ik) + dslrun(ik))
            end do
            
	  nf_hlf = nfault/2

         if(KKOUNT .eq. 25) then
         print *,' DISPLT --'
         print *, (displt(n), n=1,nf_hlf)
         print *,' '
         print *,' FRUNSL(n) --'
         print *, (frunsl(n), n=1,nf_hlf)
	 print *, ' '
         print *,' STRDIF --'
         print *, (strdif(n), n=1,160)
         KKOUNT=0
         end if
         
c
c	Print the fault slip that occurred
c
c	  print *, ' '
c	  print *,' Fault, slip, time--  ', ir, dsl(ir)

	  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
           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,NF1)=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. 1000.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: ', Nsegsl
        print *,'  ==>  Maximum Number of Segments: ', Nsegmx
        print *,'           on Time Step: ', itseg
        print *,'  ==>  Number of Sweeps: ', Nsweep
        print *,'  ==>  Number Locked in Last 1000 years: ',nzero
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
	
c
c     Print out velocities
c
      ntlo=itime + 1
      nthi=ntlo + 1000
      if(nthi .lt. NTSTP) then
        do i=ntlo,nthi
          do j=1,NF1
          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)
c
        do i=1,nfault
        pdec(i)=sngl(pdecay(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))
          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')

      write(30,'(a20)') fnmod
      write(30,'(a20)') fnstr
      write(30,*) nfault,hpl,vplx,vply,taua,tauf,amuu,
     .tminn,tstepp,itime,TBIS,bulkms,ntm,ncycle,facsdp
      write(30,'(f16.8)') 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
      if(nthi .lt. NTSTP) then
      write(30,*) ((isltot(i,j), i=1,nthi), j=1,NF1)
      else
      write(30,*) ((isltot(i,j), i=1,itime), j=1,NF1)
      end if
      write(30,*) (pdec(n), n=1,nfault)
      
      close(30)

c
	
	end if
	
c
c	End do for the time loop
c

        end do
c        
c     Close files
c
      close(10)
      close(20)
      close(40)
c
      stop
      end
c
      subroutine TSUM(t,tauast,dslip,tslip,tfact,slpvel,slpdef,
     &	timpol,slpdfd,nfault,mend,islip)
     
c
c     This subroutine calculates the multiplicative summed
c     time factors for the entire slip history of a given fault
c     
      implicit real*8(a-h,o-z)
      INCLUDE 'EQPARM_4.FOR'
      dimension tfact(10,NF), dslip(NF,NFT), 
     .tslip(NF,NFT), 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             islip is running total of number of times a
c                   fault segment has slipped
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,taufar,slpdef)
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
      subroutine STREL(tauxy,taun,dsl,dtauxy,ifltt,jp,nfault)
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), dtauxy(NF)
      dimension ifltt(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 stresses
c

	do ir = 1,nfault
	dtauxy(ir) = 0.d0
	end do
c
        do ic=1,nfault
        ind=(ic-1)*nfault
          do ir=1,nfault
          id=ind + ir
          tauxy(ir)=tauxy(ir) + dsl(ic)*Sxy(id)
          dtauxy(ir) = dtauxy(ir) + dsl(ic)*Sxy(id)
          taun(ir)=taun(ir) + dsl(ic)*Sn(id)
          end do
        end do

c
      return
      end
c
c      function sgn(x)
c
c      implicit real*8(a-h,o-z)
c      sgn=-1.d0
c      if(x .eq. 0.d0) sgn=0.d0
c      if(x .gt. 0.d0) sgn=1.d0
c      return
c      end







