c
c     Program VC_DEFORM.F
c
c     VC_DEFORM takes the output of the code VC_DEF_GREEN.f, together
c	with the output of the VC_SER.f code, and computes the
c	deformation vectors at the points on a square grid centered
c	on the fault system.
c
c     This code takes the slip histories, and computes surface
c	deformation to be plotted by the IDL-based plot routine
c	VC_VISUALIZE.pro as arrows on a background map of
c	faults, or as a set of wrapped InSAR fringes.
c
c	****** IMPORTANT NOTE ********   !!!!!!!!
c
c     In your file area, there MUST be an existing resident
c		file called:
c
c	  >>>	deform_data.d   <<<
c
c	******************************   !!!!!!!!
c
c
c	Primary modules:
c
c
c	Subroutine TSUM()	-- Calculates the multiplicative summed
c				   time factors for the entire slip
c				   history of a given fault
c
c	Subroutine DEFOR()	-- Computes the displacement at the time t
c
c	Subroutine DGRAF()	-- Computes faults and displacement
c				   arrows due to slip using
c				   kinematic displacement Greens
c				   functions and the slip history
c				   file from VC_SER.f
c
c	Function sgn()      	--  Signum function (x < 0, or x > 0)
c
c
c 	*************************************************************
c
      character fmt1*6, fmt2*13, fmt3*4
      character fname_fault*20
      INCLUDE 'DEPARM.FOR'
      dimension islip(NF)
      dimension dslip(NF,NTSTP), tslip(NF,NTSTP)
      dimension dslpT(NF,NTSTP)
      dimension deltot(NF),timpol(6)
      dimension dfric(NF),rhofac(NF),x(NPT),y(NPT)
      dimension duxp(NPT), duyp(NPT)
      dimension slpvel(NF), dstart(NF), cfric(NF)
      dimension timint(NF), taub(NF), slpdef(NF), timi(NF)
      dimension xe(NF), ye(NF), xw(NF), yw(NF)
      dimension db(NF), dt(NF)
      dimension Ux(NPT),Uy(NPT)
      dimension ptrc(NF)
      dimension Uxcov(NPF),Uycov(NPF)
      dimension Uxmall(NPMF), Uymall(NPMF)
      dimension slpav(NF)
      dimension facsdp(NF)
      dimension isltot(NTSTP,NF)
      character fname1*20, fhname*20
      character ftname*20, fdname*20, fcname*20
      character rmenu*1, resp*1, resp1*1,respfb*1
      common/allsub/hplate,vplatx,vplaty,alam,amu,tauast,
     &     islip,dslip,tslip,dslpT,x,y,
     &     nfault,mend,xe,ye,xw,yw,slpvel,nflt
      common/TimAll/timpol
      common/defr/Ux,Uy,Uxcov,Uycov,Uxmall,Uymall
      common/suball/duxp,duyp,vfact,vplate,tlater,tbase,
     &     tmin,tstep,timmin,ntime,ntm,npoint
c
c ***************************************************************
c
c     Query for slip history file from VC_SER.F
c
      Print *,' '
      print *,' Enter the name of the slip'
      print *,' history file from VC_SER.f:'
      print *, ' '
      read(*,'(a20)') fhname
c
c     Get the name of the fault coefficients
c
      print *,' '
      print *,' Enter the name of the input file'
      print *,' from VC_DEF_GREEN.f with kinematic'
      print *,' Greens functions for observation points?'
      print *, ' '
      read(*,'(a20)') fcname
c
      print *, ' '
      print *, ' !!!*** IMPORTANT NOTE ***!!!'
      print *, ' '
      print *, ' In your file area, there must be a'
      print *, ' pre-existing resident file called:'
      print *, ' '
      print *, '  	deform_data.d  '
      print *, ' '
      print *, '  !!!*********************!!!'
      print *, ' '
c
c     Open the history file and read the stuff in it
c
      open(unit=10,name=fhname,status='old')
c
c     Read in all the restart variables
c
      read(10,'(a20)') fname_fault
      read(10,'(a20)') fnstr
      read(10,*) nfault,hplate,vplatx,vplaty,tauast,taufar,amu,
     .     tmin,tstep,ntime,TBIAS,bulkm,ntm,ncycle,pdcy
      read(10,'(f16.8)') ptrade
      read(10,'(a1)') respfb
c
      read(10,*) (timi(n), n=1,nfault)
      read(10,*) (taub(n), n=1,nfault)
      read(10,*) (deltot(n), n=1,nfault)
      read(10,*) (cfric(n), n=1,nfault)
      read(10,*) (dfric(n), n=1,nfault)
      read(10,*) (slpdef(n), n=1,nfault)
      read(10,*) (slpvel(n), n=1,nfault)
      read(10,*) (rhofac(n), n=1,nfault)
      read(10,*) (islip(n), n=1,nfault)
c
      do i=1,nfault
         read(10,*) (dslip(i,k), k=1,islip(i))
      end do
c
      do i=1,nfault
         read(10,*) (tslip(i,k), k=1,islip(i))
      end do
c
c     State variables for each segment:
c     timi(i): recurrence interval (time) of segment i (corrected)
c     taub(i): background stress on segment i
c     deltot(i): time-independent differential slip deficit part of segment i
c     cfric(i): static friction coefficient of segment i
c     dfric(i): dynamic friction coefficient of segment i
c     slpdef(i): slip deficit of segment i
c     slpvel(i): slip velocity of segment i
c     rhofac(i): pressure in bars at the midpoint of segment i
c     islip(i): # of times segment i has slipped
c     dslip(i,j): amount of slip of segment i when it slipped for the jth time
c     tslip(i,j): time at which segment i slipped for the jth time
c
      nthi=ntime + 1001
c
      if(nthi .lt. NTSTP) then
         read(10,*) ((isltot(i,j), i=1,nthi), j=1,NF)
      else
         read(10,*) ((isltot(i,j), i=1,ntime), j=1,NF)
      end if
c
      read(10,*) (slpav(n), n=1,nfault)
      read(10,*) (facsdp(n), n=1,nfault)
c
c     Something new for this code...this is the INCREMENT in
c	stable + unstable slip on i-th fault segment up to
c	time step j.  This is because of aseismic slip on
c	various faults.
c
      itime = ntime - ntm
c
      read(10,*) ((dslpT(i,j), i=1,nfault), j=1,itime)
c
      close(10)	
c
c     Open the fault geometry file
c
      open(unit=40,name=fname_fault,status='old')

      read(40,*) nfault,mendog
      do i=1,nfault
         read(40,*) db(i),dt(i),xw(i),yw(i),
     .        xe(i),ye(i),slpvel(i),nsegg
      end do
c
      close(40)
c
c     Read the observation point coefficients
c
      open(unit=50,name=fcname,status='old')
      read(50,'(a20)') fdname
      read(50,*) nfault,mend,hplate,vplatx,vplaty,
     .     npoint
      read(50,*) (timpol(n), n=1,6)
      read(50,*) (x(n), n=1,npoint)
      read(50,*) (y(n), n=1,npoint)
      isum = nfault*npoint
      isumm=isum*6
      read(50,*) (Uxcov(n), n=1,isum)
      read(50,*) (Uycov(n), n=1,isum)
      read(50,*) (Uxmall(n), n=1,isumm)
      read(50,*) (Uymall(n), n=1,isumm)
      read(50,*) (Ux(n), n=1,npoint)
      read(50,*) (Uy(n), n=1,npoint)
      close(50)
c
c     x, y: coordinates of observation points
c     Uxcov, Uycov: components of displacements at all observation points of 
c                   all segments
c     Uxmall, Uymall: components of displacements for all
c     Ux, Uy: components of displacements for all observation points
c
      resp='y'
      do while(resp .eq. 'y')
c
         call DGRAF
c
         open(unit=60,name='deform_data.d',status='old')
         write(60,*) vplate, vfact, tbase, tlater, tmin, npoint
         write(60,*) (x(j), j=1,npoint)
         write(60,*) (y(j), j=1,npoint)
         write(60,*) (duxp(j), j=1,npoint)
         write(60,*) (duyp(j), j=1,npoint)
         close(60)
c
         print *, 'Compute another velocity field? (y/n)'
         read(*,'(a1)') resp
c
      end do
c
      end
c
c     SUBROUTINES
c
c
      subroutine TSUM(tfact,deltot)
c
c     INTERFACE:
c     Parameters are passed through the COMMON area
c
c     RETURN VALUE:
c     tfact(i,j): total time factor for the ith term of the Green's function 
c                 expansion for segment j
c     deltot(i): time-independent differential slip deficit part of fault i
c
c     DESCRIPTION:
c     This subroutine calculates the multiplicative summed
c     time factors for the entire slip history of a given fault
c
      INCLUDE 'DEPARM.FOR'
      character resp*1
      dimension tfactb(20,NF), tfactl(20,NF)
      dimension tfact(20,NF), dslip(NF,NTSTP)
      dimension xe(NF),ye(NF),xw(NF),yw(NF),x(NPT),y(NPT)
      dimension deltot(NF),dtotb(NF),dtotl(NF)
      dimension tslip(NF,NTSTP), islip(NF)
      dimension timpol(6), slpvel(NF)
      dimension dslpT(NF,NTSTP)
      dimension duxp(NPT), duyp(NPT)

      common/allsub/hplate,vplatx,vplaty,alam,amu,tauast,
     .     islip,dslip,tslip,dslpT,x,y,
     .     nfault,mend,xe,ye,xw,yw,slpvel,nflt

      common/TimAll/timpol

      common/suball/duxp,duyp,vfact,vplate,tlater,tbase,
     &     tmin,tstep,timmin,ntime,ntm,npoint

c
c *******************************************************
c
c     Enter patch loop
c
      do if=1,nfault
         delltt=tlater-tbase
         deltot(if)=0.
         dtotb(if)=0.
         dtotl(if)=0.
         do m=1,mend
            tfact(m,if)=0.
            tfactb(m,if)=0.
            tfactl(m,if)=0.
         end do
c
         jend=islip(if)
         resp='s'
c
         if(jend .gt. 0) resp='g'
c     
         j=0
         time = timmin - tstep
         do while(resp .eq. 'g')
            j=j+1
            time = time + tstep
c     tbarl=(tlater - tslip(if,j))/tauast
c     tbarb=(tbase - tslip(if,j))/tauast
            tbarl=(tlater - time + .001*tstep)/tauast
            tbarb=(tbase  - time + .001*tstep)/tauast
c
c	This is the OLD code
c
c        if(tbarl .le. 0.) resp='s'
c        if(tbarl .gt. 0.) then
c
c	This is the NEW code
c
            if(tbarl .lt. 0.) resp='s'
            if(tbarl .ge. 0.) then
c     
               dtotl(if)=dtotl(if) + dslpT(if,j)
               
               if(tbarb .gt. 0.)then
                  dtotb(if)=dtotb(if) + dslpT(if,j)
               end if
c
               do m=1,mend
                  tfal=1.
                  tbal=tbarl/timpol(m)
                  if(tbal .le. 10.) then
c
c	*** No viscoelasticity for now***
c
c          tfal=1. - exp(-tbal)
                  end if

                  tfactl(m,if)=tfactl(m,if) + tfal*dslpT(if,j)
                  if(tbarb .gt. 0.) then

                     tfab=1.
                     tbab=tbarb/timpol(m)
                     if(tbab .le. 10.) then
c
c	*** No viscoelasticity for now***
c
c          tfab=1. - exp(-tbab)

                     end if
                     tfactb(m,if)=tfactb(m,if) + tfab*dslpT(if,j)
                  end if
               end do
            end if
         end do
c
         do m=1,mend
c
c     This is the time-dependent differential
c		 slip deficit part
c
            tfact(m,if)=tfactl(m,if)-tfactb(m,if) - slpvel(if)*
     .           (tlater-tbase)
         end do
c
c     This is the time-independent differential
c		 slip deficit part
c
         deltot(if) = dtotl(if) - dtotb(if) - slpvel(if)*
     .        (tlater-tbase)
      end do
c
      return
      end      
c
      subroutine DEFOR(dispx,dispy,nfault,mend)
c
c     INTERFACE:
c     dispx, dispy: x and y components of the displacement of a point
c                   of the grid
c     nfault: number of fault segments
c     mend: number of terms (images) og the Green's function expansion
c     other input parameters passed through COMMON statement
c
c     RETURN VALUE:
c     dispx, dispy for all points on the grid
c
c     DESCRIPTION:
c     This subroutine computes the displacement at the time t
c
      INCLUDE 'DEPARM.FOR'
      dimension tfact(20,NF)
      dimension Dispx(NPT),Dispy(NPT),deltot(NF)
      dimension Ux(NPT),Uy(NPT)
      dimension Uxcov(NPF),Uycov(NPF)
      dimension Uxmall(NPMF), Uymall(NPMF)
      dimension duxp(NPT), duyp(NPT)

      common/defr/Ux,Uy,Uxcov,Uycov,Uxmall,Uymall

      common/suball/duxp,duyp,vfact,vplate,tlater,tbase,
     &     tmin,tstep,timmin,ntime,ntm,npoint
c
c     Compute time factors
c
      call TSUM(tfact,deltot)
c
      do i=1,npoint
         Dispx(i)=0.
         Dispy(i)=0.
      end do
c
      l=0
      k=0
      do if=1,nfault
c
         do n=1,npoint
            l=l + 1
            Dispx(n) = Dispx(n) + Uxcov(l)*deltot(if)
            Dispy(n) = Dispy(n) + Uycov(l)*deltot(if)
         end do
c
c	*** No viscoelasticity ***
c
c          do n=1,npoint
c            do m=1,mend
c            k=k + 1
c            Dispx(n)=Dispx(n) + Uxmall(k)*tfact(m,if)
c            Dispy(n)=Dispy(n) + Uymall(k)*tfact(m,if)
c            end do
c          end do
c
      end do
c
c     Add the steady state component
c
      do ip=1,npoint
         Dispx(ip)=Dispx(ip) + Ux(ip)*(tlater-tbase)
         Dispy(ip)=Dispy(ip) + Uy(ip)*(tlater-tbase)
      end do
c
      return
      end
c
      subroutine DGRAF
c
c     DESCRIPTION:
c
c     This subroutine plots faults and displacement arrows due to slip
c     on the major fault systems in California, on a map of the faults
c     in Southern California.
c
c     ******************************************************************
c
      INCLUDE 'DEPARM.FOR'
      dimension x(NPT),y(NPT),dux(NPT),duy(NPT),xp(NPT),yp(NPT)
      dimension duxp(NPT), duyp(NPT)
      dimension xw(NF),xe(NF),yw(NF),ye(NF), slpvel(NF)
      dimension dslip(NF,NTSTP),tslip(NF,NTSTP),islip(NF)
      dimension xl(NPT), yl(NPT)
      dimension dslpT(NF,NTSTP)
      character model*80,MDL*80
      character fname1*20,fname2*20,resp*1
      character*1 respns
      common/allsub/hplate,vplatx,vplaty,alam,amu,tauast,
     .     islip,dslip,tslip,dslpT,x,y,nfault,mend,
     .     xe,ye,xw,yw,slpvel,nflt
      common/suball/duxp,duyp,vfact,vplate,tlater,tbase,
     &     tmin,tstep,timmin,ntime,ntm,npoint
c
c     ******************************************************************
c
c     Get the times
c
      timmin =  tmin - float(ntime - ntm)*tstep
c
      print *, ' '
      print *, ' Minimum time (yr):  ', timmin
	t_max = timmin + float(ntime-ntm)*tstep
      print *, ' '
      print *, ' Maximum time (yr):  ', t_max
      print *,' '
      print *,' We need two dates, a later one'
      print *,' and an earlier one to calculate'
      print *,' the differential offsets.'
      print *,' '
      print *,' Enter earlier date (years)'
      read *, tbase
      print *,' '
      print *,' Enter later date (years)'
      read *, tlater
c
      vplate = sqrt(vplatx*vplatx + vplaty*vplaty)
      vpl = vplate
      vfact = vplate*(tlater-tbase)
c
c     Calculate the displacements
c
      call DEFOR(dux,duy,nfault,mend)
c
      dxmax = 0.
      dymax = 0.
      do i=1,nfault
         dxtest = sqrt(dux(i)*dux(i))
         dytest = sqrt(duy(i)*duy(i))
         if(dxtest .gt. dxmax) dxmax = dxtest
         if(dytest .gt. dymax) dymax = dytest	
      end do
c
      do i = 1, npoint
         xl(i) = x(i)/100.
         yl(i) = y(i)/100.
c	*******************************
c
c	Old definition of xp, yp
c
c        xp(i)=dux(i)/vfact + xl(i)
c        yp(i)=duy(i)/vfact + yl(i)
c
c	*******************************
c
c	New defintion of duxp,duyp
c
         duxp(i) = dux(i)/vfact
         duyp(i) = duy(i)/vfact
c
      end do
c
      return
      end
c
c     Function Sgn
c
      function sgn(x)
c
c     This is the signum function
c
      sgn=1.
      if(x .lt. 0.) sgn=-1.
      if(x .eq. 0.) sgn=0.
      return
      end