c
c	PROGRAM VC_DEF_GREEN.F
c
c	Modified July 16, 2002
c
c     	This fortran program takes the fault topology data
c     	from the data file, makes the necessary coordinate
c	transformations, and then calculates the deformation
c	(kinematic) Greens functions for a square grid of
c	surface observation points.  Only Greens functions
c	for horizontal surface deformation are computed.
c
c	The primary input to this code is the fault topology
c	data file, an example is VC_FAULTS_1999.d as given
c	in the accompanying data files.
c
c	For details about parameter definitions and other
c	information, see the header, or preamble to the code
c	VC_SER.f, the Virtual California earthquake
c	simulation code.
c
c	The output file from this code is used, together
c	with an output file from VC_SER.f, as input to
c	the code VC_DEFORM.f  That code then computes
c	the actual horizontal surface deformation difference
c	between any two times selected by the user, as
c	input for the visualization and display code
c	VC_VISUALIZE.pro With the visualization code, surface
c	deformation can be plotted either as GPS-type arrows,
c	or as wrapped InSAR fringes.
c
c	Program Modules:
c
c	Subroutine CALTRM()		-- Calculates the m-dependent
c					   terms for the Rundle-Jackson
c					   kinematic Greens function
c
c	Subroutine STEADY()		-- Calculates the steady state
c					   part of the deformation, the
c					   part that is due to plate
c					   tectonic motion in the far
c				   	   field
c	
c	Subroutine HLFSTR()		-- Evaluates the displacements
c					   parallel and perpendicular
c					   to a strike slip fault in a
c					   homogeneous, elastic half space.
c
c	Subroutine TREVNI()		-- Computes the slip variables for
c					   each fault segment
c
c	Function sgn()	    --	Signum function (x < 0, or x > 0)
c
c	Also, many assorted LINPACK subroutines that are called from
c		Subroutine TREVNI
c
c
c     ***************************************************
c
      integer recrd
      character*80 model
      character fm7*23,fm8*15,fm9*11
      character fma80*7
      character resp1*1, resp2*1, resp3*1, resp4*1
      character fmt1*8,fmt2*15,fmt3*6,fmt4*14,fma*7
      character fname*20, fname5*20
      INCLUDE 'DEPARM.FOR'
      dimension xpt(NPT),ypt(NPT),zpt(NPT)
      dimension timpol(10), amatm(6,6), bmatm(6,6), Uxvec(6),
     .     Uyvec(6), ansx(6), ansy(6)
      dimension xptt(20,NPT),yptt(20,NPT)
      dimension x(NPT),y(NPT),alen(NF),xxpp(NPT),yypp(NPT)
      dimension dbp(NF),dtp(NF),xpwn(NF),ypwn(NF),slpvel(NF),
     .     xpen(NF),ypen(NF),unitx(NF),unity(NF)
      dimension xfwn(NF),yfwn(NF),xfen(NF),yfen(NF)
      dimension dbn(NF),dtn(NF)
      dimension thetan(NF),xmidn(NF),ymidn(NF)
      dimension Ux(NPT),Uy(NPT),Uxm(20,NPT),Uym(20,NPT)
      dimension Uxmp(10,NPT), Uymp(10,NPT)
      dimension unx(NF),uny(NF),Uxcov(NPF),Uycov(NPF)
      dimension Uxmall(NPMF), Uymall(NPMF)
      dimension AMAT(NF,NF),IPVT(NF),WORK(NF), ZVEC(NF)
      dimension upar(NPT),uper(NPT)
      dimension uparp(NPT,20),uperp(NPT,20)
      
      dimension Uxcov2(NPF), Uycov2(NPF)
      dimension Uxmal2(NPMF), Uymal2(NPMF)
c
c     Get the fault data file name
c
      print *, ' '
      print *,' Enter the name of the fault topology data file:'
      read(*,'(a20)') fname
c
c     Get the number of faults and number of terms in the Green function
c     expansion
c
      open(unit=20,name=fname,status='old')
      fname5=fname
      read(20,*) nfault,mend
c
c     nfault: number of fault segments
c     mend: number of terms (images) in Green's function expansion
c
      do i=1,nfault
         read(20,*) dbn(i),dtn(i),xfwn(i),yfwn(i),
     .        xfen(i),yfen(i),slpvel(i),nsegg
      end do
c
c     dbn, dtn: depth of bottom and top of segment
c     xfwn, yfwn: x and y coordinates of west end of segment
c     xfen, yfen: ...of east end
c
c     Change sign of y so that +=south
c
      do i=1,nfault
         yfen(i)=-yfen(i)
         yfwn(i)=-yfwn(i)
      end do
        
      read(20,*) vplatx, vplaty
c
c     vplatx, vplaty: x and y components of plate velocity
c
c     Change sign of vplate y to coincide with south facing y-axis
c
      vplaty=-vplaty
c
      pole1=.5
      pmult=2.
      timpol(1)=pole1
      do i=2,6
         timpol(i)=timpol(i-1)*pmult
      end do
c
      do i=1,6
         si=1./timpol(i)
         do j=1,6
            sj=1./timpol(j)
            amatm(i,j)=sj/(si*(si+sj))
            bmatm(i,j)=amatm(i,j)
         end do
      end do
c
c     Read the slip history file
c
      print *, ' '
      print *,' Enter the name of the slip history file:'
      print *,' (output from VC_SER.f)'
      read(*,'(a20)') fname
c
      open(unit=50,name=fname,status='old')
      
      read(50,'(a20)') fname
      read(50,'(a20)') fname
      read(50,*) nfault,hplate,vplx,vply,taua,tauf,amuu,
     .     tminn,tstepp,ntime,TBIS,bulkm,ntm,ncycle,pdcy
      close(50)
c
c     Open the output file
c
      print *,' '
      print *,' Enter the name for the output file:'
      read(*,'(a20)') fname
      open(unit=90,name=fname,status='new')
c
      indn=0
      indnm=0
      do i=1,nfault
c
c     Define variables
c
         xfe=xfen(i)
         yfe=yfen(i)
         xfw=xfwn(i)
         yfw=yfwn(i)
         dt=dtn(i)
         db=dbn(i)
c
c     Compute the fault midpoint
c
         xmidn(i)=.5*(xfw+xfe)
         ymidn(i)=.5*(yfw+yfe)
c
c     Compute the rotation angle
c
         testy=yfw-yfe
         if(abs(testy) .le. .01) then
            xmax=amax1(xfw,xfe)
            xmin=amin1(xfw,xfe)
            xfw=xmax
            xfe=xmin
            thetan(i)=0.
c
         else
            delx=xfe-xmidn(i)
            dely=yfe-ymidn(i)
            if(delx .le. 1.) delx=(1.e-10)*sgn(delx)
            thetan(i)=atan(dely/delx)
         end if
c
      end do
c
c	Find center point of fault grid to add to each point
c
      xmax = -10000.
      xmin =  10000.
      ymax = -10000.
      ymin =  10000.
c	
      do i = 1,nfault
         xmax1 = amax1(xfen(i),xfwn(i))
         if (xmax1 .gt. xmax) xmax = xmax1
         xmin1 = amin1(xfen(i),xfwn(i))
         if (xmin1 .lt. xmin) xmin = xmin1
         ymax1 = amax1(yfen(i),yfwn(i))
         if (ymax1 .gt. ymax) ymax = ymax1
         ymin1 = amin1(yfen(i),yfwn(i))
         if (ymin1 .lt. ymin) ymin = ymin1
      end do
c	
      xmidf = .5 *(xmax + xmin)
      ymidf = .5 *(ymax + ymin)
c
c     Define grid of observation points
c
      npoint = 400
c
c     Define points with respect to center of fault grid
c	then add mid point to the coordinates
c
      l = 0
      do i = 1,20
         do j = 1,20
            l = l+1
            xpt(l) = -400. + 40.*float(i-1) + xmidf
            ypt(l) = -400. + 40.*float(j-1) + ymidf
         end do
      end do
c
c     Rotate observation points
c
      write(90,'(a20)') fname5
c
c     Turn vplaty around 
c
      ntterm = 6
      vpy = -vplaty
      write(90,*) nfault,ntterm,hplate,vplatx,vpy,
     .     npoint
      write(90,*) (timpol(n), n=1,6)
      write(90,*) (xpt(n), n=1,npoint)
c
c     Write y's with y axis north
c
      do iy=1,npoint
         y(iy) = -ypt(iy)
      end do
c
      write(90,*) (y(n), n=1,npoint)
c
      print *, ' '
      htest=3.*hplate
      do i=1,nfault
         st=sin(thetan(i))
         ct=cos(thetan(i))
         xmid=xmidn(i)
         ymid=ymidn(i)
c
c     All points are affected by all faults
c   
         do n=1,npoint
            xpt(n)=xpt(n)-xmid
            ypt(n)=ypt(n)-ymid
            x(n)= xpt(n)*ct + ypt(n)*st
            y(n)=-xpt(n)*st + ypt(n)*ct
            xxpp(n)=x(n)
            yypp(n)=y(n)
            xpt(n)=xpt(n)+xmid
            ypt(n)=ypt(n)+ymid
         end do
c
         sl=.5*sqrt((xfwn(i)-xfen(i))**2+(yfwn(i)-yfen(i))**2)
         dt=dtn(i)
         db=dbn(i)
c
c     Calculate displacements
c
         amu = 1.
         alam = bulkm - 2./3.
c
c     Change sign so that right lateral slip is positive
c
         uds=-1.
c
         print *, ' Doing Calculations for Fault Segment: ', i
c
         call CALTRM(upar,uparp,uper,uperp,xxpp,yypp,zpt,
     &        hplate,uds,sl,db,dt,alam,amu,npoint,mend)

c
c     Rotate the displacements back into the correct 
c     reference frame
c
         do np=1,npoint
            Ux(np)= upar(np)*ct - uper(np)*st
            Uy(np)= upar(np)*st + uper(np)*ct
         end do
c
         do np=1,npoint
            do m=1,mend
               Uxm(m,np)= uparp(np,m)*ct - uperp(np,m)*st
               Uym(m,np)= uparp(np,m)*st + uperp(np,m)*ct
c	******************************************
c	Set Uxm, Uym = 0 since for now we dont use
c		viscoelasticity
c	******************************************
               Uxm(m,np)  =  0.
               Uym(m,np)  =  0.
            end do
         end do
c
c     Adjust the m-dependent terms
c
         do ip=1,npoint
            mindx=mend
            if(mindx .gt. 10) mindx=10
            do ii=1,6
               si=1./timpol(ii)
               fact=1.
               Uxvec(ii)=0.
               Uyvec(ii)=0.
               do m=1,mindx
                  fact=(1.d0/(1.d0+si))*fact
                  if(abs(fact) .le. 1.e-15) fact=0.d0
                  Uxvec(ii)=fact*Uxm(m,ip) + Uxvec(ii)
                  Uyvec(ii)=fact*Uym(m,ip) + Uyvec(ii)
               end do
               Uxvec(ii)=Uxvec(ii)/si
               Uyvec(ii)=Uyvec(ii)/si 
            end do
c	************************************
c	Comment out next statements since no
c		viscoelasticity
c
c          itask=1
c          call TREVNI(amatm,Uxvec,ansx,itask)
c          call TREVNI(bmatm,Uyvec,ansy,itask)
c	************************************
            do m=1,6
               Uxmp(m,ip)=ansx(m)
               Uymp(m,ip)=ansy(m)
c	************************************
c	Uxmp, Uymp = 0 since no viscoelasticity
c	************************************
               Uxmp(m,ip)=0.
               Uymp(m,ip)=0.
            end do
         end do
c
c     Remember to turn the y-axis around again
c
         do np=1,npoint
            indn=indn + 1
            Uxcov(indn) = Ux(np)
            Uycov(indn) = -Uy(np)
         end do
c
         do np=1,npoint
            do m=1,6
               indnm=indnm + 1
               Uxmall(indnm) = Uxmp(m,np)
               Uymall(indnm) = -Uymp(m,np)
            end do
         end do
c
c      This * end do * goes with the loop over i above
c
      end do
c
c     Write the output into file 90
c
      isum = nfault*npoint
      isumm=isum*6
c         
      do ik = 1,isum
         Uxcov2(ik) = 0.
         Uycov2(ik) = 0.
         if(abs(Uxcov(ik)) .ge. 1.e-10) Uxcov2(ik) = Uxcov(ik)
         if(abs(Uycov(ik)) .ge. 1.e-10) Uycov2(ik) = Uycov(ik)
      end do
	
      do ik = 1,isumm
         Uxmal2(ik) = 0.
         Uymal2(ik) = 0.
         if(abs(Uxmall(ik)) .ge. 1.e-10) Uxmal2(ik) = Uxmall(ik)
         if(abs(Uymall(ik)) .ge. 1.e-10) Uymal2(ik) = Uymall(ik)
      end do

      write(90,*) (Uxcov2(n), n=1,isum)
      write(90,*) (Uycov2(n), n=1,isum)
      write(90,*) (Uxmal2(n), n=1,isumm)
      write(90,*) (Uymal2(n), n=1,isumm)
c
c     ************************************************
c
c     Calculate the steady state part
c
c     ************************************************
c
c     Find northernmost end
c
      xnorth=0.
      ynorth=0.
      do i=1,nfault
         if(yfwn(i) .le. ynorth) then
            inorth=i
            ynorth=yfwn(i)
            xnorth=xfwn(i)
         end if
      end do
c
c     Find southernmost end
c
      xsouth=xnorth
      ysouth=ynorth
      do i=1,nfault
         if(yfen(i) .ge. ysouth) then
            isouth=i
            ysouth=yfen(i)
            xsouth=xfen(i)
         end if
      end do
c
c     Add two new faults
c
      vmag=sqrt(vplatx**2 + vplaty**2)
      inew=nfault + 1
      yfen(inew)=ynorth
      xfen(inew)=xnorth
      yfwn(inew)=-abs(vplaty*3000./vmag) + ynorth
      xfwn(inew)=-abs(vplatx*3000./vmag) + xnorth
c     
      slpvel(inew)=vmag
c
      inew=nfault + 2
      yfwn(inew)=ysouth
      xfwn(inew)=xsouth
      yfen(inew)=abs(vplaty*3000./vmag) + ysouth
      xfen(inew)=abs(vplatx*3000./vmag) + xsouth
c     
      slpvel(inew)=vmag
c      
      nflt=nfault + 2
c
      do i=1,NPT
         Ux(i)=0.
         Uy(i)=0.
      end do
c
c     Define variables
c
      do i=1,nflt
         xfe=xfen(i)
         yfe=yfen(i)
         xfw=xfwn(i)
         yfw=yfwn(i)
c
c     Compute the fault midpoints
c
         xmidn(i)=.5*(xfw+xfe)
         ymidn(i)=.5*(yfw+yfe)
c
c     Compute the rotation angle
c
         testy=yfw-yfe
         if(abs(testy) .le. 1.) then
            xmax=amax1(xfw,xfe)
            xmin=amin1(xfw,xfe)
            xfw=xmax
            xfe=xmin
            thetan(i)=0.
c
         else 
            delx=xfe-xmidn(i)
            dely=yfe-ymidn(i)
            if(delx .le. 1.) delx=(1.e-10)*sgn(delx)
            thetan(i)=atan(dely/delx)
         end if
      end do
c
c     Rotate observation points
c
      do i=1,nflt
         st=sin(thetan(i))
         ct=cos(thetan(i))
         xmid=xmidn(i)
         ymid=ymidn(i)
c
         do n=1,npoint
            xpt(n)=xpt(n)-xmid
            ypt(n)=ypt(n)-ymid
            x(n)= xpt(n)*ct + ypt(n)*st
            y(n)=-xpt(n)*st + ypt(n)*ct
            xpt(n)=xpt(n)+xmid
            ypt(n)=ypt(n)+ymid
         end do
c
         sl=.5*sqrt((xfwn(i)-xfen(i))**2+(yfwn(i)-yfen(i))**2)
         dt=0.
         if (i .eq. nflt) dt = 20.
         db=3000.
c
c     Calculate displacements
c
c	********************************
c	amu = 1. since no viscoelasticity
c	********************************
c          amu=.001
         amu = 1.
         alam = 1.
c     
c	Turn uds off except for north
c
         uds = slpvel(i)
c
c     Change sign so that right lateral slip is positive
c
         uds=-1.*uds
c
         call STEADY(upar,uper,x,y,zpt,uds,sl,db,dt,alam,
     .        amu,npoint)
c
c     Rotate the displacements back into the correct reference frame.  
c
         do n=1,npoint
            Ux(n)= Ux(n) + upar(n)*ct - uper(n)*st
            Uy(n)= Uy(n) + upar(n)*st + uper(n)*ct
         end do
c
c      This * end do * goes with the loop over i above
c
      end do
c
c     Remember to turn the y-axis around
c
      do n=1,npoint
         Uy(n) = -Uy(n)
      end do
c
c     Write the output into file 90
c
      write(90,*) (Ux(n), n=1,npoint)
      write(90,*) (Uy(n), n=1,npoint)
c
      close(20)
      close(30)
      close(40)
      close(90)
      stop
      end
c
      function sgn(x)
c
c     DESCRIPTION:
c     This function computes the sign of the supplied 
c     quantity x
c
      sgn=1.
      if(x .lt. 0.) sgn=-1.
      if(x .eq. 0.) sgn=0.
      return
      end
c
c     ******************************************************
c
      subroutine CALTRM(upa,uparp,upe,uperp,
     .xpt,ypt,zpt,h,uds,sl,db,dt,alam,amu,npoint,mend)
c
c     INTERFACE:
c     upa, uparp: parallel displacements (m)
c     upe, uperp: perpendicular displacements (m)
c     xpt,ypt,zpt: coordinates of the observation
c               point in chinnery's coordinates (km)
c     h: layer thickness (km)
c     uds: fault slip (m)
c     sl: the fault semilength (km)
c     db: depth to the fault bottom (km)
c     dt: depth to the fault top (km)
c     alam: first lame constant (any unit)
c     amu: shear modulus (same unit as alam)
c     npoint:the number of observation points
c     mend: number of terms in the Green's function expansion
c           (kept in the summation)
c
c     RETURN VALUE:
c     uparp, uperp
c
c     DESCRIPTION:
c     This subroutine calculates the m-dependent
c     terms for the rundle-jackson green function
c          
c     ******************************************************
c
      INCLUDE 'DEPARM.FOR'
      dimension xpt(NPT),ypt(NPT),zpt(NPT)
      dimension upa(NPT),upe(NPT)
      dimension uparp(NPT,20),uperp(NPT,20)
      dimension zmh(NPT),uv(NPT),upap(NPT),upep(NPT)
c
c     *******************************************************
c
      if(npoint .gt. NPT) pause 'number of observation 
     .     points too big'
c
      do ipoint=1,NPT
         upe(ipoint)=0.
         upa(ipoint)=0.
         upap(ipoint)=0.
         upep(ipoint)=0.
      end do
c  
      call HLFSTR(db,dt,sl,xpt,ypt,zpt,uds,upa,upe,
     .     uv,alam,amu,npoint)
c
      do m=1,mend
         do ipoint=1,npoint
            zmh(ipoint)=2.*float(m)*h 
         end do
c
         call HLFSTR(db,dt,sl,xpt,ypt,zmh,uds,upap,upep,
     .        uv,alam,amu,npoint)
c     
c     We don't call hlfstr to compute upam, upem, because
c     the answer is the same as for upap, upep for 
c     surface observation points 
c
c     Define the m-dependent terms
c
         do ipoint=1,npoint
            uparp(ipoint,m)=upap(ipoint)
            uperp(ipoint,m)=upep(ipoint)
         end do
      end do
      return
      end
c
      subroutine STEADY(upa,upe,xpt,ypt,zpt,
     .     uds,sl,db,dt,alam,amu,npoint)
c
c     INTERFACE:
c     upa: parallel displacements (m)
c     upe: perpendicular displacements (m)
c     xpt,ypt,zpt: coordinates of the observation
c               point in chinnery's coordinates (km)
c     uds: fault slip (m)
c     sl: the fault semilength (km)
c     db: depth to the fault bottom (km)
c     dt: depth to the fault top (km)
c     alam: first lame constant (any unit)
c     amu: shear modulus (same unit as alam)
c     npoint:the number of observation points
c
c     RETURN VALUE:
c     upa, upe
c
c     DESCRIPTION:
c     This subroutine calculates the steady state
c     part of the solution
c          
c     ******************************************************
c
      INCLUDE 'DEPARM.FOR'
      dimension xpt(NPT),ypt(NPT),zpt(NPT)
      dimension upa(NPT),upe(NPT),uv(NPT)
c
c     *******************************************************
c
      do ipoint=1,npoint
         upe(ipoint)=0.
         upa(ipoint)=0.
      end do
c  
      do ipoint=1,npoint
         zpt(ipoint)=0.
      end do
c
c     Call the half space fault calculations
c
      call HLFSTR(db,dt,sl,xpt,ypt,zpt,uds,upa,upe,
     .     uv,alam,amu,npoint)
c
      return
      end
c
      subroutine HLFSTR(db,dt,sl,xpt,ypt,zpt,uds,upara,
     .     uperp,uvert,alam,amu,npoint)
c
c     INTERFACE:
c     db: depth to the fault bottom (km)
c     dt: depth to the fault top (km)
c     sl: the fault semilength (km)
c     xpt,ypt,zpt: coordinates of the observation
c               point in chinnery's coordinates (km)
c     uds: fault slip (m)
c     upara: displacement parallel to strike (m)
c     uperp: displacement perpendicular to strike (m)
c     uvert: displacement in vertical direction (m)
c     alam: first lame constant (any unit)
c     amu: shear modulus (same unit as alam)
c     npoint:the number of observation points
c
c     RETURN VALUE:
c     upara, uperp
c
c     DESCRIPTION:
c     This subroutine evaluates the displacements parallel and 
c     perpendicular to a strike slip fault in a homogeneous, 
c     elastic half space.
c    
c     ***************************************************
      implicit real*8(a-h,o-z)
      real*4 uperp,upara,uvert,xpt,ypt,zpt,db,dt,sl,uds,alam,amu
      INCLUDE 'DEPARM.FOR'
      dimension upp(NPT),upa(NPT),uvt(NPT)
      dimension uperp(NPT),upara(NPT),uvert(NPT),
     &	xpt(NPT),ypt(NPT),zpt(NPT)
c
c     ********************************************************
c
      fourpi=12.5663706140d0
      u=dble(uds)/fourpi
      al=dble(sl)
      ddb=dble(db)
      ddt=dble(dt)
      dalam=dble(alam)
      damu=dble(amu)
      alph=(dalam+damu)/(dalam+2.d0*damu)
      beta=(dalam-damu)/(dalam+damu)
      cee=damu*dalam/((dalam+damu)**2)
c
      do jk=1,npoint
         x=dble(xpt(jk))
         y=dble(ypt(jk))
         z=dble(zpt(jk))
c     
         upp(jk)=0.d0
         upa(jk)=0.d0
         uvt(jk)=0.d0
c
c     Define parameters
c
         do j=1,4
            if(j .eq. 1) then
               x1=al
               x3=ddb
               sign=1.d0
            else if(j .eq. 2) then
               x1=-al
               x3=ddb
               sign=-1.d0
            else if(j .eq. 3) then    
               x1=al
               x3=ddt
               sign=-1.d0
            else if(j .eq. 4) then
               x1=-al
               x3=ddt
               sign=1.d0
            end if
c
            p=x3+z
            q=x3-z
            t=x1-x
            s1=dsqrt(t*t+y*y+q*q)
            s2=dsqrt(t*t+y*y+p*p)
            s2sq=t*t+y*y+p*p
c
c     Calculate displacements
c
            upa(jk)=upa(jk) + (y*t*alph*(1.d0/(s1*(s1+q))
     .           + ((1.d0+cee)*s2+(1.d0-beta)*p+q)/(s2*(s2+p)*(s2+p))
     .           - ((p*p-q*q)*(2.d0*s2+p))/(2.d0*s2*s2sq*(s2+p)*(s2+p))
     .           )+ datan(q*t/(y*s1)) + datan(p*t/(y*s2)))*sign
c     
            upp(jk)=upp(jk) + ((1.d0-alph)*dlog(s1+q)+
     .           (1.d0-alph-alph*cee)*dlog(s2+p) -
     .           alph*(((beta+cee)*p-q)/(s2+p) + (p*p-q*q)/(2.d0*s2*
     .           (s2+p))) - y*y*alph*(1.d0/(s1*(s1+q))
     .           +((1.d0+cee)*s2+(1.d0-beta)*p+q)/(s2*(s2+p)**2)
     .           -((p*p-q*q)*(2.d0*s2+p))/(2.d0*s2*s2sq*(s2+p)**2)))
     .           *sign
c     
            uvt(jk)=uvt(jk) + (y*(alph*(1.d0/s1-1.d0/s2) +
     .           (s2*(1.d0-beta)+2.d0*alph*p + 2.d0*q*(1.d0-alph))/
     .           (s2*(s2+p)) + (alph*(p*p-q*q))/(2.d0*s2*s2sq)))
     .           *sign
c     
         end do
c
         upa(jk)=upa(jk)*u
         upp(jk)=upp(jk)*u
         uvt(jk)=uvt(jk)*u
         upara(jk)=sngl(upa(jk))
         uperp(jk)=sngl(upp(jk))
         uvert(jk)=sngl(uvt(jk))
      end do
c
      return
      end
c

      subroutine TREVNI(Amats,sdrops,deltas,itask)
c
c     This subroutine computes the slip variables for each 
c     fault segment
c
      implicit real*8(a-h,o-z)
      real*4 Amats,sdrops,deltas
      
      dimension Amats(6,6),sdrops(6),deltas(6)
      
      dimension Amat(6,6),sdrop(6),delta(6)
      dimension work(6),iwork(6)
c
c     Invert the system of equations
c
      ijob = 0
	
      do i = 1,6
         do j = 1,6
            Amat(i,j) = dble(Amats(i,j))
         end do
      end do
      
      do i = 1,6
         sdrop(i) = dble(sdrops(i))
      end do

      call dgeco(Amat,6,6,iwork,rcond,work)
      call dgesl(Amat,6,6,iwork,sdrop,ijob)
c
      if(ind .lt. 0) then
         print *,' ** IND in DGEFS .lt. 1 **'
         print *,' IND =', IND
      end if
c
      do i=1,6
         deltas(i)=sdrop(i)
      end do
c
      return
      end
c
      double precision function dasum(n,dx,incx)
c ------------------------------------------------------------------c
c     takes the sum of the absolute values.
c     jack dongarra, linpack, 3/11/78.
c ------------------------------------------------------------------c
      double precision dx(1),dtemp
      integer i,incx,m,mp1,n,nincx
 
      dasum = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c ---------------------------------
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
         dtemp = dtemp + dabs(dx(i))
 10   continue
      dasum = dtemp
      return
c
c code for increment equal to 1
c
c clean-up loop
c -----------------------------
c
 20   m = mod(n,6)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
         dtemp = dtemp + dabs(dx(i))
 30   continue
      if( n .lt. 6 ) go to 60
 40   mp1 = m + 1
      do 50 i = mp1,n,6
         dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2))
     *        + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5))
 50   continue
 60   dasum = dtemp      
      return
      end
c
      subroutine daxpy(n,da,dx,incx,dy,incy)
c---------------------------------------------------------------c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c---------------------------------------------------------------c
      double precision dx(1),dy(1),da
      integer i,incx,incy,m,mp1,n
 
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c -----------------------------------------------
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
         dy(iy) = dy(iy) + da*dx(ix)
         ix = ix + incx
         iy = iy + incy
 10   continue
      return
c
c code for both increments equal to 1
c
c clean-up loop
c -----------------------------------
c
 20   m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
         dy(i) = dy(i) + da*dx(i)
 30   continue
      if( n .lt. 4 ) return
 40   mp1 = m + 1
      do 50 i = mp1,n,4
         dy(i) = dy(i) + da*dx(i)
         dy(i + 1) = dy(i + 1) + da*dx(i + 1)
         dy(i + 2) = dy(i + 2) + da*dx(i + 2)
         dy(i + 3) = dy(i + 3) + da*dx(i + 3)
 50   continue
      return
      end
c
      double precision function ddot(n,dx,incx,dy,incy)
c-------------------------------------------------------------c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c-------------------------------------------------------------c
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
 
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c -----------------------------------------------
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
         dtemp = dtemp + dx(ix)*dy(iy)
         ix = ix + incx
         iy = iy + incy
 10   continue
      ddot = dtemp
      return
c
c code for both increments equal to 1
c
c clean-up loop
c -----------------------------------
c
 20   m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
         dtemp = dtemp + dx(i)*dy(i)
 30   continue
      if( n .lt. 5 ) go to 60
 40   mp1 = m + 1
      do 50 i = mp1,n,5
         dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *        dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + 
     *        dx(i + 4)*dy(i + 4)
 50   continue
 60   ddot = dtemp

      return
      end
c
      subroutine dgeco(a,lda,n,ipvt,rcond,z)
      integer lda,n,ipvt(1)
      double precision a(lda,1),z(1)
      double precision rcond
c
c     INTERFACE:
c     a       double precision(lda, n)
c                the matrix to be factored.
c     lda     integer
c                the leading dimension of the array  a .
c     n       integer
c                the order of the matrix  a .
c
c     RETURN VALUE:
c     a       an upper triangular matrix and the multipliers
c             which were used to obtain it.
c             the factorization can be written  a = l*u  where
c             l  is a product of permutation and unit lower
c             triangular matrices and  u  is upper triangular.
c     ipvt    integer(n)
c             an integer vector of pivot indices.
c     rcond   double precision
c             an estimate of the reciprocal condition of  a .
c             for the system  a*x = b , relative perturbations
c             in  a  and  b  of size  epsilon  may cause
c             relative perturbations in  x  of size  epsilon/rcond .
c             if  rcond  is so small that the logical expression
c                           1.0 + rcond .eq. 1.0
c             is true, then  a  may be singular to working
c             precision.  in particular,  rcond  is zero  if
c             exact singularity is detected or the estimate
c             underflows.
c     z       double precision(n)
c             a work vector whose contents are usually unimportant.
c             if  a  is close to a singular matrix, then  z  is
c             an approximate null vector in the sense that
c             norm(a*z) = rcond*norm(a)*norm(z) .
c
c     DESCRIPTION:
c     dgeco factors a double precision matrix by gaussian elimination
c     and estimates the condition of the matrix.
c
c     if  rcond  is not needed, dgefa is slightly faster.
c     to solve  a*x = b , follow dgeco by dgesl.
c     to compute  inverse(a)*c , follow dgeco by dgesl.
c     to compute  determinant(a) , follow dgeco by dgedi.
c     to compute  inverse(a) , follow dgeco by dgedi.
c
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     linpack dgefa
c     blas daxpy,ddot,dscal,dasum
c     fortran dabs,dmax1,dsign
c
c     internal variables
c
      double precision ddot,ek,t,wk,wkm
      double precision anorm,s,dasum,sm,ynorm
      integer info,j,k,kb,kp1,l
c
c
c     compute 1-norm of a
c
      anorm = 0.0d0
      do 10 j = 1, n
         anorm = dmax1(anorm,dasum(n,a(1,j),1))
 10   continue
c
c     factor
c
      call dgefa(a,lda,n,ipvt,info)
c
c     rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
c     estimate = norm(z)/norm(y) where  a*z = y  and  trans(a)*y = e .
c     trans(a)  is the transpose of a .  the components of  e  are
c     chosen to cause maximum local growth in the elements of w  where
c     trans(u)*w = e .  the vectors are frequently rescaled to avoid
c     overflow.
c
c     solve trans(u)*w = e
c
      ek = 1.0d0
      do 20 j = 1, n
         z(j) = 0.0d0
 20   continue
      do 100 k = 1, n
         if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
         if (dabs(ek-z(k)) .le. dabs(a(k,k))) go to 30
         s = dabs(a(k,k))/dabs(ek-z(k))
         call dscal(n,s,z,1)
         ek = s*ek
 30      continue
         wk = ek - z(k)
         wkm = -ek - z(k)
         s = dabs(wk)
         sm = dabs(wkm)
         if (a(k,k) .eq. 0.0d0) go to 40
         wk = wk/a(k,k)
         wkm = wkm/a(k,k)
         go to 50
 40      continue
         wk = 1.0d0
         wkm = 1.0d0
 50      continue
         kp1 = k + 1
         if (kp1 .gt. n) go to 90
         do 60 j = kp1, n
            sm = sm + dabs(z(j)+wkm*a(k,j))
            z(j) = z(j) + wk*a(k,j)
            s = s + dabs(z(j))
 60      continue
         if (s .ge. sm) go to 80
         t = wkm - wk
         wk = wkm
         do 70 j = kp1, n
            z(j) = z(j) + t*a(k,j)
 70      continue
 80      continue
 90      continue
         z(k) = wk
 100  continue
      s = 1.0d0/dasum(n,z,1)
      call dscal(n,s,z,1)
c
c     solve trans(l)*y = w
c
      do 120 kb = 1, n
         k = n + 1 - kb
         if (k .lt. n) z(k) = z(k) + ddot(n-k,a(k+1,k),1,z(k+1),1)
         if (dabs(z(k)) .le. 1.0d0) go to 110
         s = 1.0d0/dabs(z(k))
         call dscal(n,s,z,1)
 110     continue
         l = ipvt(k)
         t = z(l)
         z(l) = z(k)
         z(k) = t
 120  continue
      s = 1.0d0/dasum(n,z,1)
      call dscal(n,s,z,1)
c
      ynorm = 1.0d0
c
c     solve l*v = y
c
      do 140 k = 1, n
         l = ipvt(k)
         t = z(l)
         z(l) = z(k)
         z(k) = t
         if (k .lt. n) call daxpy(n-k,t,a(k+1,k),1,z(k+1),1)
         if (dabs(z(k)) .le. 1.0d0) go to 130
         s = 1.0d0/dabs(z(k))
         call dscal(n,s,z,1)
         ynorm = s*ynorm
 130     continue
 140  continue
      s = 1.0d0/dasum(n,z,1)
      call dscal(n,s,z,1)
      ynorm = s*ynorm
c
c     solve  u*z = v
c
      do 160 kb = 1, n
         k = n + 1 - kb
         if (dabs(z(k)) .le. dabs(a(k,k))) go to 150
         s = dabs(a(k,k))/dabs(z(k))
         call dscal(n,s,z,1)
         ynorm = s*ynorm
 150     continue
         if (a(k,k) .ne. 0.0d0) z(k) = z(k)/a(k,k)
         if (a(k,k) .eq. 0.0d0) z(k) = 1.0d0
         t = -z(k)
         call daxpy(k-1,t,a(1,k),1,z(1),1)
 160  continue
c     make znorm = 1.0
      s = 1.0d0/dasum(n,z,1)
      call dscal(n,s,z,1)
      ynorm = s*ynorm
c
      if (anorm .ne. 0.0d0) rcond = ynorm/anorm
      if (anorm .eq. 0.0d0) rcond = 0.0d0
      return
      end
c
      subroutine dgedi(a,lda,n,ipvt,det,work,job)
      integer lda,n,ipvt(1),job
      double precision a(lda,1),det(2),work(1)
c
c     INTERFACE:
c     a       double precision(lda, n)
c             the output from dgeco or dgefa.
c     lda     integer
c             the leading dimension of the array  a .
c     n       integer
c             the order of the matrix  a .
c     ipvt    integer(n)
c             the pivot vector from dgeco or dgefa.
c     work    double precision(n)
c             work vector.  contents destroyed.
c     job     integer
c                = 11   both determinant and inverse.
c                = 01   inverse only.
c                = 10   determinant only.
c
c     RETURN VALUE:
c     a       inverse of original matrix if requested.
c             otherwise unchanged.
c     det     double precision(2)
c             determinant of original matrix if requested.
c             otherwise not referenced.
c             determinant = det(1) * 10.0**det(2)
c             with  1.0 .le. dabs(det(1)) .lt. 10.0
c             or  det(1) .eq. 0.0 .
c
c     ERROR CONDITION:
c     a division by zero will occur if the input factor contains
c        a zero on the diagonal and the inverse is requested.
c        it will not occur if the subroutines are called correctly
c        and if dgeco has set rcond .gt. 0.0 or dgefa has set
c        info .eq. 0 .
c
c     DESCRIPTION:
c     dgedi computes the determinant and inverse of a matrix
c     using the factors computed by dgeco or dgefa.
C
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,dswap
c     fortran dabs,mod
c
c     internal variables
c
      double precision t
      double precision ten
      integer i,j,k,kb,kp1,l,nm1
c
c
c     compute determinant
c
      if (job/10 .eq. 0) go to 70
      det(1) = 1.0d0
      det(2) = 0.0d0
      ten = 10.0d0
      do 50 i = 1, n
         if (ipvt(i) .ne. i) det(1) = -det(1)
         det(1) = a(i,i)*det(1)
c     ...exit
         if (det(1) .eq. 0.0d0) go to 60
 10      if (dabs(det(1)) .ge. 1.0d0) go to 20
         det(1) = ten*det(1)
         det(2) = det(2) - 1.0d0
         go to 10
 20      continue
 30      if (dabs(det(1)) .lt. ten) go to 40
         det(1) = det(1)/ten
         det(2) = det(2) + 1.0d0
         go to 30
 40      continue
 50   continue
 60   continue
 70   continue
c     
c     compute inverse(u)
c
      if (mod(job,10) .eq. 0) go to 150
      do 100 k = 1, n
         a(k,k) = 1.0d0/a(k,k)
         t = -a(k,k)
         call dscal(k-1,t,a(1,k),1)
         kp1 = k + 1
         if (n .lt. kp1) go to 90
         do 80 j = kp1, n
            t = a(k,j)
            a(k,j) = 0.0d0
            call daxpy(k,t,a(1,k),1,a(1,j),1)
 80      continue
 90      continue
 100  continue
c
c        form inverse(u)*inverse(l)
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 140
      do 130 kb = 1, nm1
         k = n - kb
         kp1 = k + 1
         do 110 i = kp1, n
            work(i) = a(i,k)
            a(i,k) = 0.0d0
 110     continue
         do 120 j = kp1, n
            t = work(j)
            call daxpy(n,t,a(1,j),1,a(1,k),1)
 120     continue
         l = ipvt(k)
         if (l .ne. k) call dswap(n,a(1,k),1,a(1,l),1)
 130  continue
 140  continue
 150  continue
      return
      end
c
      subroutine dgefa(a,lda,n,ipvt,info)
c
c     INTERFACE:
c     a       double precision(lda, n)
c             the matrix to be factored.
c     lda     integer
c             the leading dimension of the array  a .
c     n       integer
c             the order of the matrix  a .
c
c     RETURN VALUE:
c     a       an upper triangular matrix and the multipliers
c             which were used to obtain it.
c             the factorization can be written  a = l*u  where
c             l  is a product of permutation and unit lower
c             triangular matrices and  u  is upper triangular.
c     ipvt    integer(n)
c             an integer vector of pivot indices.
c     info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
c
c     DESCRIPTION:
c     dgefa factors a double precision matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     REFERENCES:
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c
c     internal variables
c
      integer lda,n,ipvt(1),info
      double precision a(lda,1)
      double precision t
      integer idamax,j,k,kp1,l,nm1
c
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0d0) go to 40
c
c           interchange if necessary
c
         if (l .eq. k) go to 10
         t = a(l,k)
         a(l,k) = a(k,k)
         a(k,k) = t
 10      continue
c
c           compute multipliers
c
         t = -1.0d0/a(k,k)
         call dscal(n-k,t,a(k+1,k),1)
c     
c           row elimination with column indexing
c
         do 30 j = kp1, n
            t = a(l,j)
            if (l .eq. k) go to 20
            a(l,j) = a(k,j)
            a(k,j) = t
 20         continue
            call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
 30      continue
         go to 50
 40      continue
         info = k
 50      continue
 60   continue
 70   continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0d0) info = n
      return
      end
c
      subroutine dgesl(a,lda,n,ipvt,b,job)
c
c     INTERFACE:
c     a       double precision(lda, n)
c             the output from dgeco or dgefa.
c     lda     integer
c             the leading dimension of the array  a .
c     n       integer
c             the order of the matrix  a .
c     ipvt    integer(n)
c             the pivot vector from dgeco or dgefa.
c     b       double precision(n)
c             the right hand side vector.
c     job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     RETURN VALUE:
c     b       the solution vector  x .
c
c     ERROR CONDITION:
c     a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgeco has set rcond .gt. 0.0
c        or dgefa has set info .eq. 0 .
c
c
c     DESCRIPTION:
c     dgesl solves the double precision system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgeco or dgefa.
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     REFERENCES:
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c
c     internal variables
c
      integer lda,n,ipvt(1),job
      double precision a(lda,1),b(1)
      double precision ddot,t
      integer k,kb,l,nm1
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
      if (nm1 .lt. 1) go to 30
      do 20 k = 1, nm1
         l = ipvt(k)
         t = b(l)
         if (l .eq. k) go to 10
         b(l) = b(k)
         b(k) = t
 10      continue
         call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
 20   continue
 30   continue
c
c        now solve  u*x = y
c
      do 40 kb = 1, n
         k = n + 1 - kb
         b(k) = b(k)/a(k,k)
         t = -b(k)
         call daxpy(k-1,t,a(1,k),1,b(1),1)
 40   continue
      go to 100
 50   continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
      do 60 k = 1, n
         t = ddot(k-1,a(1,k),1,b(1),1)
         b(k) = (b(k) - t)/a(k,k)
 60   continue
c
c        now solve trans(l)*x = y
c
      if (nm1 .lt. 1) go to 90
      do 80 kb = 1, nm1
         k = n - kb
         b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
         l = ipvt(k)
         if (l .eq. k) go to 70
         t = b(l)
         b(l) = b(k)
         b(k) = t
 70      continue
 80   continue
 90   continue
 100  continue
      return
      end
c
      subroutine  dscal(n,da,dx,incx)
c----------------------------------------------------------------c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c----------------------------------------------------------------c
      double precision da,dx(1)
      integer i,incx,m,mp1,n,nincx
 
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c ---------------------------------
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
         dx(i) = da*dx(i)
 10   continue
      return
c
c code for increment equal to 1
c
c clean-up loop
c -----------------------------
c
 20   m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
         dx(i) = da*dx(i)
 30   continue
      if( n .lt. 5 ) return
 40   mp1 = m + 1
      do 50 i = mp1,n,5
         dx(i) = da*dx(i)
         dx(i + 1) = da*dx(i + 1)
         dx(i + 2) = da*dx(i + 2)
         dx(i + 3) = da*dx(i + 3)
         dx(i + 4) = da*dx(i + 4)
 50   continue

      return
      end
      subroutine  dswap (n,dx,incx,dy,incy)
c
c     interchanges two vectors.
c     uses unrolled loops for increments equal one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c       code for unequal increments or equal increments not equal
c         to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
         dtemp = dx(ix)
         dx(ix) = dy(iy)
         dy(iy) = dtemp
         ix = ix + incx
         iy = iy + incy
 10   continue
      return
c
c       code for both increments equal to 1
c
c
c       clean-up loop
c
 20   m = mod(n,3)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
         dtemp = dx(i)
         dx(i) = dy(i)
         dy(i) = dtemp
 30   continue
      if( n .lt. 3 ) return
 40   mp1 = m + 1
      do 50 i = mp1,n,3
         dtemp = dx(i)
         dx(i) = dy(i)
         dy(i) = dtemp
         dtemp = dx(i + 1)
         dx(i + 1) = dy(i + 1)
         dy(i + 1) = dtemp
         dtemp = dx(i + 2)
         dx(i + 2) = dy(i + 2)
         dy(i + 2) = dtemp
 50   continue
      return
      end

      integer function idamax(n,dx,incx)
c-----------------------------------------------------------------c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c-----------------------------------------------------------------c
      double precision dx(1),dmax
      integer i,incx,ix,n
c
      idamax = 0
      if( n .lt. 1 ) return
      idamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c ---------------------------------
c
      ix = 1
      dmax = dabs(dx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(dabs(dx(ix)).le.dmax) go to 5
         idamax = i
         dmax = dabs(dx(ix))
    5    ix = ix + incx
 10   continue
      return
c
c code for increment equal to 1
c -----------------------------
c
 20   dmax = dabs(dx(1))
      do 30 i = 2,n
         if(dabs(dx(i)).le.dmax) go to 30
         idamax = i
         dmax = dabs(dx(i))
 30   continue
      return
      end