PROGRAM INV_GEN_VEC
c
c
C     PROGRAMME D'INVERSION DE DONNEES DE DEFORMATION DU SOL
C     PIERRE BRIOLE
C     30/8/87
C     Sur sun le 11 Mars 1992
c     modifie 1997 par JBC 
c
C ------------------------------------------------
C     CE PROGRAMME UTILISE LES VARIABLES SUIVANTES
C ------------------------------------------------
c
c    possibilite de travailler en UTM ou en geographique (degres)
c
C ****************************************
C    UN FICHIER DE PARAMETRE DE FAILLES avec les incertitude sur chaque parametre (si incertitude=0 on n'inverse pas le parametre)
C ****************************************
c      2 format possibles
c
c*** format type 1 : invers3 :
c       NF,NINV, 0, 0  : nb de failles, nb d'inversions (maxi 5), 0 , 0
C       XO,YO: PROJECTION EN SURFACE DU CENTRE DE L'ARETE SUPERIEURE DE
C              LA FAILLE (en m)
C       PHI: AZIMUT DE LA TRACE EN SURFACE DE LA FAILLE COMPTE POSITIF
C            DANS LE SENS HORAIRE A PARTIR DU NORD (en deg)
C       H: PROFONDEUR DU CENTRE DE L'ARETE SUPERIEURE DE LA FAILLE (en m)
C       L,D: DEMI-LONGUEUR DE LA FAILLE, LARGEUR DE LA FAILLE (en m)
C       TETA: PENDAGE DE LA FAILLE (en deg)
c       SLIP, RAKE, OUV : MODULE DU GLISSEMENT (cm), RAKE (deg), ouverture (cm)
c             slip toujours positif, 
c             rake <0 failles normales, >0 pour inverses
c             -90<rake<90 pour failles senestres 
c            -180<=rake<-90 et 90<rake<=180 pour dextres
c            purement dextre rake = +-180
c            purement senestre : rake= 0 
c
c**** format de type 2 : jbc :
c       NF,NINV, 0, 0  	:nb de failles, nb d'inversions (maxi 5), 0 , 0 
c       N0         	:numero de la faille
c      XS, YS 		:intersection avec la surface du centre de l'arrete superieure de la faille (km)
c      PHI, TETA 	:azimut et dip de la faille (deg)
c      H1, H2, L, 	:profondeur minimum et maximum de la dislocation, demi longueur de la faille (km)
c      SLIP, RAKE, OUV 	:Idem que le precedent format (cm, deg, cm)
c
C ************************************************
C    DONNEES : 
C ************************************************
c      sont lues dans les fichiers de donnees 
c      on peut lire les incertitude sur les donnees en ligne (meme incertitude pour toutes les donnees d'un meme fichier)
c     ou alors dans le fichier (si l'incertitude = 999)
c
c    X Y DONNEE (INCERTITUDE_DONNEE)   : unite : les positions doivent arriver en km (ou degres), et les donnees en cm
c                                        on donne un facteur multiplicatif si les donnees ne sont pas de la bonne unite
c
c    incertitude sur les donnees : meme unite que les donnees ( si = 999, on lit les incertitudes en 4e colonne du fichier)
c    reference :   0 ou 1 si ref=0 on considere les mesures absolues, si ref=1 les mesures sont relatives au premier point
c                     du fichier
c    vecteur unitaire vx vy vz : sur lequel sera projete le vecteur deplacement (ex : pour donnnees verticales 0 0 1)
c    
c
C ------------------------------------------------------
C     CE PROGRAMME CONTIENT LES SOUS-PROGRAMMES SUIVANTS
C ------------------------------------------------------
c
C     INI: PROGRAMME D'ENTREE DES DONNEES
C     SORT: SORTIE SUR FICHIER DES RESULTATS
C     DISLO: CALCULE (DEPL) LES DEPLACEMENTS DES POINTS (POINT) AVEC DES
C            PARAMETRES DE FAILLE DONNES (FAILLE)
C     FAULT: SOUS-PROGRAMME APPELE PAR DISLO
C     DISLOK: SOUS-PROGRAMME APPELE PAR FAULT
C     INVERSE: INVERSE LES PARAMETRES DES FAILLES PAR LA METHODE DE
C              TARANTOLA-VALETTE (FAILLE) ET CALCULE LES DEPLACEMENTS
C              CORRESPONDANTS (DEPL)
C     CHOL: SOUS-PROGRAMME D'INVERSION DE MATRICE
c
C ----------------------
C     REMARQUES DIVERSES
C ----------------------
c
C      LE NOMBRE MAXIMUM DE POINTS : NPM
c                          FAILLES : NFM
c                        INVERSIONS : NINVM (NINVM-1)
C      L'INVERSION DES PARAMETRES EST EFFECTUEE UNIQUEMENT SUR CEUX POUR
C      LESQUELS ERPA EST DIFFERENT DE 0, LES AUTRES PARAMETRES SONT
C      CONSIDERES COMME CONSTANTS, LE PROGRAMME NE FONCTIONNE PAS SI
C      TOUS LES ERPA SONT A 0, MEME POUR LE CALCUL DIRECT
c
        IMPLICIT REAL*8 (A-H,O-Z)
        IMPLICIT REAL*8 (M)
        PARAMETER (NPM=25000,NFM=50,NFMX9=400,NINVM=15,nfi=9)
        CHARACTER*80 POINTV(nfi),FICHV
        COMMON /DIS/ POINT(NPM,2),DEPL(NPM,2),FAILLE(NFM,10),NV(nfi,2),
     &  IRV(nfi),tsa(nfi,3),tsb(nfi,3),tsc(nfi,3),dref(nfi,2)
        COMMON /INV/ ERPA(NFM,10),ERD(NPM),MATV(NFMX9,NFMX9),NPAR,IV,IMP
        COMMON /SOR/ DEPLS(NINVM,NPM),FAILS(NINVM,NFMX9)
        COMMON /ES/ POINTV,FICHV
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
        COMMON/TYP/ityp,igeo,facs,facd
c
c--- changement coordonnees UTM GEO
c
        PI= 3.14159265359
        cc=6399593.6258000
        AK=0.9996d0
        eprim=0.08209443779496
        EP2=EPRIM*EPRIM
        E2=EP2
        ifu=19
C ----------Initialisation
        open(99,file='INVGEN.LOG')
        CALL INI(NPOINTS,NF,NINV,nff)
c       do k=1,nf
c          write(*,*)(faille(k,i),i=1,10)
c       enddo
c       do k=1,npoints
c          print*,depl(k,1)
c       enddo
c
C ----------Boucle d'inversion (NINV=0 pour le calcul direct)
      DO 1, I=1,NINV+1
        IF(I.EQ.1) THEN
          
          CALL DISLO(NPOINTS,1,NF,nff)
        ELSE
          WRITE(*,*) 'Inversion numero ',I-1
          write(99,*)'Inversion numero ',I-1
          CALL INVERSE(NPOINTS,NF,nff)
C ----------Ouvre un fichier pour MATV lorsque IV different de 0
          IF(IV.NE.0.AND.I.EQ.NINV+1) THEN
C ----------Ne garde que le dernier tableau
c ?????? 
c         READ(16,'(A)')  FICHV
c         print*,'***** 16 ********'
c         OPEN(9,STATUS='NEW',FILE=FICHV)
c         DO 3 J=1,NPAR
c   3     WRITE(9,11) (MATV(J,L),L=1,NPAR)
   11   FORMAT(20F10.5)
c         CLOSE(9)
          ENDIF
        ENDIF

C ----------Remplissage du tableau DEPLS
        DO 2 K=1,NPOINTS
        DEPLS(I,K)=DEPL(K,2)
  2     continue
C ----------Remplissage du tableau FAILS
        NUM=1
        DO 5 K=1,NF
        DO 5 L=1,10
        FAILS(I,NUM)=FAILLE(K,L)
    5   NUM=NUM+1

    1 CONTINUE
       CALL SORT(NINV,NPOINTS,NF,nff)

      close(99)
      STOP
      END

C -----------------------------------
C     SOUS-PROGRAMME D'INITIALISATION
C -----------------------------------

      SUBROUTINE INI(K,NFT,NINV,nff)

      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT REAL*8 (M)
      PARAMETER (NPM=25000,NFM=50,NFMX9=400,NINVM=15,nfi=9)
      CHARACTER*80 POINTV(nfi),FICHV,PARAM,fivec
      logical stopp
        COMMON /DIS/ POINT(NPM,2),DEPL(NPM,2),FAILLE(NFM,10),NV(nfi,2),
     &  IRV(nfi),tsa(nfi,3),tsb(nfi,3),tsc(nfi,3),dref(nfi,2)
      COMMON /INV/ ERPA(NFM,10),ERD(NPM),MATV(NFMX9,NFMX9),NPAR,IV,IMP
      COMMON /ES/ POINTV,FICHV
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
        COMMON/TYP/ityp,igeo,facs,facd
c
        print*,'type de parametrisation en entree'
        print*,'1 ---> type invers3 '
        print*,'2 ---> type jbc' 
        read(*,*)ityp
        print*,'ityp',ityp
c
        print*,'coordonnees et donnees en ----> UTM : 1'
        print*,'                        geographique :2'
        read(*,*)igeo
c
        if(igeo.eq.2) then
          print*,'fuseau utm, et hemisphere (1:nord, -1:sud)'
          read(*,*)ifu,kh
          print*,'coordonnees en geographiques'
          print*,'fuseau hemisphere',ifu,kh
          write(99,*)'geographique, fuseau, hemisphere',ifu,kh
        else
          print*,'coordonnees en UTM'
          write(99,*)'coordonnees en UTM'
        endif
c
C ----------Lecture des parametres sur fichier
      WRITE(*,*) 'Nom du fichier de parametres'
          write(99,*)'Nom du fichier de parametres'
      READ(*,'(A)') PARAM
      print*,param
          write(99,*)param
      OPEN(16,STATUS='OLD',FILE=PARAM)

C ----------Lecture du nombre de failles et d'inversions
C           et de IV (pour avoir les variances a-posteriori)
      READ(16,*) NFT,NINV,IV,IMP
          write(99,*)'nb failles, inversions ',nft,ninv
c
C ----------Introduction des 10 parametres de chaque faille dans
C           l'ordre indique dans les formules d'introduction et sous la
C           forme "P EP" (Parametre et ecart type)
c
        if(ityp.eq.1) then
          DO 1, NF=1,NFT
           write(99,*)'faille no',nf
           READ(16,*) (FAILLE(NF,K),ERPA(NF,K),K=1,10)
           write(99,*) (sngl(FAILLE(NF,K)),sngl(ERPA(NF,K)),K=1,10)
           if(igeo.eq.2) then
      call geoutm(faille(nf,1),faille(nf,2),faille(nf,1),faille(nf,2))
           endif
    1     CONTINUE
        else
          do 5 NF=1,NFT
           read(16,*)innu
           read(16,*)faille(nf,1),erpa(nf,1),faille(nf,2),erpa(nf,2)
           read(16,*)faille(nf,3),erpa(nf,3),faille(nf,7),erpa(nf,7)
           read(16,*)faille(nf,4),erpa(nf,4),faille(nf,6),erpa(nf,6)
     &              ,faille(nf,5),erpa(nf,5) 
           read(16,*)faille(nf,8),erpa(nf,8),faille(nf,9),erpa(nf,9)
     &              ,faille(nf,10),erpa(nf,10) 
           write(99,*) (sngl(FAILLE(NF,K)),sngl(ERPA(NF,K)),K=1,10)
           if(igeo.eq.2) then
        call geoutm(faille(nf,1),faille(nf,2),faille(nf,1),faille(nf,2)) 
           erpa(nf,1)=erpa(nf,1)*1000.*100.
           erpa(nf,2)=erpa(nf,2)*1000.*100.
           else
           faille(nf,1)=faille(nf,1)*1000. 
           faille(nf,2)=faille(nf,2)*1000.
           erpa(nf,1)=erpa(nf,1)*1000. 
           erpa(nf,2)=erpa(nf,2)*1000.
           endif
           faille(nf,4)=faille(nf,4)*1000. 
           faille(nf,6)=faille(nf,6)*1000.
           faille(nf,5)=faille(nf,5)*1000.
           erpa(nf,4)=erpa(nf,4)*1000.
           erpa(nf,5)=erpa(nf,5)*1000. 
           erpa(nf,6)=erpa(nf,6)*1000. 
    5     continue
        endif
        close (16)
c         do nf=1,nft
c         write(*,*)(faille(nf,i),i=1,10)
c         write(*,*)(erpa(nf,i),i=1,10)   
c         enddo
c
C ----------Intoduction des fichiers de donnees
C           Les donnees sont contenues dans des fichiers differents
C           selon leur nature (trois natures differentes possibles)
C           Le programme demande successivement le nom des trois
C           fichiers et l'incertitude sur les donnees associee (on
C           suppose que des donnees de meme nature ont la meme
C           incertitude)
C           Si pas de donnees d'un type appeler le fichier NULL.PTS
c
c   lecture du nombre de fichier de donnees
      print*,'nb de fichiers de donnees'
      read(*,*)nff
      print*,' '
      print*,nff,' fichier(s)'
      write(99,*)' '
      write(99,*)nff,' fichier(s)'
c
      K=1
c boucle sur les fichiers
       do 2 i=1,nff
       nv(i,1)=k
c
C ----------Lecture du nom du fichier
        print*,' '
        print*,'****fichier no ',i
        write(99,*)' '
        write(99,*)'****fichier no ',i
        READ(*,'(A)') POINTV(I)
        print*,POINTV(I)
        write(99,*)POINTV(I)

C ----------Lecture des incertitudes sur les donnees
C           Si les incertitudes sur les donnees varient pour un meme
C           type de donnees, il faut modifier le programme et les
C           introduire dans le fichier des points
      print*,'facteurs multiplicatifs '
      print*,'pour la position des points et les donnees' 
      print*,' doivent etre en km (ou deg) et cm (1 1)'
      read(*,*)facs,facd
      print*,'****facteurs****',facs,facd
        write(99,*)'****facteurs****',facs,facd
      if(igeo.eq.1) then
         facs=facs*1000d0
      else
         facs=1d0
      endif
      print*,'incertitudes sur les donnees (si lu dans fichier:999)'
      read(*,*)erda
      print*,'****incertitude :****',erda
        write(99,*)'****incertitude :****',erda
      print*, 'reference : '
      print*,'si 0 donnees absolues'
      print*,'si n (different de 0) reference au nieme point du fichier'
      read(*,*)IRV(i)
      if(irv(i).eq.0) then
        print*,'****donnees absolues****'
      else
        print*,'***donnees relatives***'
        print*,' au ',irv(i),'eme point du fichier'
      endif
c     print*,'vecteur unitaire x y z'
c     read(*,*)ts(i,1),ts(i,2),ts(i,3)
c     print*,'****vecteur***',ts(i,1),ts(i,2),ts(i,3)
          print*,'fichier de 4 vecteurs de projection'
          read(*,'(A)')fivec
          call planv(igeo,fivec,i)
          do ia=1,3
            print*,tsa(i,ia),tsb(i,ia),tsc(i,ia)
          enddo
c 
C ----------Lecture des points sur fichier
      OPEN(7,STATUS='OLD',FILE=POINTV(I))
c
c
c       if(igeo.eq.1) then
        if(erda.eq.999.d0) then
         DO 3, N=0,NPM
            READ(7,*,END=4) X,Y,DAM,erda
            POINT(K,1)=X*facs
            POINT(K,2)=Y*facs
            DEPL(K,1)=DAM*facd
            ERD(K)=ERDA*facd
            K=K+1
    3    CONTINUE
        else
         DO 6, N=0,NPM
            READ(7,*,END=4) X,Y,DAM
            POINT(K,1)=X*facs
            POINT(K,2)=Y*facs
            DEPL(K,1)=DAM*facd
            ERD(K)=ERDA*facd
            K=K+1
    6    CONTINUE
        endif
c
    4 CLOSE(7)
      NV(I,2)=k-1
      if(igeo.eq.2) then
        do ig=nv(i,1),nv(i,2)
        call geoutm(point(ig,1),point(ig,2),point(ig,1),point(ig,2))
        enddo
      endif
c-- reference eventuelle au nieme point de la serie irv()
        IF(IRV(i).NE.0) THEN
          irel=nv(i,1)+irv(i)-1
          rell=DEPL(irel,1)
          dref(i,1)=DEPL(irel,1)
          write(99,*)'donnees relative ',irv(i),sngl(dref(i,1))
          DO  kk=nv(i,1),nv(i,2)
            DEPL(kk,1)=DEPL(kk,1)-dref(i,1)
          enddo
        else
            print*,'mesures absolues'
          write(99,*)'donnees absolues'
            dref(i,1)=0.0
        ENDIF
      print*,' '
    2 CONTINUE
      K=K-1

C ----------Controle des dimensions
      nmes=k-1
      stopp=.false.
      IF(NFT.GT.nfm) then
        WRITE(*,*) 'Trop de failles ',NFT
        print*,'MAXI= ',nfm
        stopp=.true.
      endif
      if(NINV.GT.ninvm) THEN
        WRITE(*,*) 'Trop d"inversions :',ninv
        print*,'maxi=',ninvm
        stopp=.true.
      endif
      IF(NMES.GT.npm) THEN
        WRITE(*,*) 'Trop de points',nmes
        print*,'maxi= ',npm
        STOPp=.true.
      ENDIF
      if (stopp) then
       print*,'changer les dimensions'
       print*,'ARRET'
       stop
      endif
      print*,' '
       print*,'--------------------------'
       print*,nff,'  fichiers de donnees'
      WRITE(*,*) K,'    Points lus'
      print*,' '
      do if=1,nff
        print*,pointv(if)
        print*,nv(if,2)-nv(if,1)+1,' points'
        print*,' '
      enddo
      WRITE(*,*) NFT,'    Failles'
      WRITE(*,*) NINV,'    Inversions'
       print*,'--------------------------'
c
      RETURN
      END

C ------------------------------
C     SOUS-PROGRAMME D'INVERSION
C ------------------------------

      SUBROUTINE INVERSE(NPOINTS,NF,nff)

      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT REAL*8 (M)
      INTEGER M
      PARAMETER (NPM=25000,NFM=50,NFMX9=400,NPMX2=NPM*2,nfi=9)
        COMMON /DIS/ POINT(NPM,2),DEPL(NPM,2),FAILLE(NFM,10),NV(nfi,2),
     &  IRV(nfi),tsa(nfi,3),tsb(nfi,3),tsc(nfi,3),dref(nfi,2)
      COMMON /INV/ ERPA(NFM,10),ERD(NPM),MATV(NFMX9,NFMX9),NPAR,IV,IMP
        COMMON/TYP/ityp,igeo,facs,facd

C ----------Matrice des derivees partielles, matrice MAT a inverser et
C           vecteurs
      DIMENSION MATRICE(NPMX2,NFMX9),MAT(NFMX9,NFMX9),X(NFMX9)
c      NMES=NV(1)+2*(NV(2)+NV(3))
       nmes=nv(nff,2)

C ----------Vidage des matrices
      DO 11, I=1,NFMX9
      DO 11, J=1,NFMX9
      MAT(I,J)=0
   11 CONTINUE
      DO 12, I=1,NMES
      DO 12, J=1,NFMX9
      MATRICE(I,J)=0
   12 CONTINUE

C ----------Boucle realisant:
C             remplissage de la matrice
C             remplissage de la premiere partie de MAT (termes lies aux
C             incertitudes sur les parametres)
C             calcul du nombre (NPAR) de parametres a inverser
      NPAR=1
      DO 1, I=1,NF
      DO 1, J=1,10

C ----------Ne calcule pas lorsque ERPA=0 (parametre constant)
        IF(ERPA(I,J).NE.0) THEN
          FAILLE(I,J)=FAILLE(I,J)+0.01
          CALL DISLO(NPOINTS,I,I,nff)
          KB=1

          DO 3, K=1,NPOINTS
            MATRICE(KB,NPAR)=DEPL(K,2)
            KB=KB+1
    3     CONTINUE

          FAILLE(I,J)=FAILLE(I,J)-0.02
          CALL DISLO(NPOINTS,I,I,nff)
          KB=1

          DO 4, K=1,NPOINTS
          MATRICE(KB,NPAR)=(MATRICE(KB,NPAR)-DEPL(K,2))*100
          KB=KB+1
    4     CONTINUE

          FAILLE(I,J)=FAILLE(I,J)+0.01
          MAT(NPAR,NPAR)=1/ERPA(I,J)**2
          NPAR=NPAR+1
        ENDIF

    1 CONTINUE
      NPAR=NPAR-1

C ----------Remplissage de la seconde partie de MAT (termes lies aux
C           derivees partielles)
      DO 5, I=1,NPAR
      DO 5, L=1,NPAR
      K=1
      DO 5, J=1,NPOINTS
      MAT(I,L)=MAT(I,L)+MATRICE(K,I)*MATRICE(K,L)/(ERD(J)**2)
      K=K+1
    5 CONTINUE

C ----------Calcul facultatif des variances a-posteriori
      IF(IV.NE.0) THEN
      DO 9, J=1,NPAR
      DO 10, I=1,NPAR
      X(I)=0.
   10 CONTINUE
      X(J)=1.
      CALL CHOL(MAT,X,NPAR)
      DO 8, L=1,NPAR
      MATV(L,J)=X(L)
    8 CONTINUE
    9 CONTINUE
      ENDIF
C ----------Calcule les valeurs theoriques des deplacements utilises
C           par l'inversion
      CALL DISLO(NPOINTS,1,NF,nff)

C ----------Calcule le vecteur d'entree
      DO 2, I=1,NPAR
      X(I)=0
      K=1
      DO 2, J=1,NPOINTS
      X(I)=X(I)+MATRICE(K,I)*(DEPL(J,1)-DEPL(J,2))/(ERD(J)**2)
      K=K+1
    2 CONTINUE

C ----------Calcul du vecteur des variations
      CALL CHOL(MAT,X,NPAR)

C ----------Changement des parametres des failles
      NPAR=1
      DO 17, I=1,NF
      DO 17, J=1,10
      IF(ERPA(I,J).NE.0) THEN
      AQUO=faille(i,j)/(faille(i,j)+x(npar))
      if(dabs(aquo).ne.aquo) then
        faille(i,j)=dabs(faille(i,j))/faille(i,j)*0.0001
      else
      FAILLE(I,J)=FAILLE(I,J)+X(NPAR)
      endif
C ----------Changement facultatif des incertitudes sur les parametres
      IF(IV.NE.0.AND.IMP.NE.0) THEN
      ERPA(I,J)=DSQRT(MATV(NPAR,NPAR))
      ENDIF
      NPAR=NPAR+1
      ENDIF
   17 CONTINUE
      NPAR=NPAR-1

C ----------Controle des nouvelles valeurs
      DO 18 I=1,NF

      IF(FAILLE(I,3).GT.180.) THEN
      FAILLE(I,3)=FAILLE(I,3)-360.
      ENDIF
      IF(FAILLE(I,3).LT.-180.) THEN
      FAILLE(I,3)=FAILLE(I,3)+360.
      ENDIF
c
c cas type 2 jbc controle que h2>=h1
      if(faille(i,6).lt.faille(i,4).and.ityp.eq.2) then
       faille(i,6)=faille(i,4)+0.01
      endif

      DO 20 J=4,6
      IF(FAILLE(I,J).LT.0.) THEN
      WRITE(*,*) 'Faille nø',I
      IF(J.EQ.4) THEN
      WRITE(*,*) 'Probleme: H negatif - Remis a zero arbitrairement'
      ELSE IF(J.EQ.5) THEN
      WRITE(*,*) 'Probleme: L negatif - Remis a zero arbitrairement'
      ELSE
      WRITE(*,*) 'Probleme: D negatif - Remis a zero arbitrairement'
      ENDIF
      FAILLE(I,J)=0.
      ENDIF
   20 CONTINUE

      FI7=FAILLE(I,7)
      IF(FI7.LT.0.OR.FI7.GT.90.) THEN
      FI9=FAILLE(I,9)
      FI3=FAILLE(I,3)
      FAILLE(I,9)=-FI9
      IF(FAILLE(I,3).GT.0.) THEN
      FAILLE(I,3)=FAILLE(I,3)-180.
      ELSE
      FAILLE(I,3)=FAILLE(I,3)+180.
      ENDIF
      FAILLE(I,7)=180.-FI7
      IF(FI7.LT.0.) THEN
      FAILLE(I,7)=-FI7
      FAILLE(I,4)=FAILLE(I,4)-FAILLE(I,6)*DSIN(FI7)
      DLL=FAILLE(I,6)*DCOS(FI7)
      FAILLE(I,1)=FAILLE(I,1)+DLL*DCOS(FI3)
      FAILLE(I,2)=FAILLE(I,2)-DLL*DSIN(FI3)
      IF(FAILLE(I,4).LT.0.) THEN
      FAILLE(I,4)=0.
      ENDIF
      ENDIF
      ENDIF
      IF(FAILLE(I,10).LT.0.) THEN
      WRITE(*,*) 'Faille nø',I
      WRITE(*,*) 'Probleme: OUV negatif - Remis a zero arbitrairement'
      FAILLE(I,10)=0.
      ENDIF
   18 CONTINUE

C ----------Calcul des deplacements pour la nouvelle faille
      CALL DISLO(NPOINTS,1,NF,nff)

      RETURN
      END

C ------------------------------------------------------------------
C     SOUS-PROGRAMME D'INVERSION DE MATRICE (CALCUL DE X AVEX Y=A*X)
C ------------------------------------------------------------------

      SUBROUTINE CHOL(MAT,X,NPAR)

      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT REAL*8 (M)
      PARAMETER (NFMX9=400)
      DIMENSION MAT(NFMX9,NFMX9),X(NFMX9),DIAG(NFMX9)

C ----------Inversion de la matrice MAT par la methode de CHOLEWSKI
C           Le vecteur de sortie est place dans X(I)
      DIAG(1)=DSQRT(MAT(1,1))
      DO 6 I=2,NPAR
    6 MAT(1,I)=MAT(I,1)/DIAG(1)
      DO 10 J=2,NPAR
      S=MAT(J,J)
      KMX=J-1
      DO 7 K=1,KMX
    7 S=S-MAT(K,J)**2
      DIAG(J)=DSQRT(S)
      IF(J.EQ.NPAR) GOTO 10
      IMI=J+1
      DO 8 I=IMI,NPAR
      S=MAT(I,J)
      DO 9 K=1,KMX
    9 S=S-MAT(K,I)*MAT(K,J)
      MAT(J,I)=S/DIAG(J)
    8 CONTINUE
   10 CONTINUE
      X(1)=X(1)/DIAG(1)
      DO 14 I=2,NPAR
      S=X(I)
      JMX=I-1
      DO 13 J=1,JMX
   13 S=S-MAT(J,I)*X(J)
   14 X(I)=S/DIAG(I)
      X(NPAR)=X(NPAR)/DIAG(NPAR)
      DO 16 K=2,NPAR
      I=NPAR+1-K
      JMI=I+1
      S=X(I)
      DO 15 J=JMI,NPAR
   15 S=S-MAT(I,J)*X(J)
   16 X(I)=S/DIAG(I)

      RETURN
      END

C -----------------------------------------------------------
      SUBROUTINE FAULT(XX,YY,XP,YP,AL,W,D,AX,AY,AZ,IND)
C -----------------------------------------------------------

C ----------Calcul du deplacement dans le repere geographique
C           J.C.Ruegg-1986 d'apres Okada(1985,BSSA)
C           Modifie le 30/8/87

C           Faille en tension  I= 1
C           Faille en coulissage  I= 2
C           Faille plongeante  I= 3
C           Coordonnes initiales XX,YY dans repere geographique UTM
C           Repere local de faille X1,Y1
C           Repere Okada  X2,Y2
C           Sortie : Deplacements A,B,C surface libre

      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COMOB/Q,ST,CT,CTT,SF,CF,ELAS,US(3)
      PI= 3.14159265359

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*CTT

C ----------Coins de la faille dans le repere Okada
      QS1= X2+AL
      QS2= X2-AL
      P1= Y2*CT +D*ST
      P2= P1 - W
      Q= Y2*ST - D*CT

C ----------Boucle de calcul aux quatre coins de la faille
      CALL DISLOK (QS1,P1,UX1,UY1,UZ1,IND)
      CALL DISLOK (QS1,P2,UX2,UY2,UZ2,IND)
      CALL DISLOK (QS2,P1,UX3,UY3,UZ3,IND)
      CALL DISLOK (QS2,P2,UX4,UY4,UZ4,IND)

      A2 = UX1-UX2-UX3+UX4
      B2 = UY1-UY2-UY3+UY4
      C2 = UZ1-UZ2-UZ3+UZ4

C ----------Deformation dans le repere geographique
      AX=(A2*SF - B2*CF)/2./PI
      AY=(A2*CF + B2*SF)/2./PI
      AZ=C2/2./PI

      RETURN
      END

C -------------------------
C     SOUS-PROGRAMME DISLOK
C -------------------------

C ----------Calcul des deplacements ux,uy,uz
C           associes a une dislocation
C           FORMULES DE OKADA (1985, BSSA)
C           D'apres J.C. RUEGG, 1986

      SUBROUTINE DISLOK(QSI,ETA,FUX,FUY,FUZ,IND)

      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COMOB/Q,ST,CT,CTT,SF,CF,ELAS,US(3)

      BY=ETA*CT+Q*ST
      BD=ETA*ST-Q*CT
      R=DSQRT(QSI**2+ETA**2+Q**2)
      XX=DSQRT(QSI**2+Q**2)

      RE=R+ETA
      RQ=R+QSI
      QRE=Q/R/RE
      QRQ=Q/R/RQ
      ATQE=DATAN(QSI*ETA/Q/R)

      AT= ETA*(XX+Q*CT)+XX*(R+XX)*ST
      BT= QSI*(R+XX)*CT

      FUX=0.
      FUY=0.
      FUZ=0.

      AI4= ELAS*(DLOG(R+BD)-ST*DLOG(R+ETA))/CT
      AI5= 2*ELAS*DATAN(AT/BT)/CT
      AI3=ELAS*(BY/CT/(R+BD)-DLOG(R+ETA))+AI4/CTT
      AI1= -ELAS*QSI/CT/(R+BD)-AI5/CTT
      AI2= -ELAS*DLOG(R+ETA)-AI3


C ----------Faille en tension
      IF(US(1).NE.0) THEN
      FUX=(Q*QRE - AI3*ST*ST)*US(1)
      FUY=-(BD*QRQ+ST*QRE*QSI-ST*ATQE+AI1*ST*ST)*US(1)
      FUZ=(BY*QRQ+CT*QSI*QRE-CT*ATQE-AI5*ST*ST)*US(1)
      ENDIF

C ----------Faille de coulissage
      IF(US(2).NE.0) THEN
      FUZ=FUZ-(BD*QRE+Q*ST/RE+AI4*ST)*US(2)
      FUX=FUX-(QSI*QRE+ATQE+AI1*ST)*US(2)
      FUY=FUY-(BY*QRE+Q*CT/RE+AI2*ST)*US(2)
      ENDIF

C ----------Faille plongeante
      IF(US(3).NE.0) THEN
      FUZ=FUZ+(-BD*QRQ-ST*ATQE+AI5*ST*CT)*US(3)
      FUX=FUX+(-Q/R+AI3*ST*CT)*US(3)
      FUY=FUY-(BY*QRQ+CT*ATQE-AI1*ST*CT)*US(3)
      ENDIF

      RETURN
      END

C -------------------------------------------
C     SOUS-PROGRAMME DE CALCUL DE DISLOCATION
C -------------------------------------------

      SUBROUTINE DISLO(NPOINTS,NFAIL1,NFAIL2,nff)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NPM=25000,NFM=50,nfi=9)
        COMMON /DIS/ POINT(NPM,2),DEPL(NPM,2),FAILLE(NFM,10),NV(nfi,2),
     &  IRV(nfi),tsa(nfi,3),tsb(nfi,3),tsc(nfi,3),dref(nfi,2)
      COMMON /COMOB/ Q,ST,CT,CTT,SF,CF,ELAS,US(3)
        COMMON/TYP/ityp,igeo,facs,facd
      PI=3.14159265339
c
C ----------Le coefficient d'elasticit‚ est choisi egal a 0.5
      ELAS=0.5
c     NVH=NV(1)+NV(2)
c
      do kf=nfail1,nfail2
      if(ityp.eq.2)then
         call sr2usd(faille(kf,8),faille(kf,9))
         call new2jcr(faille(kf,1),faille(kf,2),faille(kf,4),
     & faille(kf,6),faille(kf,3),faille(kf,7))
      else
         call sr2usd(faille(kf,8),faille(kf,9))
      endif
c     write(*,*)(faille(kf,jj),jj=1,10)
      enddo
c
      do if=1,nff
       dref(if,2)=0.d0 
      enddo
     
c
C ----------Vidage de DEPL
      DO 1, I=1,NPOINTS
      DEPL(I,2)=0
    1 CONTINUE

C ----------Boucle de calcul par faille
      DO 2, I=NFAIL1,NFAIL2
c
C ----------Transformation des parametres contenus dans FAILLE aux
C           parametres utilises par le sous-programme FAULT
      FA=FAILLE(I,3)*PI/180
      US(1)=FAILLE(I,10)
      US(2)=FAILLE(I,8)
      US(3)=-FAILLE(I,9)
      TA=FAILLE(I,7)*PI/180
      CT=DCOS(TA)
      ST=DSIN(TA)
      CTT=CT/ST
      CF=DCOS(FA)
      SF=DSIN(FA)
      W=FAILLE(I,6)
      RA=FAILLE(I,5)
      H=FAILLE(I,4)
      D=H+W*ST
      XPA=FAILLE(I,1)-H*CF*CTT
      YPA=FAILLE(I,2)+H*SF*CTT
c     print*,'PARAMETRES'
c     print*,faille(i,1),faille(i,2),fa,ta,h,w,ra,us(1),us(2),us(3)
c
c----------Boucle de calcul par fichier
c
      do 5 if=1,nff
       do 4 j=nv(if,1),nv(if,2)
        XX=POINT(J,1)
        YY=POINT(J,2)
        CALL FAULT(XX,YY,XPA,YPA,RA,W,D,AX,AY,AZ,1)
        tsx=tsa(if,1)*xx+tsb(if,1)*yy+tsc(if,1)
        tsy=tsa(if,2)*xx+tsb(if,2)*yy+tsc(if,2)
        tsz=tsa(if,3)*xx+tsb(if,3)*yy+tsc(if,3)
        DEPL(J,2)=DEPL(J,2)+tsz*AZ+tsx*AX+tsy*AY
c       print*,DEPL(J,2)
 4     continue
c
C ----------Reference eventuelle aux premieres donnees des listes de
C           mesures verticales et horizontales
        IF(IRV(if).NE.0) THEN
          irel=nv(if,1)+irv(if)-1
          rell=DEPL(irel,2)
          dref(if,2)=dref(if,2)+rell
          DO 3, kf=nv(if,1),nv(if,2)
            DEPL(kf,2)=DEPL(kf,2)-rell
 3        enddo
        ENDIF
c
 5      continue
c
      if(ityp.eq.2)then
         call usd2sr(faille(i,8),faille(i,9))
         call jcr2new(faille(i,1),faille(i,2),faille(i,4),
     & faille(i,6),faille(i,3),faille(i,7))
      else
         call usd2sr(faille(i,8),faille(i,9))
      endif
 2    continue
c
      END

C ----------------------------
C     SOUS PROGRAMME DE SORTIE
C ----------------------------

      SUBROUTINE SORT(NINV,NPOINTS,NF,nff)

      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NPM=25000,NINVM=15,NFM=50,NFMX9=400,nfi=9)
      character*80ficout1,ficout2,ficout3,nom1,nom2,nom3,form
        COMMON /DIS/ POINT(NPM,2),DEPL(NPM,2),FAILLE(NFM,10),NV(nfi,2),
     &  IRV(nfi),tsa(nfi,3),tsb(nfi,3),tsc(nfi,3),dref(nfi,2)
        COMMON /SOR/ DEPLS(NINVM,NPM),FAILS(NINVM,NFMX9)
        COMMON/TYP/ityp,igeo,facs,facd
        COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
C ----------Tableaux
      OPEN(7,FILE='FAIL')

C ----------Remplissage de FAIL
c
      ninv=ninv+1
      do j=1,ninv
       do i=1,nf*10,10
       if(ityp.eq.2)then
          if(igeo.eq.2) then
            CALL UTMGEO(fails(j,i),fails(j,i+1)
     #      ,fails(j,i),fails(j,i+1))
          else
          fails(j,i)=fails(j,i)/1000.
          fails(j,i+1)=fails(j,i+1)/1000.
          endif
          fails(j,i+3)=fails(j,i+3)/1000.
          fails(j,i+4)=fails(j,i+4)/1000.
          fails(j,i+5)=fails(j,i+5)/1000.
       else
          if(igeo.eq.2) then
            CALL UTMGEO(fails(j,i),fails(j,i+1)
     #      ,fails(j,i),fails(j,i+1))
          endif
       endif
       enddo
      enddo
      kl=0 
      if(ityp.eq.1) then
         DO 2, I=1,NF*10,10
         kl=kl+1
         write(7,*)kl
         WRITE(7,110) (FAILS(J,I),J=1,NINV)
         write(7,111)(fails(j,i+1),j=1,ninv)
         write(7,112)(fails(j,i+2),j=1,ninv)
         write(7,114)(fails(j,i+3),j=1,ninv)
         write(7,116)(fails(j,i+4),j=1,ninv)
         write(7,120)(fails(j,i+5),j=1,ninv)
         write(7,113)(fails(j,i+6),j=1,ninv)
         write(7,117)(fails(j,i+7),j=1,ninv)
         write(7,118)(fails(j,i+8),j=1,ninv)
         write(7,119)(fails(j,i+9),j=1,ninv)
 2       continue
      else
         do i=1,nf*10,10
         kl=kl+1 
         write(7,*)kl
         write(7,110)(fails(j,i),j=1,ninv)
         write(7,111)(fails(j,i+1),j=1,ninv) 
         write(7,112)(fails(j,i+2),j=1,ninv) 
         write(7,113)(fails(j,i+6),j=1,ninv) 
         write(7,114)(fails(j,i+3),j=1,ninv) 
         write(7,115)(fails(j,i+5),j=1,ninv) 
         write(7,116)(fails(j,i+4),j=1,ninv) 
         write(7,117)(fails(j,i+7),j=1,ninv) 
         write(7,118)(fails(j,i+8),j=1,ninv) 
         write(7,119)(fails(j,i+9),j=1,ninv) 
         enddo
      endif
      close(7)
 110  format(' x position',25f16.6)
 111  format(' y position',25f16.6)
 112  format(' strike    ',25f16.6) 
 113  format(' dip       ',25f16.6) 
 114  format(' h1        ',25f16.6) 
 115  format(' h2        ',25f16.6) 
 116  format(' demi-long ',25f16.6) 
 117  format(' slip      ',25f16.6) 
 118  format(' rake      ',25f16.6) 
 119  format(' ouverture ',25f16.6) 
 120  format(' w largeur ',25f16.6)
c
c
      open (56,file='out.fault')
       ij=ninv
       write(56,*)nf,ninv-1,' 0 0'
c
       write(99,*)'*********** RESU *************'
       if(ityp.eq.2) then
       do 22, i=1,nf*10,10
        write(56,*)(i-1)/10+1
        write(56,*)fails(ij,i),wq,fails(ij,i+1),wq
        write(56,*)fails(ij,i+2),wq,fails(ij,i+6),wq
        write(56,*)fails(ij,i+3),wq,fails(ij,i+5),wq,fails(ij,i+4),wq
        write(56,*)fails(ij,i+7),wq,fails(ij,i+8),wq,fails(ij,i+9),wq
        write(99,*)(i-1)/10+1
        write(99,*)fails(ij,i),wq,fails(ij,i+1),wq
        write(99,*)fails(ij,i+2),wq,fails(ij,i+6),wq
        write(99,*)fails(ij,i+3),wq,fails(ij,i+5),wq,fails(ij,i+4),wq
        write(99,*)fails(ij,i+7),wq,fails(ij,i+8),wq,fails(ij,i+9),wq
 22    continue
       else
       do 23 i=1,nf*10,10
        write(56,*)fails(ij,i),wq,fails(ij,i+1),wq
        write(56,*)fails(ij,i+2),wq
        write(56,*)fails(ij,i+3),wq
        write(56,*)fails(ij,i+4),wq,fails(ij,i+5),wq
        write(56,*)fails(ij,i+6),wq
        write(56,*)fails(ij,i+7),wq,fails(ij,i+8),wq,fails(ij,i+9),wq
        write(99,*)fails(ij,i),wq,fails(ij,i+1),wq
        write(99,*)fails(ij,i+2),wq
        write(99,*)fails(ij,i+3),wq
        write(99,*)fails(ij,i+4),wq,fails(ij,i+5),wq
        write(99,*)fails(ij,i+6),wq
        write(99,*)fails(ij,i+7),wq,fails(ij,i+8),wq,fails(ij,i+9),wq
 23    enddo
       endif
       close(56)
c
C ----------Remplissage de  *.out
c
         rmst=0d0
        do jj=1,nff
         write(99,*)JJ
         nom1='data'
         nom2='modl'
         nom3='diff'
         form='(a4,i1,a4)'
         write(ficout1,form)nom1,jj,'.out'
         write(ficout2,form)nom2,jj,'.out'
         write(ficout3,form)nom3,jj,'.out'
         open(89,file=ficout3)
         open(88,file=ficout1)
         open(87,file=ficout2)
         rms=0 
         nndd=0
         print*,'references ',dref(jj,1),dref(jj,2)
         write(99,*)' reference donnees modl',dref(jj,1),dref(jj,2)
         do kk=nv(jj,1),nv(jj,2)
          xg=point(kk,1)
          yg=point(kk,2)
          nndd=nndd+1
          deff=(depl(kk,2)+dref(jj,2))
          ddaa=(depl(kk,1)+dref(jj,2))
          adiff=ddaa-deff
          if (igeo.eq.2) then
            call utmgeo(xg,yg,xg,yg)
          else
            xg=xg/facs
            yg=yg/facs
          endif
          rms=adiff**2+rms
          rmst=adiff**2+rmst
          write(87,*)xg,yg,sngl(deff)
          write(88,*)xg,yg,sngl(ddaa)
          write(89,*)xg,yg,sngl(adiff)
         enddo
         rms=dsqrt(rms/dble(nndd))
         write(89,*)' RMS= ',rms
         write(99,*)'ndata RMS= ',nndd,rms
         print*,' nb data, RMS = ',nndd,rms
         close(89)
         close(88)
         close(87)
        enddo
        rmst=dsqrt(rmst/dble(npoints))
         print*,npoints
        print*,' RMS =',rmst
c
        write(99,*)''
        write(99,*)' TOUT ndata rms ',npoints,rmst
c
c---- remplissage de DIFF

      OPEN(15,FILE='DEPL')
      DO 1, I=1,NPOINTS
      WRITE(15,11) ((DEPLS(K,I)/facd),K=1,NINV)
 1    continue
      CLOSE(15)
   11 FORMAT(100F16.6)
   12 FORMAT(I6,100F16.6)
      END
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 METRES
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 ----------------------------------
      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 RETURN
      END
C ------------------------------------------------
      SUBROUTINE UTMGEO(XX,YY,xLON,yLAT)
C --------------------------------------
C         TRANSFORMATION  UTM - GEOGRAPHIQUE
C           ENTREE :  COORD. UTM  EN  METRES
C           SORTIE :  COORD. GEOGRAPH. EN degres
C
C           HEMISPHERE NORD : KH=1
C                      SUD  : KH=-1
C
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
      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
      ylat=ylat*0.9
      xlon=xlon*0.9
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 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----------------------------------------
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 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
c
        slip=dsqrt(us**2+ud**2)
        rake=-datan2d(ud,us)
c
        return
        end
c
c-----------------------------------------
c
c--------------------------------------------
          subroutine planv(igeo,fivec,kk)
c--------------------------------------------
c
        IMPLICIT REAL*8(A-H,O-Z)
        PARAMETER (NPM=25000,NFM=50,NFMX9=400,NINVM=15,nfi=9)
       character*80 fivec,nul*1
       dimension px(10),py(10),tx(10),ty(10),tz(10)
c
c       common/data2/nv,tsa,tsb,tsc,irv,dref
        COMMON /DIS/ POINT(NPM,2),DEPL(NPM,2),FAILLE(NFM,10),NV(nfi,2),
     &  IRV(nfi),tsa(nfi,3),tsb(nfi,3),tsc(nfi,3),dref(nfi,2)
         COMMON/COMOG/CC,E2,PI,IFU,AK,KH
c
        amil=1000.d0
c
        open(35,file=fivec)
c         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)
          print*,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))
c          px(i)=px(i)*amil
c          py(i)=py(i)*amil
          endif
        enddo
        close(35)
        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
        do k=1,3
          print*,sngl(tsa(kk,k)),sngl(tsb(kk,k)),sngl(tsc(kk,k))
        enddo
c
       do i=1,nd
          ttx=tsa(kk,1)*px(i)+tsb(kk,2)*py(i)+tsc(kk,1)
          tty=tsa(kk,2)*px(i)+tsb(kk,2)*py(i)+tsc(kk,2) 
          ttz=tsa(kk,3)*px(i)+tsb(kk,3)*py(i)+tsc(kk,3) 
          print*,px(i),py(i),ttx,tty,ttz
       enddo
       return
       end