c
c programme inversion par gradient
c avec pregradient integre et ponderation des fichiers de donnnees
c
c janv2000 - essaie de fixer des patchs a des valeurs fixees par modele a priori
c mars 2000 - augmentation de delta tous les 10 iterations
c 31 mars 2000 modif pour fichier de sortie lorsque nrec!=1  (*alea(i,j))
c 
       implicit real*8 (a-h,o-z)
       parameter (mxd=50000,mxs=200,maxrec=10,mxf=9)
cc mxs=maxi de sources, mxf=maxi de fichiers de donnees, mxd=maxi de donnees
       character*100 fires,form,formf,formo,fida,fimo
       character*100 nom1,nom2,fisou,fault,nom3,nom4,nom5
c
       dimension a(mxd,mxs),grad(mxs)
       dimension bs(mxd),x(mxs)
       dimension alea(maxrec,mxs),io(maxrec)
       dimension bests(maxrec),bestx(maxrec,mxs)
       dimension at(mxs,mxd),tapar(mxs,10),libre(mxs)
       dimension xdata(mxd),ydata(mxd),ddata(mxd)
       dimension weight(mxf),edata(mxd)
       dimension nv(mxf,2),irv(mxf)
       dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2)
       dimension ixxx(mxs,2),iyyy(mxs,2)
       dimension xf(mxs),yf(mxs),ff(mxs),h(mxs)
       dimension w(mxs),r(mxs),tf(mxs)
       dimension us(mxs),ud(mxs),ut(mxs)
       dimension t_momu(mxs)
       real*4, rand_unif
c
c
        common/param/tapar,numsour,libre
        common/data1/xdata,ydata,ddata,edata,weight
        common/data2/nv,tsa,tsb,tsc,irv,dref
        common/data3/nfic
        common/defo/a
        COMMON/F1/XF,YF,FF,H,W,R,TF
        COMMON/F2/US,UD,UT
        COMMON/COMOB/Q,ST,CT,ELAS,U0
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
        common/grid_mom/ixxx,iyyy
c
       pi=4*datan(1.d0)
       smu  = 3.0d10
       elas=0.5d0
       cc=6399593.6258000d0
       AK=0.9996d0
       eprim=0.08209443779496d0
       EP2=EPRIM*EPRIM
       E2=EP2
c
       idimax=100
        idim=0
c
        open(99,file='GRAD.LOG')
        call pregradient (ndat,nstri,ndip,igeo)
        print*,'FIN CALCUL DIRECT'
c
c
        numx=nstri
        numy=ndip
        npar=numsour
c
        do 50 jj=1,nfic 
        do 51 kk=nv(jj,1),nv(jj,2)
              xx=xdata(kk)/1000.
              yy=ydata(kk)/1000.
           if (igeo.eq.2) then
                call utmgeo (xx,yy,xx,yy)
           endif
              do ll=1,npar
                 x(ll)=1.d0
              enddo
 51     enddo
 50     enddo
c
c
          call mvmult(mxd,mxs,ndat,npar,a,x,bs)
        do k=1,npar
           t_momu(k)=tapar(k,8)
        enddo
        print*,ndat,' donnees ; ',npar,' parametres'
        write(99,*)ndat,' donnees ; ',npar,' parametres'
c
        write(99,*)' SOURCES' 
        print*,'SOURCES'
        do kd=1,numy
           k1=(kd-1)*numx+1
           k2=k1+numx-1   
           write(99,*)(real(tapar(j,8)),j=k1,k2)
           write(*,*)(real(tapar(j,8)),j=k1,k2)
        enddo
c
       print*,'nb de recherches aleatoires <',maxrec
       read(*,*)nrec
       write(99,*)'nb de recherches aleatoires <',maxrec
       nrec=nrec+1
c
       print*,' glissement min, max, increment (cm)'
       read(*,*)smin,smax,sinc
       write(99,*)' glissement min, max, increment (cm)'
       write(99,*)sngl(smin),sngl(smax),sngl(sinc)
c
       print*,'nb d"iterations pour gradient'
       read(*,*)niter
       write(99,*)' nb d"iterations pour gradient',niter
c
      k=1
      do j=1,npar
         alea(k,j)=1.d0
      enddo
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      irec=0
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ BOUCLE SUR les REC
 300  continue
      irec=irec+1
      write(99,*)''
      write(99,*)'' 
      write(99,*)'+++++++++++++++++++'
      write(99,*)' IREC ',irec
      idim=0
      delta=1

c
c      print*,'ALEA'
      do 13 j=1,npar
c         print*,alea(irec,j)
         x(j)=alea(irec,j)
         bestx(irec,j)=x(j)
 13   enddo
c
c calcul de l'erreur quadratique
c
      call mvmult(mxd,mxs,ndat,npar,a,x,bs)
      ssr=0.d0
      sw=0.d0
c       call calcssr(ndat,ddata,bs,ssr)
      write(99,*)''
      write(99,*)' CALCUL DE L ERREUR INITIALE'
      write(99,*)'no fich,  rms fich,          rms tout'
      do if=1,nfic
        n1=nv(if,1)
        n2=nv(if,2)
        call calcssrt(ddata,bs,n1,n2,ssrs)
        ssr=ssr+ssrs*weight(if)
        sw=sw+weight(if)*(n2-n1+1)
        ww=weight(if)
        print*,if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/(sw))
      write(99,*)if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/(sw))
      enddo
      print*,'RMS TOT',dsqrt(ssr/sw)
      write(99,*)'RMS TOT ',dsqrt(ssr/sw)
c
      bests(irec)=ssr
c
      write(99,*) ' '
      print*,'DELTA INI ',delta
      write(99,*)'DELTA INI ',real(delta)
      write(*,'(1000f8.1)')(x(j)*tapar(j,8),j=1,npar)
      write(99,'(1000f8.1)')(real(x(j)*tapar(j,8)),j=1,npar)
      print*,'PAR INI'
      do j=1,npar
           print*,x(j)
      enddo
      
c
c
      call mtxtrp(mxd,mxs,ndat,npar,a,at)
c
 12   continue
c
      increase_delta=0
c
c****************** 
c*** ITER *******
c******************
      do 100 iter=1,niter
        increase_delta=increase_delta+1
        if (increase_delta.gt.10) then
            print*,'INCREASE DELTA ', delta*2.
            delta=delta*2.
            increase_delta=0
        endif            
c
        print*,''
        print*,'ITER',iter
        write(99,*)''
        write(99,*)'ITER',iter
        call findgrad(mxd,mxs,ndat,npar,a,at,ddata,x,grad)
c
        do 200 j=1,npar
          if (libre(j).ne.0) then
              x(j)=x(j)+delta*grad(j)
          endif
          slipr=x(j)*tapar(j,8)
          if(slipr.lt.0.0) x(j)=0.d0
 200    continue
        call mvmult(mxd,mxs,ndat,npar,a,x,bs)
      ssr=0.d0
      sw=0.d0
      write(99,*)'no fich,   rms fich,      rms tout'
      do if=1,nfic
        n1=nv(if,1)
        n2=nv(if,2)
        call calcssrt(ddata,bs,n1,n2,ssrs)
        ssr=ssr+ssrs*weight(if)
        sw=sw+weight(if)*(n2-n1+1)
        ww=weight(if)  
c       print*,if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/sw)
      write(99,*)if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/sw)
      enddo
      print*,'RMS TOT',dsqrt(ssr/sw)
      write(99,*)'RMS TOT ',dsqrt(ssr/sw)
c
c       write(99,*)' SOURCES'
c       do kd=1,numy 
c          k1=(kd-1)*numx+1
c          k2=k1+numx-1
c          write(99,*)(real(x(j)*tapar(j,8)),j=k1,k2) 
c          write(99,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2)
c          write(*,*)(real(x(j)*tapar(j,8)),j=k1,k2)
c          write(*,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2)
c       enddo
c
        if (ssr.lt.bests(irec)) then
          brms=dsqrt(ssr/sw)
          bests(irec)=ssr
          do j=1,npar
             bestx(irec,j)=x(j)
          enddo
        else
          idim=idim+1
          print*,'idim DIMIN DELTA',idim,delta*.5
          write(99,*)'idim,DIMIN DELTA', idim,delta*.5
          increase_delta=0
          ssr=bests(irec)
          delta=delta*.5
          do j=1,npar
            x(j)=bestx(irec,j)
          enddo
          iter=iter-1
        endif
c
c
        if(idim.gt.idimax) then
         write(99,*)' MAXIMUM DIMINUTION'
         write(99,*)' idimax idim',idimax,idim
         print*,' MAXIMUM DIMINUTION'
         print*,'idimax ',idim
         goto 101
        endif
 100    continue
 101    continue
        print*,iter, ' iterations'
        write (99,*)''
        write (99,*)' NB ITERATIONS EFFECTUEES', iter
c
cccc
c998   continue
c       brms=bests(irec)/sw
c       print*,' essais autres valeurs de slip'
c       do j=1,npar
c          print*,real(bestx(irec,j)*tapar(j,8))
c          print*,' slip ?'
c          read(*,*)slip
c          x(j)=slip/tapar(j,8)
c       enddo
c       call mvmult(mxd,mxs,ndat,npar,a,x,bs)
c       do if=1,nfic
c         n1=nv(if,1)
c         n2=nv(if,2)
c         call calcssrt(ddata,bs,n1,n2,ssrs)
c         ssr=ssr+ssrs*weight(if)
c         sw=sw+weight(if)*(n2-n1+1)
c         ww=weight(if)
c       enddo
c       print*,'slip RMS ',slip, dsqrt(ssr/sw)
c       print*,'BEST', (bestx(irec,1)*tapar(1,8)),brms
c       if(ssr.lt.bests(irec)) then
c         bests(irec)=ssr
c         brms=dsqrt(ssr/sw)
c         do j=1,npar
c            bestx(irec,j)=x(j)
c         enddo
c       endif
c       goto 998
c          
c
c----------------------meilleur modele
c
        print*,'MEILLEUR MODELE'
        print*,bests(irec),dsqrt(bests(irec)/sw)
        write(*,'(1000f8.1)')(bestx(irec,j)*tapar(j,8),j=1,npar)
c
        write(99,*)''
        write(99,*)'meilleur modele'
        write(99,*)'RMS ',real(dsqrt(bests(irec)/sw))
        write(99,*)' SOURCES' 
        print*,' SOURCES'
        do kd=1,numy
           k1=(kd-1)*numx+1
           k2=k1+numx-1
           write(99,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2)
           write(*,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2) 
        enddo
c
c
c-----------------------recherche modele aleatoire
      print*,''  
        if(irec.ge.nrec) then
           goto 400
        endif
        do k=1,npar
         tmax=smax/tapar(k,8)
         tmin=smin/tapar(k,8)
         tdif=tmax-tmin
         tinc=sinc/tapar(k,8)
         ns=idnint(tdif/tinc)+1
         alea(irec+1,k)=tmin+dble(nint(ns*rand_unif(0)))*tinc
        enddo
        goto 300
c
 400    continue
        print*,'nrec ',nrec
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ BOUCLE SUR les REC 
c

c-------------------------------------------------
c-------------tri et ecriture des modeles aleatoires recherches
c
        do k=1,nrec
           io(k)=k
        enddo
c
        call tri (bests,io,nrec)
c
c
        open(51,file='resu')
        nw=npar*4
        open(52,file='resu.bin',access='direct',recl=nw)
c
c
       write(51,*)'nb modeles sources ',nrec,npar
       do k=1,nrec
         write(51,*)bests(k),dsqrt(bests(irec)/ndat)
       write(51,'(1000f8.1)')(bestx(k,j)*tapar(j,8),j=1,npar)
c      write(*,'(1000f8.1)')(bestx(k,j)*tapar(j,8),j=1,npar)
       write(52,rec=k)(sngl(bestx(k,j)*tapar(j,8)),j=1,npar)
        enddo
        close(51)
        close(52)
c
        nom1='resi'
        nom2='sour'
        nom3='out'
        nom4='data'
        nom5='modl'
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
        write(99,*)''
        write(99,*)'*********** SORTIE ****************' 
c
        do 60 k=1,nrec
          do j=1,npar
             tapar(j,8)=t_momu(j)*bestx(k,j)*alea(k,j)
          enddo
c
          call direct 
           do j=1,npar
              x(j)=1.0
           enddo
           call mvmult(mxd,mxs,ndat,npar,a,x,bs)
c
         if(k.lt.10) then
          form='(a4,i1,a4)'
          formf='(a3,i1,a6)'
          formo='(a4,i1,i1,a4)'
         elseif(k.lt.100) then
          form='(a4,i2,a4)' 
          formf='(a3,i2,a6)' 
          formo='(a4,i2,i1,a4)'  
         else
          form='(a4,i3,a4)'
          formf='(a3,i3,a6)' 
          formo='(a4,i2,i1,a4)'  
          endif
c
c
         rmst=0.d0
         rmsp=0.d0
         sw=0.d0
         do 61 jj=1,nfic
           write(99,*)''
           print*,'FICHIER',jj
           write(99,*)'FICHIER',jj
           write(fida,formo)nom4,k,jj,'.out'
           write(fimo,formo)nom5,k,jj,'.out'
           write(fires,formo)nom1,k,jj,'.out'
           open(56,file=fida)
           open(57,file=fimo)
           open(53,file=fires)
c
           rms=0.0
           nndd=0
           write(99,*)' reference donnees modl',dref(jj,1),dref(jj,2)
c           
           do 62 kk=nv(jj,1),nv(jj,2)
              nndd=nndd+1
              if(igeo.eq.2) then
                call utmgeo (xdata(kk)/1000.,ydata(kk)/1000.,xx,yy)
              else
                xx=xdata(kk)/1000.
                yy=ydata(kk)/1000.
              endif 
            deff=bs(kk)+dref(jj,2)
            ddaa=ddata(kk)+dref(jj,2)
            adiff=ddaa-deff
            write(56,*)xx,yy,float(ddaa)
            write(57,*)xx,yy,float(deff)
            write(53,*)xx,yy,float(adiff)
            rms=adiff**2+rms
 62        enddo
c
c
           rmsp=rms*weight(jj)+rmsp
           sw=sw+weight(jj)*(nv(jj,2)-nv(jj,1)+1)
           rmst=rms+rmst
c
c
           rms=dsqrt(rms/dble(nndd))
           write(53,*)' RMS= ',rms
           write(99,*)'ndata, RMS, poids ',nndd,rms,weight(jj)
           print*,'ndata RMS= ',nndd,rms
c
c
           close(53)
           close(56)
           close(57)
 61     enddo
c
c
        rmst=dsqrt(rmst/dble(ndat))
        rmsp=dsqrt(rmsp/sw)
        print*,'rms tot',rmst
        write(99,*)''
        write(99,*)' ndat RMS TOTAL ',ndat,rmst
        write(99,*)' RMS PONDERE    ',rmsp
c
c------- fichier sour_k_.out (source en binaire)
 
         write(fisou,form)nom2,k,'.out'
         print*,fisou
         nw=4*numy
         open(54,file=fisou,access='direct',recl=nw)
         l=0
         do 73 j=numx,1,-1
           l=l+1
c       write(*,'(1000f7.1)')(*tapar(i,8),
c    &  i=j,j+numx*(numy-1),numx)
        write(*,*)(real(tapar(i,8))
     &  ,i=j,j+numx*(numy-1),numx)
        write(54,rec=l)(real(tapar(i,8))
     &  ,i=j,j+numx*(numy-1),numx)
 73     enddo
         close(54)
 60   enddo
c
c----------------------------Ecriture sources----
c
         do 72 i=1,npar
         call jcr2new(tapar(i,1),tapar(i,2),tapar(i,5),
     &    tapar(i,6),tapar(i,3),tapar(i,4))
           if(igeo.eq.2) then
             call utmgeo(tapar(i,1),tapar(i,2),tapar(i,1),tapar(i,2))
           endif    
 72     enddo
c
        do 71 k=1,nrec
c          print*,"RECHERCHE ",K
c
         do 74 i=1,npar
             tapar(i,8)=t_momu(i)*bestx(k,i)*alea(k,i)
c              print*,t_momu(i),bestx(k,i),tapar(i,8)
c           print*,tapar(i,1),tapar(i,2)
 74     enddo
c
         write(fault,formf)nom3,k,'.fault'
         print*,fault
         open(55,file=fault)
         qq=0.d0
         write(55,*)npar,' 5 0 0 '
         do i=1,npar
c
         write(55,*)i
         write(55,*)tapar(i,1),qq,tapar(i,2),qq
         write(55,*)tapar(i,3),qq,tapar(i,4),qq
         write(55,*)tapar(i,5),qq,tapar(i,6),qq,tapar(i,7),qq
         write(55,*)tapar(i,8)+0.001,qq,tapar(i,9),qq,tapar(i,10),qq
         enddo
         close(55)
c
c------- fichier sour_k_.out (source en binaire)
c
c        write(fisou,form)nom2,k,'.out'
c        print*,fisou
c        nw=4*numy
c        open(54,file=fisou,access='direct',recl=nw)
c        l=0
c        do 73 j=numx,1,-1
c          l=l+1
c       write(*,'(1000f7.1)')(tapar(i,8),
c    &  i=j,j+numx*(numy-1),numx)
c       write(*,*)(tapar(i,8)
c    &  ,i=j,j+numx*(numy-1),numx)
c       write(54,rec=l)(real(tapar(i,8))
c    &  ,i=j,j+numx*(numy-1),numx)
c73     enddo
c        close(54)
c

 71   enddo
c
c-----------------------
       stop
       end
c
c---------------------------------
      subroutine  mtxtrp(mxd,mxs,ndat,npar,a,at) 
c
       implicit real*8 (a-h,o-z)
       dimension at(mxs,mxd)
       dimension a(mxd,mxs)
c
       do  i=1,ndat
         do   j=1,npar
           at(j,i)=a(i,j)
         enddo
       enddo
       return
       end
c------------------------------------
c
        subroutine findgrad(mxd,mxs,ndat,
     &  npar,a,at,ddata,x,grad)
c
       implicit real*8 (a-h,o-z)
        dimension a(mxd,mxs),at(mxs,mxd)
        dimension x(mxs),ddata(mxd)
        dimension grad(mxs)
        dimension passa(100000)
c
        do 10 i=1,ndat
          sum=0.d0
           do 20 j=1,npar
             sum=sum+a(i,j)*x(j)
  20       continue
           passa(i)=ddata(i)-sum
c          print*,'passa ',passa(i)
  10    continue
        call mvmult(mxs,mxd,npar,ndat,at,passa,grad)
        return
        end
c------------------------------------
c
       subroutine mvmult(max1,max2,n1,n2,a,v,r)
c
       implicit real*8 (a-h,o-z)
       parameter (mxd=50000,mxs=200)
c
       dimension a(max1,max2),v(1),r(1)
c
       do 10 i=1,n1
         sum=0.d0
         do 20 j=1,n2
            sum=sum+a(i,j)*v(j)
 20      continue
         r(i)=sum
 10      continue
       return
       end
c
c-----------------------------------------
c
       subroutine calcssr(ndat,b,bs,ssr)
c
       implicit real*8 (a-h,o-z)
       dimension b(1),bs(1)
c
       ssr=0.d0
       do 10 i=1,ndat
          dif=bs(i)-b(i)
          ssr=ssr+dif**2
 10    continue
       return
       end
c

c-----------------------------------------
c
       subroutine calcssrt(b,bs,n1,n2,ssr)
c
       implicit real*8 (a-h,o-z) 
c
       dimension bs(1),b(1)
c
c
       ssr=0.d0
       do 10 i=n1,n2
          dif=bs(i)-b(i)
c         PRINT*,bs(i),b(i)
          ssr=ssr+dif**2
 10    continue
       return
       end
c
c------------------------------------------------
c
c       subroutine ecr_par(igeo,numx,numy)
c
c      implicit real*8 (a-h,o-z)
c
c      parameter (mxs=200)
c       parameter (mxf=9)
c
c       character*100 fpar
c      dimension tapar(mxs,10)
c
c       common/param/tapar,numsour
c
c       print*,'nom du fichier sortie parametre de faille'
c       read(*,'(a)')fpar
c       open (40,file=fpar)
c
c       nnuummss=0
c       write(40,*)numsour,' 5  0  0 '
c       do k=1,numsour
c          nnuummss=nnuummss+1
c          write(40,*)nnuummss
c         write(40,*)tapar(k,1),zz,tapar(k,2),zz
c         write(40,*)tapar(k,3),zz,tapar(k,4),zz
c         write(40,*)tapar(k,5),zz,tapar(k,6),zz,tapar(k,7),zz
c         write(40,*)tapar(k,8),zz,tapar(k,9),zz,tapar(k,10),zz
c       enddo
c
c
c       close(40)
c       return
c       end
c
c--------------------------------------------------------
c
        subroutine  maxv (mxs,npar,grad,vmin,vmax) 
c
c
       implicit real*8 (a-h,o-z)
       dimension grad(mxs)
c
       vmin=1e50
       vmax=-1e50
       do j=1,npar
         vmin=dmin1(vmin,grad(j))
         vmax=dmax1(vmax,grad(j))  
       enddo
       return
       end

c
c
C****************************************************************************
c
       SUBROUTINE RAND_INIT
c
       integer time, isec
c
       real*4 dtime,tarray(2),adum
c
       INTEGER*4 I,J
c
       COMMON /RANDCOM/ INITFLAG,RANTAB,INDEX1,INDEX2,INDEX3
       LOGICAL INITFLAG
       INTEGER*4 RANTAB(0:54)
       INTEGER*4 INDEX1,INDEX2,INDEX3
       DATA INITFLAG /.FALSE./
c
c   get the year,month,day,hour,min,sec from time() and calculate a seed
c   based on this value [secs] (relative to 1970) plus the time elapsed for
c   executing a 'waiting' loop - to make sure that two subsequent calls to
c   RAND_INIT don't yield the same seed
c
      isec = int(dtime(tarray))
      isec = time()
      adum = 0.01
      do j = 1,1000000
        adum = atan(atan(atan(atan(atan(adum)))))
      enddo
      isec = isec + int(dtime(tarray)*1.e6)
      RANTAB(0) = isec
c
C
C   Use a linear congruential psuedo random number generator to
C   initialize the table for the RAND_$UNIF function. Note that
C   Knuth's rules for 'good' generators requires a multiplier of
C   the form ending in ...x21 where 'x' is an even digit. The
C   choice of 31415821 is taken from one of his examples.
C
       DO 100 I = 1,54
         RANTAB(I) = RANTAB(I-1)*31415821+1
 100   CONTINUE
c
       INITFLAG = .TRUE.
       INDEX1 = 0
       INDEX2 = 23
       INDEX3 = 54
 
       RETURN
       END

C****************************************************************************
****************************************************************************
C
C           Uniform Random Number Generator.
C           Generates a psuedo-random number in the range 0.0 to 1.0 with
C           a uniform distribution. The first time RAND_$UNIF is called the
C           argument is used to seed the random number generator if the
C           random number generator has not already been seeded by a call
C           to the RAND_$INIT routine. If the seed passed on the first call
C           is greater than 0, then the seed is used to initialize the
C           random number generator, otherwise a built in seed number is
C           used to set up the initial table. Since IEEE standard floating
C           point numbers have a 24 bit mantissa (one bit being the hidden bit)
C           we generate random integers in the range 0 to 2**24-1 and then
C           divide them by 2**24-1 to get a real number in the range 0 to
C           1.0, inclusive.
C
C****************************************************************************
 
      REAL*4 FUNCTION RAND_UNIF (ISEED)
 
      COMMON /RANDCOM/ INITFLAG,RANTAB,INDEX1,INDEX2,INDEX3
      LOGICAL INITFLAG
      INTEGER*4 RANTAB(0:54)
      INTEGER*4 INDEX1,INDEX2,INDEX3
 
      INTEGER*4 ISEED
 
C
C   If random number generator has not already been seeded
C   by calling the RAND_$INIT routine, then set up the initial
C   table for the additive congruential method used here. If
C   the seed number supplied by the caller is 0, then use our
C   own seed number.
C
      IF (INITFLAG) GOTO 100
 
      IF (ISEED.GT.0) THEN
        RANTAB(0) = ISEED
      ELSE
        RANTAB(0) = 31415926
      ENDIF
C
C   Use a linear congruential psuedo random number generator to
C   initialize the table for the RAND_$UNIF function. Note that
C   Knuth's rules for 'good' generators requires a multiplier of
C   the form ending in ...x21 where 'x' is an even digit. The
C   choice of 31415821 is taken from one of his examples.
C
 
      DO 10 I = 1,54
        RANTAB(I) = RANTAB(I-1)*31415821+1
10    CONTINUE
 
      INITFLAG = .TRUE.
      INDEX1 = 0
      INDEX2 = 23
      INDEX3 = 54
 
C
C   Generate the next random integer in the range 0 to 2**24-1
C   and then convert it into a real number in the range 0.0 to 1.0
C
100   INDEX1 = MOD(INDEX1+1,55)
      INDEX2 = MOD(INDEX2+1,55)
      INDEX3 = MOD(INDEX3+1,55)
      RANTAB(INDEX1) = RANTAB(INDEX2)+RANTAB(INDEX3)
 
      RAND_UNIF = FLOAT(IAND(RANTAB(INDEX1),2**24-1))/FLOAT(2**24-1)
      RETURN
      END
c
c------------------------------------
       subroutine tri(tab,io,n)
c
       implicit real*8(a-h,o-z)
       logical drap
       dimension tab(1),io(1)
c
 200   drap=.true.
c
       do 201 i=1,n
        if(i+1.gt.n)goto 201
        if(tab(i).gt.tab(i+1)) then
c
          bidu=tab(i)
          tab(i)=tab(i+1)
          tab(i+1)=bidu
c
          ibi=io(i)
          io(i)=io(i+1)
          io(i+1)=ibi
c
          drap=.false.
        endif
 201    continue
        if(.not.drap) then
         goto 200
        endif
c
        return
        end
      
c-----------------------------------------------------------
       subroutine pregradient(ndat,numx,numy,igeo)
c
       implicit real*8 (a-h,o-z)
c
       parameter (mxs=200,mxd=50000)
       parameter (mxf=9)
c
c mxs=maxi de sources, mxf=maxi de fichiers de donnees, mxd=maxi de donnees
c
c
       dimension xdata(mxd),ydata(mxd),ddata(mxd)
       dimension weight(mxf),edata(mxd)
       dimension nv(mxf,2),irv(mxf)
       dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2)
       dimension a(mxd,mxs)
       dimension tapar(mxs,10),libre(mxs)
       dimension ixxx(mxs,2),iyyy(mxs,2)
       dimension xf(mxs),yf(mxs),ff(mxs),h(mxs)
       dimension w(mxs),r(mxs),tf(mxs)
       dimension us(mxs),ud(mxs),ut(mxs)
c
        common/param/tapar,numsour,libre
        common/data1/xdata,ydata,ddata,edata,weight
        common/data2/nv,tsa,tsb,tsc,irv,dref
         common/data3/nfic 
        common/defo/a
        COMMON/F1/XF,YF,FF,H,W,R,TF
        COMMON/F2/US,UD,UT
        COMMON/COMOB/Q,ST,CT,ELAS,U0
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
        common/grid_mom/ixxx,iyyy
c
c
c
c ouverture fichier LOG
c
       call lec_par (igeo,numx,numy)
c
c
       print*,'nombre de fichiers de donnees'
       read(*,*)nfic
       write(99,*)nfic,' fichiers de donnees'
c
       call lec_data (igeo,numx,numy)
c      
       ndat =nv(nfic,2)
c
          print*,'DIRECT'
       call direct
c
c
       sigma2=0.d0
       do i=1,ndat 
           sigma2=sigma2+(ddata(i))**2
       enddo
       sigma2=dsqrt(sigma2)
       sigma2=1.0
       print*,' COV= ',sigma2
c
c
cc...print
       write(99,*)''
      print*,' '
        print*,'DATA no fich , n1, n2, nn, reference rel, abs, poids'
        write(99,*)'DATA no , n1, n2, nn, reference, poids'
      do ii=1,nfic
        n1=nv(ii,1)
        n2=nv(ii,2)
        nn=n2-n1+1
        if(irv(ii).ne.0) then
          nref=irv(ii)
          mref=nv(ii,1)+irv(ii)-1
        else
          nref=0
          mref=0
        endif
        write(99,*)ii,n1,n2,nn,nref,mref,sngl(weight(ii))
        print*,ii,n1,n2,nn,nref,mref,sngl(weight(ii))
      enddo
c
        return
        end
c
c----------------------------------------
        subroutine lec_par (igeo,numx,numy)
c-------------------------------
c
        implicit real*8 (a-h,o-z)
c
        parameter (mxs=200,mxd=50000)
        parameter (mxf=9)
c
        character*100 fpar
       dimension tapar(mxs,10),libre(mxs)
c
        common/param/tapar,numsour,libre
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
        print*,'nom du fichier parametre de faille'
        read(*,'(a)')fpar
        write(99,*)fpar
c
        print*,' 2 : geographique' 
        print*,' 1 : UTM' 
        read(*,*)igeo
        if(igeo.eq.2) then
          print*,' fuseau'
          read(*,*)ifu
          write(99,*)'geographique , fuseau ',ifu
        endif
c
        print*,'nombre de failles along strike'
        read(*,*)numx
        print*,'nombre de failles along dip'
        read(*,*)numy
        write(99,*)'nb failles en X et Y ',numx,numy
c
        print*,' voulez vous normaliser le slip 1/0'
       read(*,*)lnorm
        if(lnorm.eq.1) then
             print*,'valeur du slip (en cm)'
             read(*,*)slipnorm
             write(99,*)'slip normalise ', slipnorm
         else
             write(99,*)'slip non normalise'
        endif
c
c
        open (30,file=fpar)
c
        read(30,*)numsour
        nnppaarr=0
        do k=1,numsour
           read(30,*)nnuummss
           write(99,*)''
           write(99,*)'source no ',nnuummss
           read(30,*)tapar(k,1),zz,tapar(k,2),zz
           read(30,*)tapar(k,3),zz,tapar(k,4),zz
           read(30,*)tapar(k,5),zz,tapar(k,6),zz,tapar(k,7),zz
           read(30,*)tapar(k,8),wli,tapar(k,9),zz,tapar(k,10),zz 
           libre(k)=nint(wli)
           if(lnorm.eq.1.and.libre(k).ne.0) tapar(k,8)=slipnorm
           write(99,*)tapar(k,1),tapar(k,2)
           write(99,*)tapar(k,3),tapar(k,4)
           write(99,*)tapar(k,5),tapar(k,6),tapar(k,7)
           write(99,*)tapar(k,8),tapar(k,9),tapar(k,10) 
           if(libre(k).ne.0) then
             write(99,*)' SLIP LIBRE',libre(k)
             nnppaarr=nnppaarr+1
           else
             write(99,*)' SLIP FIXE',libre(k)
           endif
c
          if(igeo.eq.2) then
            if(tapar(k,2).lt.0.0d0) then 
               kh=-1
            else
               kh=1
            endif
            call geoutm(tapar(k,1),tapar(k,2),tapar(k,1),tapar(k,2))
          endif
          call new2jcr(tapar(k,1),tapar(k,2),tapar(k,5),
     &    tapar(k,6),tapar(k,3),tapar(k,4))
c        write(*,*)(tapar(k,il),il=1,10) 
        enddo
        close(30)
          write(99,*)'nb de parametres libre = ',nnppaarr
          write(*,*)'nb de parametres libre = ',nnppaarr
         if(nnppaarr.eq.0) then
            print*,' ARRET PROGRAMME'
            write(99,*)' ARRET PROGRAMME'
            stop
         endif
c
        return
        end
c
c------------------------------------------------
c
        subroutine ecr_par(igeo,numx,numy)
c
       implicit real*8 (a-h,o-z)
c
        parameter (mxs=200,mxd=50000)
        parameter (mxf=9)
c
        character*100 fpar
       dimension tapar(mxs,10),libre(mxs)

c
        common/param/tapar,numsour,libre
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
        print*,'nom du fichier sortie parametre de faille'
        read(*,'(a)')fpar
        write(99,*)fpar
        open (40,file=fpar)
c
        nnuummss=0
        write(40,*)numsour,' 5  0  0 '
        do k=1,numsour
           nnuummss=nnuummss+1
           write(40,*)nnuummss
           call jcr2new(tapar(k,1),tapar(k,2),tapar(k,5),
     &    tapar(k,6),tapar(k,3),tapar(k,4))
           if(igeo.eq.2) then
             call utmgeo(tapar(k,1),tapar(k,2),tapar(k,1),tapar(k,2))
           endif
           write(99,*)'source no ',nnuummss
          write(40,*)tapar(k,1),zz,tapar(k,2),zz
          write(40,*)tapar(k,3),zz,tapar(k,4),zz
          write(40,*)tapar(k,5),zz,tapar(k,6),zz,tapar(k,7),zz
          write(40,*)tapar(k,8),zz,tapar(k,9),zz,tapar(k,10),zz
        enddo
c
c
        close(40)
        return
        end
c
c-----------------------------------------------
        subroutine lec_data(igeo,numx,numy)
c
       implicit real*8 (a-h,o-z)
c
         parameter (mxs=200,mxd=50000)
         parameter (mxf=9)
c
         character*100 ficdata,fivec
c
         dimension xdata(mxd),ydata(mxd),ddata(mxd)
         dimension weight(mxf),edata(mxd)
         dimension nv(mxf,2),irv(mxf)
       dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2)
c
         common/data1/xdata,ydata,ddata,edata,weight
        common/data2/nv,tsa,tsb,tsc,irv,dref
         common/data3/nfic
         COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
         amil=1000.d0
         nt=1
        do kk=1,nfic
          nv(kk,1)=nt
          print*,''
          print*,'nom du fichier no ',kk
          write(99,*)''
          write(99,*)'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
          write(99,*)'nom du fichier no ',kk
          read(*,'(a)')ficdata
          write(99,*)ficdata
c
          print*,'incertitude donnees (si lu dans fichier : 999)'
          read(*,*)erda
          write(99,*)'incertitude ',erda
          print*,'no du point de reference (=0 si absolu)'
          read(*,*)irv(kk)
          print*,'fichier de 4 vecteurs de projection'
          read(*,'(A)')fivec
          call planv(igeo,fivec,kk)
          print*,'poids sur les donnees'
          read(*,*) weight(kk)
          write(99,*)'poids ',weight(kk)
c
          open(31,file=ficdata)
          if(erda.eq.999.d0) then
 2          read(31,*,end=4) x,y,dam,erda
            if(igeo.eq.2) then
              call geoutm(x,y,x,y)
            endif
            xdata(nt)=x*amil
            ydata(nt)=y*amil
            ddata(nt)=dam
            edata(nt)=erda
            nt=nt+1
            goto 2
          else
 3          read(31,*,end=4) x,y,dam
            if(igeo.eq.2) then 
              call geoutm(x,y,x,y) 
            endif
            xdata(nt)=x*amil 
            ydata(nt)=y*amil 
            ddata(nt)=dam
            edata(nt)=erda
            nt=nt+1
            goto 3
          endif
 4        close(31)
c
          nv(kk,2)=nt-1
          print*,nv(kk,1),nv(kk,2)
c
c  difference par rapport a une reference eventuelle
c
          print*,'reference ?',irv(kk)
c
          if(irv(kk).ne.0) then
            nref=nv(kk,1)+irv(kk)-1
            dref(kk,1)=ddata(nref)
            print*,'point de reference'
            print*,'no, no, x , y, def '
      write(99,*)'mesures relatives au point : no, x y def'
      write(99,*)irv(kk),sngl(xdata(nref)),sngl(ydata(nref)),
     *sngl(dref(kk,1))
      call utmgeo(xdata(nref)/1000.,ydata(nref)/1000.,xdd,ydd)
      write(99,*)irv(kk),sngl(xdd),sngl(ydd),
     *sngl(dref(kk,1))
      print*,irv(kk),nref,xdata(nref),ydata(nref),dref(kk,1)
            do jk=nv(kk,1),nv(kk,2)
               ddata(jk)=ddata(jk)-dref(kk,1)
            enddo
          else
            print*,'mesures absolues'
            write(99,*)'mesures absolues'
            dref(kk,1)=0.0
          endif
            
        enddo
c
        return
          end
c
c--------------------------------------------
          subroutine planv(igeo,fivec,kk)
c--------------------------------------------
c
        IMPLICIT REAL*8(A-H,O-Z)
       parameter (mxf=9)
       character*100 fivec,nul*1
       dimension nv(mxf,2),irv(mxf)   
       dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2)
       dimension px(10),py(10),tx(10),ty(10),tz(10)
c
        common/data2/nv,tsa,tsb,tsc,irv,dref
         COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
        amil=1000.d0
c
        open(35,file=fivec)
          write(99,*)fivec
c
        read(35,'(a1)')nul
        read(35,*)nd
        do i=1,nd
 1        read(35,*)px(i),py(i),tx(i),ty(i),tz(i)
          write(99,*)sngl(px(i)),sngl(py(i)),sngl(tx(i)),
     #     sngl(ty(i)),sngl(tz(i))
          if(igeo.eq.2) then
          call geoutm(px(i),py(i),px(i),py(i))
          px(i)=px(i)*amil
          py(i)=py(i)*amil
          endif
        enddo
        close(35)
c
        sax=0.
        sbx=0.
        scx=0.
        say=0.
        sby=0.
        scy=0.
        saz=0.
        sbz=0.
        scz=0.
        nc=0
        do i=1,nd
          do j=i+1,nd
            do k=j+1,nd
c x
      bx=((tx(i)-tx(k))-(((tx(j)-tx(k))*(px(i)-px(k)))
     &/(px(j)-px(k))))/(((py(i)-py(k))*(px(j)-px(k))-(py(j)-py(k))*
     &(px(i)-px(k)))/(px(j)-px(k)))
          ax=((tx(i)-tx(k))-bx*(py(i)-py(k)))/(px(i)-px(k))
          cx=tx(i)-ax*px(i)-bx*py(i) 
          sax=sax+ax
          sbx=sbx+bx 
          scx=scx+cx 
c y
      by=((ty(i)-ty(k))-(((ty(j)-ty(k))*(px(i)-px(k))) 
     &/(px(j)-px(k))))/(((py(i)-py(k))*(px(j)-px(k))-(py(j)-py(k))* 
     &(px(i)-px(k)))/(px(j)-px(k))) 
          ay=((ty(i)-ty(k))-by*(py(i)-py(k)))/(px(i)-px(k)) 
          cy=ty(i)-ay*px(i)-by*py(i) 
          say=say+ay 
          sby=sby+by  
          scy=scy+cy
c z
      bz=((tz(i)-tz(k))-(((tz(j)-tz(k))*(px(i)-px(k)))  
     &/(px(j)-px(k))))/(((py(i)-py(k))*(px(j)-px(k))-(py(j)-py(k))* 
     &(px(i)-px(k)))/(px(j)-px(k)))  
          az=((tz(i)-tz(k))-bz*(py(i)-py(k)))/(px(i)-px(k))  
          cz=tz(i)-az*px(i)-bz*py(i)  
          saz=saz+az  
          sbz=sbz+bz   
          scz=scz+cz
c         print*,ax,bx,cx
c         print*,ay,by,cy 
c         print*,az,bz,cz 
          nc=nc+1
            enddo
          enddo
        enddo
        tsa(kk,1)=sax/nc
        tsb(kk,1)=sbx/nc
        tsc(kk,1)=scx/nc
        tsa(kk,2)=say/nc 
        tsb(kk,2)=sby/nc 
        tsc(kk,2)=scy/nc
        tsa(kk,3)=saz/nc 
        tsb(kk,3)=sbz/nc 
        tsc(kk,3)=scz/nc
         write(99,*)'plan a b c '
          do k=1,3
          write(99,*)sngl(tsa(kk,k)),sngl(tsb(kk,k)),sngl(tsc(kk,k))
          enddo
       return
       end

        
c-------------------------------------
          subroutine direct 
c
c
        IMPLICIT REAL*8(A-H,O-Z)
c
        parameter (mxs=200,mxd=50000)
       parameter (mxf=9)
c
        dimension xf(mxs),yf(mxs),ff(mxs),h(mxs)
        dimension w(mxs),r(mxs),tf(mxs)
        dimension us(mxs),ud(mxs),ut(mxs)
        dimension tapar(mxs,10) ,libre(mxs)
        dimension a(mxd,mxs)  
       dimension xdata(mxd),ydata(mxd),ddata(mxd)
       dimension weight(mxf),edata(mxd)
       dimension nv(mxf,2),irv(mxf)
       dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2)
c
        COMMON/F1/XF,YF,FF,H,W,R,TF
        COMMON/F2/US,UD,UT
        COMMON/COMOB/Q,ST,CT,ELAS,U0
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
        common/param/tapar,numsour,libre
        common/defo/a
        common/data1/xdata,ydata,ddata,edata,weight
        common/data2/nv,tsa,tsb,tsc,irv,dref
         common/data3/nfic 
c
        trad=pi/180.d0
        amil=1000.d0
       smu  = 3.3d10
c
c
        do i=1,numsour
             xf(i)=tapar(i,1)*amil
             yf(i)=tapar(i,2)*amil
             ff(i)=tapar(i,3)*trad
             tf(i)=tapar(i,4)*trad
             h(i)=tapar(i,5)*amil
             w(i)=tapar(i,6)*amil
             r(i)=tapar(i,7)*amil
             slip=tapar(i,8)
c
             usa=slip*dcosd(tapar(i,9))
             uda=-slip*dsind(tapar(i,9))
c            print*,usa,uda
             us(i)=usa
             ud(i)=uda
             ut(i)=0.0000001d0
c       write(*,*)(tapar(i,k),k=1,10) 
c       print*,'PARAMETRES',slip
c       print*,xf(i),yf(i),ff(i),tf(i),h(i),w(i),r(i),us(i),ud(i),ut(i)
        enddo
        call def_unit(numsour)
        return
        end
c
c-------------------------------------
          subroutine direct1
c
c
        IMPLICIT REAL*8(A-H,O-Z)
c
        parameter (mxs=200,mxd=50000)
       parameter (mxf=9)
c
        dimension xf(mxs),yf(mxs),ff(mxs),h(mxs)
        dimension w(mxs),r(mxs),tf(mxs)
        dimension us(mxs),ud(mxs),ut(mxs)
        dimension tapar(mxs,10),libre(mxs)
        dimension a(mxd,mxs)
       dimension xdata(mxd),ydata(mxd),ddata(mxd)
       dimension weight(mxf),edata(mxd)
       dimension nv(mxf,2),irv(mxf)
       dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2)
c
        COMMON/F1/XF,YF,FF,H,W,R,TF
        COMMON/F2/US,UD,UT
        COMMON/COMOB/Q,ST,CT,ELAS,U0
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
        common/param/tapar,numsour,libre
        common/defo/a
        common/data1/xdata,ydata,ddata,edata,weight
        common/data2/nv,tsa,tsb,tsc,irv,dref
         common/data3/nfic
c
        trad=pi/180.d0
        amil=1000.d0 
       smu  = 3.0d10
c
c
        do i=1,numsour
             xf(i)=tapar(i,1)*amil
             yf(i)=tapar(i,2)*amil
             ff(i)=tapar(i,3)*trad
             tf(i)=tapar(i,4)*trad
             h(i)=tapar(i,5)*amil
             w(i)=tapar(i,6)*amil
             r(i)=tapar(i,7)*amil
c            slip=1.d20/(w(i)*r(i)*2*smu)
             usa=(tapar(i,8))
             uda=(tapar(i,9))
             us(i)=usa
             ud(i)=uda
             ut(i)=0.0000001d0
        print*,'PAR',xf(i),yf(i),ff(i),tf(i),h(i),w(i),r(i)
        print*,usa,uda
c       write(*,*)(tapar(i,k),k=1,10)
        enddo
        call def_unit(numsour)
        return
        end
C--------------------------------------------------------------------------
      subroutine def_unit(numsour)
C-----------------------------------------
C           Calcul de la deformation due a un ensemble de failles
c
       IMPLICIT REAL*8(A-H,L,O-Z)
c
       parameter (mxs=200,mxd=50000)
       parameter (mxf=9)
       dimension xf(mxs),yf(mxs),ff(mxs),h(mxs)
       dimension w(mxs),r(mxs),tf(mxs)
       dimension xdata(mxd),ydata(mxd),ddata(mxd)
       dimension weight(mxf),edata(mxd)
       dimension us(mxs),ud(mxs),ut(mxs)
       dimension nv(mxf,2),irv(mxf)
       dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2)
       dimension a(mxd,mxs)
c
        COMMON/F1/XF,YF,FF,H,W,R,TF
        COMMON/F2/US,UD,UT
        common/defo/a
        common/data1/xdata,ydata,ddata,edata,weight
        common/data2/nv,tsa,tsb,tsc,irv,dref
         common/data3/nfic 
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
        COMMON/COMOB/Q,ST,CT,ELAS,U0
c
c
cc...initialize
c 
      do jj=1,nfic
       dref(jj,2)=0.d0 
      enddo
      do kk = 1,nv(nfic,2)
        do i = 1,numsour
          a(kk,i) = 0.
        enddo
      enddo
c
c boucle sur les failles
c
      DO 200 I=1,numsour
       FA=FF(I)
       TA=TF(I)
       HA=H(I)
      XPA=XF(I)-HA*DCOS(FA)/DTAN(TA)
      YPA=YF(I)+HA*DSIN(FA)/DTAN(TA)
       USA=US(I)
       UDA=UD(I)
       UTA=UT(I)
       WA = W(I)
       DA1=HA/DSIN(TA)
       DA2=DA1+WA
       RA=R(I) 
       HD= DA2*DSIN(TA)

C -------- boucle sur les points
c
      do 201 jj=1,nfic
      do 202 kk=nv(jj,1),nv(jj,2)
       dxunit=0.d0
       dyunit=0.d0
       dzunit=0.d0
       xx=xdata(kk)
       yy=ydata(kk)
C                 Coulissage :
      IF(DABS(USA).LE.0.001) GOTO 51
       UF=USA  
       IM=2 
       CALL FAULT( XX,YY,XPA,YPA,FA,TA,RA,WA,HD,UF,IM,A1,B1,C1)
         dxunit = dxunit+A1
         dyunit = dyunit+B1
         dzunit = dzunit+C1
c                  Plongeante :
   51 IF(DABS(UDA).LE.0.001) GOTO 52
       IM=3
       UF=-UDA
      CALL FAULT(XX,YY,XPA,YPA,FA,TA,RA,WA,HD,UF,IM,A1,B1,C1)
         dxunit = dxunit+A1
         dyunit = dyunit+B1
         dzunit = dzunit+C1
c                  Tensile :
   52 IF(UTA.LE.0.001) GOTO 53
       UF=UTA
       IM=1
      CALL FAULT(XX,YY,XPA,YPA,FA,TA,RA,WA,HD,UF,IM,A1,B1,C1)
         dxunit = dxunit+A1 
         dyunit = dyunit+B1 
         dzunit = dzunit+C1 
   53 CONTINUE
      tsx=tsa(jj,1)*xx+tsb(jj,1)*yy+tsc(jj,1)
      tsy=tsa(jj,2)*xx+tsb(jj,2)*yy+tsc(jj,2)
      tsz=tsa(jj,3)*xx+tsb(jj,3)*yy+tsc(jj,3)
      a(kk,i)=dxunit*tsx+dyunit*tsy+dzunit*tsz
c     print*,a(kk,i),ddata(kk)
c
 202  continue
c
c  difference a une reference eventuelle
c
      if(irv(jj).ne.0) then 
c       print*,'******** reference ********'
        nref=nv(jj,1)+irv(jj)-1
c       print*,nref
        ref=a(nref,i)
        dref(jj,2)=ref+dref(jj,2)
c       print*,'ref sommeref',ref,dref(jj,2)
        do kk=nv(jj,1),nv(jj,2)
          a(kk,i)=a(kk,i)-ref
        enddo
      else
       dref(jj,2)=0.0
      endif
 201  continue
c
 200  CONTINUE
      write(99,*)''
      write(99,*)'calcul direct'
      do jj=1,nfic
      write(99,*)'fichier no ',jj
      refd=(dref(jj,1))
      refc=(dref(jj,2))
      write(99,*)'no point,ref mesure, calcule, dif',irv(jj),sngl(refd),
     &sngl(refc),sngl(refc-refd)
      print*,'no point, ref mesure, calcule dif ',irv(jj),sngl(refd),
     &sngl(refc),sngl(refc-refd)
      enddo
c  
      RETURN
      END
c
cc----------------------------------------------------------------
C------------------------------------------------
      SUBROUTINE FAULT(XX,YY,XP,YP,FI,TE,AL,W,D,U,IMOD,A,B,C)
C-----------------------------
C
C          Calcul du deplacement dans le repere geographique
C          J.C.Ruegg-1986 d'apres Okada(1985,BSSA)
c
c                  Faille en tension  IMOD= 1
c                  Faille en coulissage  IMOD= 2
c                  Faille plongeante  IMOD= 3
C        Coordonnes initiales XX,YY dans repere geographique UTM
C        Repere local de faille,Y1
c        Repere Okada  X2,Y2
c        Sortie : Deplacements A,B,C surface libre
c
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COMOB/Q,ST,CT,ELAS,U0
c
      PI= 3.14159265359d0
      SF=DSIN(FI)
      CF=DCOS(FI)
      CT=DCOS(TE)
      ST=DSIN(TE)
C---------- Passage repere geographique au repere local X1,Y1
      X1= (XX-XP)*SF + (YY-YP)*CF
      Y1=-(XX-XP)*CF + (YY-YP)*SF
C---------- Passage au repere Okada
      X2= X1
      Y2= Y1 + D*CT/ST
C---------- Deformation dans le repere Okada
      QS1= X2+AL
      QS2= X2-AL
      P1= Y2*CT +D*ST
      P2= P1 - W
      Q= Y2*ST - D*CT
      U0 = U
c
C     WRITE(*,*) U0,CT,ST,CF,SF
c
      CALL DISLOK (QS1,P1,UX1,UY1,UZ1,IMOD)
      CALL DISLOK (QS1,P2,UX2,UY2,UZ2,IMOD)
      CALL DISLOK (QS2,P1,UX3,UY3,UZ3,IMOD)
      CALL DISLOK (QS2,P2,UX4,UY4,UZ4,IMOD)
c
      A2 = UX1-UX2-UX3+UX4
      B2 = UY1-UY2-UY3+UY4
      C2 = UZ1-UZ2-UZ3+UZ4
c
C--------- Deformation dans le repere geographique
c
      A = A2*SF - B2*CF
      B = A2*CF + B2*SF
      C = C2
      RETURN
      END
C--------------------------------------------------
      SUBROUTINE DISLOK(QSI,ETA,FUX,FUY,FUZ,IMOD)
C--------------------------------------------------
C     Calcul des deplacements ux,uy,uz
C           associes a une dislocation
C           FORMULES DE OKADA (1985, BSSA)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COMOB/Q,ST,CT,ELAS,U0
c
      PI= 3.14159265359d0
      PIP=PI+PI
      TT=ST/CT
c
      BY=ETA*CT+Q*ST
      BD=ETA*ST-Q*CT
      R= DSQRT(QSI**2+ETA**2+Q**2)
      XX=DSQRT(QSI**2+Q**2)
c
      RE= R+ETA
      RQ= R+QSI
      QRE= Q/R/RE
      QRQ= Q/R/RQ
c
      ATQE=DATAN(QSI*ETA/Q/R)
c
      AH= ETA*(XX+Q*CT)+XX*(R+XX)*ST
      BH= QSI*(R+XX)*CT
      TH= AH/BH
c
      AJ1= -ELAS*QSI/(R+BD)/CT
      AJ2= -ELAS*DLOG(R+ETA)
      AJ3=  ELAS*(BY/(R+BD)/CT-DLOG(R+ETA))
      AJ4=  ELAS*(DLOG(R+BD)-ST*DLOG(R+ETA))/CT
      AJ5=2.d0*ELAS*DATAN(TH)/CT
c
      IF (IMOD.NE.1) GOTO 10
C----------Faille en tension :
c
      FUX= Q*QRE - AJ3*ST*ST - AJ4*ST*ST*TT
      FUY=-BD*QRQ - ST*QRE*QSI + ST*ATQE - AJ1*ST*ST + AJ5*ST*ST*TT
      FUZ= BY*QRQ + CT*QSI*QRE - CT*ATQE - AJ5*ST*ST
      GOTO 20
c
   10 IF (IMOD.NE.2) GOTO 11
C------------Faille de coulissage pur (strike-slip) :
c
      FUX= -QSI*QRE - ATQE - AJ1*ST + AJ5*ST*TT
      FUY= -BY*QRE - Q*CT/RE - AJ2*ST + AJ3*ST + AJ4*ST*TT
      FUZ= -BD*QRE - Q*ST/RE - AJ4*ST
      GOTO 20
c
   11 IF(IMOD.NE.3) GOTO 12
C----------Faille plongeante pure ( dip-slip)
c
      FUX= -Q/R + AJ3*ST*CT + AJ4*ST*ST
      FUY= -BY*QRQ - CT*ATQE + AJ1*ST*CT - AJ5*ST*ST
      FUZ= -BD*QRQ - ST*ATQE + AJ5*ST*CT
      GOTO 20
c
  20  FUX= FUX*U0/PIP
      FUY= FUY*U0/PIP
      FUZ= FUZ*U0/PIP
c
  12  RETURN
      END
c
c--------------------------------------------
C ------------------------------------------------
      SUBROUTINE UTMGEO(XE,YE,xLON,yLAT)
C --------------------------------------
C         TRANSFORMATION  UTM - GEOGRAPHIQUE
C           ENTREE :  COORD. UTM  EN  KILOMETRES
C           SORTIE :  COORD. GEOGRAPH. EN DEGRES
C
C           HEMISPHERE NORD : KH=1
C                      SUD  : KH=-1
C
      IMPLICIT REAL*8 (A-H,L,O-Z)
      COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
c      print*,'utmgeo i ',xe,ye


      XX=XE*1000.
      YY=YE*1000.
       MI=IFU*6-183
      AMI=FLOAT(MI)/0.9D0
      RD=PI/200.D0
      AM0=AMI*RD 
      YR=YY
      if(KH.EQ.1) goto 11
      YR=10.D6-YR
   11 YR=YR/AK  
      P=KH*YR*PI*0.5D-7
      CO=DCOS(P) 
      zL=P*KH  
      S=E2*CO*CO
      R=CC/DSQRT(1.d0+S)
      PP=P*KH    
      YR=YR-ARCME(zL)
      U=(XX-5.D5)/AK/R
      V=S*U*U
      Q=YR/R*(1.D0-V/2.D0)
      U=U*(1.D0-V/6.D0)
      P=zL+Q
      W=DEXP(U)
      AM=DATAN((W*W-1.D0)/(2.D0*W*DCOS(P)))
      V=DSIN(P)/DCOS(P)
      G=DATAN(V*COS(AM))
      V=1.D0+S-1.5D0*DSIN(zL)*E2*(G-zL)*DCOS(zL)
      xLON=AM+AM0
      yLAT=zL+V*(G-zL)
      xLON=(AM+AM0)/RD
      yLAT=yLAT/RD*KH
      xlon=xlon*0.9d0
      ylat=ylat*0.9d0
     

c      print*,'utmgeo o ',xlon,ylat
c
      RETURN
      END
C ------------------------------------------------
C ------------------------------------------------
       FUNCTION ARCME(AL)
C--------------------------
C         Fonction Arc de m\002ridien en fonction de la latitude :
 
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COMOG/CC,E2,PI,IFU,AK,KH
 
      AQ=E2*3.d0/4.d0
      BQ=E2*E2*15.d0/16.d0
      CQ=(E2**3)*35.d0/64.d0
      DQ=(E2**4)*315.d0/512.d0
      U=DSIN(AL)*DCOS(AL)
      V=DCOS(AL)**2
      AJ2=AL+U
      AJ4=(AJ2*3.d0+U*V*2.d0)/4.d0
      AJ6=(AJ4*5.d0+2.d0*U*V**2)/3.d0
      AJ8=(AJ6*7.d0+4.d0*U*V**3)/8.d0
      ARCME=CC*(AL-AQ*AJ2+BQ*AJ4-CQ*AJ6+DQ*AJ8)
      RETURN
      END
c
c-----------------------------------------
C-----------------------------------------------------------
      SUBROUTINE GEOUTM(LON,LAT,XX,YY)
C---------------------------------------------
C          TRANSFORMATION GEOGRAPHIQUE - U.T.M.
C            ENTREE : LONGITUDE ET LATITUDE EN degres
C            SORTIE :  COORD. UTM  EN KILOMETRES
C
      IMPLICIT REAL*8 (A-H,L,O-Z)
      COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c-----------------------------------
c      PI=4.d0*DATAN(1.d0)
       eps = 0.00000001
       trg=200.d0/180.d0
c      CC= 6399593.6257585d0
c      EPRIM=  0.08209443794984d0
c      E2=EPRIM**2
c      IFU= 54
       AK= 0.9996d0
c ----------------------------------
c      print*,'geoutm i ',lon,lat
      lo=lon*trg
      la=lat*trg
      MI=IFU*6-183
      AMI=DBLE(MI)/0.9D0
      RD=PI/200.D0
      AM=LO*RD
      AL=LA*RD
      AM0=AMI*RD
      CO=DCOS(AL)
      AMU=AM-AM0
      XI=CO*DSIN(AMU)
      XI=0.5D0*DLOG((1.D0+XI)/(1.D0-XI))
c      TA=DATAN(DSIN(AL)/(CO*DCOS(AMU)))-AL
      TA=DATAN2(DSIN(AL),CO*DCOS(AMU))-AL
      GN=CC/DSQRT(1.D0+E2*CO*CO)
      CX=E2*(CO*XI)**2
      DX=GN*XI*(1.D0+CX/6.D0)
      DY=GN*TA*(1.D0+CX/2.D0)
      XX=500000.D0+AK*DX
      YY=AK*(DY+ARCME(AL))
      if(LAT) 1,2,2
    1 YY=(YY+10000000.d0)
    2 CONTINUE
      xx=xx/1000.
      yy=yy/1000.
c      print*,'geoutm o ',xx,yy
      RETURN
      END
c-----------------------------------------
c--------------------------------------------
c      subroutine de transformation de :
c          us,ud en slip, rake
c
        subroutine usd2sr (slip,rake)
c
c       us>0 senestre
c       us<0 dextre
c       ud>0 normal
c       ud<0 inverse
c
c       -180 <rake< 180
c       rake>0 inverse
c       rake<0 normal
c       90< |rake| <180 senestre
c       0 < |rake| <90 dextre
c
      IMPLICIT REAL*8(A-H,O-Z)
c
 1      continue
        us=slip
        ud=rake
c
        slip=dsqrt(us**2+ud**2)
        rake=-datan2d(ud,us)
c
        return
        end
c
c-----------------------------------------
c
c      subroutine de transformation de :
c          slip,rake  ---> en us ud
c
        subroutine sr2usd (slip,rake)
c----------------------------------------
c
c       us>0 senestre
c       us<0 dextre
c       ud>0 normal
c       ud<0 inverse
c
c       -180 <rake< 180
c       rake>0 inverse
c       rake<0 normal
c       90< |rake| <180 senestre
c       0 < |rake| <90 dextre
c
       IMPLICIT REAL*8(A-H,O-Z)
c
 1      continue
c
        us=slip*dcosd(rake)
        ud=-slip*dsind(rake)
c
        slip=us
        rake=ud
c
        return
        end
c
c--------------------------------------------
c
        subroutine jcr2new (xs,ys,h,w,str,dip)
c
c   transformation de xs ys h (prof) w (largeur)
c  en xq yq (en surface) h1 h2 (deux profondeurs de la faille)
c
      IMPLICIT REAL*8(A-H,O-Z)
c
        wh=h/dtand(dip)
        xq=xs-wh*dcosd(str)
        yq=ys+wh*dsind(str)
        h2=w*dsind(dip)+h
c
        xs=xq
        ys=yq
        w=h2
c
        return
        end
c
c
c-----------------------------------------
c
        subroutine new2jcr (xq,yq,h1,h2,str,dip)
c
c   transformation de xq yq (en surface) h1 h2 (deux profondeurs de la faille)
c   en  xs ys h (prof) w (largeur)
c
      IMPLICIT REAL*8(A-H,O-Z)
c
      wh=h1/dtand(dip)
      xs=xq+wh*dcosd(str)
      ys=yq-wh*dsind(str)
      wa=(h2-h1)/dsind(dip)
c      
      xq=xs
      yq=ys
      h1=h1
      h2=wa
      return
      end
c      
c----------------------------------------
