c------------------------------------------------
c modification jbc pour gmt et traitement image 8000xn
c   creation de fichiers temporaires
c  modif avec defgentvec --> possibilite de demander plusieurs projts
c fabrique les fichiers temporaires dans le directory courant
c
C     Programme   D E F O R M XYZ  - option DEFORM-SS
c     adaptÇ de DEFORMXY pour sortie ligne sur fichier grapher
c     sans appel a biblio graphique calhpgl
c     
c     SORTIE : grapher : x,y, dZ
c                      : x,y, mod.(dX,dY,dZ)
c                      : x,y, proj. sur vecteur Satellite- Terre
c                      : direct image SAR synthetique...
C------------------------------------------------
c     ATTENTION : ce progr. utilise le s/p SUBDEFO -> memes dimensions de common
C       CALCUL DES DEFORMATIONS ET DESSINS dans le plan et sur coupe
c       -- calcul direct des deplacements a la surface d'un demi espace
c       elastique associÇs Ö une sÇrie de failles rectangulaires
c       modelisees par des dislocations (F. en coulissage, dip-slip, tension)
c       Formules analytiques d'Okada (1985, BSSA) adaptees par :
c -----------------
c        ( d'apres DEFORMP avec P pour imPrimante - XY pour Horizontal)
c          Version revisee de avril et sept 94 d'apres version decembre 90
c          J.C. Ruegg, Institut de Physique du Globe, Paris - 1995
c          -------  originale :  04 nov. 1988 - rev. jan 89 et sept 89 et 01/95.
c       Adapte de DESDEFOR et YDEFORM, utilisable comme tel --
c          Modifications par rapport aux versions anterieures (Desdefor,Coupe)
c          sur la definition de la faille et les entrees sorties 
c---------------------------------------------------------
c       Options :
c          - Calcul des deformations a partir d'un fichier de points connus
c          - Definition d'une grille
c          - lecture d'un fichier topo
c       Sorties :
c          - numeriques ecran et fichier edition
c---------------------------------------------------------
C ----- Elements de definition des failles :
C          - XS, YS : coordonnees du point milieu superieur de la faille
C          - PHI : Azimut de la faille, la direction
C            de plongement etant a droite de la direction definie
C            par PHI  - choix des unites (VANG): grades (1) ou degres (0.9) 
C          - US  : amplitude de la dislocation strike-slip ( m )
C                  US > 0 - faille SENESTRE
C                  US < 0 - faille DEXTRE
C          - UD  : amplitude de la dislocation dip-slip ( m )
C                  UD > ,0 - faille NORMALE
C                  UD < 0 - faille INVERSE
c          - UT : amplitude de l'ouverture d'une fissure en tension
c          - H   : profondeur du pt milieu de l'arete superieure
c          - D ou W : largeur de la faille (le long ligne de + grande pente
C          - D1, D2 : distances entre limites superieure et inferieure
C            de la faille et la trace en surface de la faille (km)
C          - R : demi-longueur de la faille (km)
C          - TETA : angle de plongement ( 0 - 100 gr.(90¯))
C-------------------------------------------------------------------------
      parameter (nbpc=3000,nbf=200,nbc=256)
      IMPLICIT REAL*8(A-H,L,O-Z)
      CHARACTER*80 FICH1,FICH2,FICH4,nom*80
      character*80 ficht1,ficht2,ficht3,ficvec
      character*1 nomout*80,nul
      CHARACTER*80 TEXT1
      character*20 xg,yg,vg,form,v20,w20,wg
      character*300 script
      real*4 DX(nbpc),DY(nbpc),DZ(nbpc)
      dimension v(nbpc),u(nbpc)
      dimension ta(4),tb(4),tc(4)
      dimension idep(10),isor(5)
      logical geo
      COMMON/F1/XF(nbf),YF(nbf),FF(nbf),H(nbf),W(nbf),R(nbf),TF(nbf)
      COMMON/F2/US(nbf),UD(nbf),UT(nbf)
      COMMON/COMOB/Q,ST,CT,ELAS,U0
      COMMON/COMOU/CC,E2,PI,IFU,AK,KH

      common/plan/px(4),py(4),tsx(4),tsy(4),tsz(4),npp
      common/exit/npx,npy,vmin,vmax,xc1,yc1,pas
C---------------------------------------------
c   VARIABLES A DEFINIR
c
c  defac =  facteur multiplicatif de la deformation
c      defac=3.d0/5.d0
       defac=1.d0
c     hemisphere kh=1 (nord), kh=-1 (sud)
      kh=-1
c
c     Vecteur moyen Terre-Satellite pour les orbites montantes de l'Etna
c     ts_x = 0.3791
c     ts_y = -0.09103
c     ts_z = 0.9208
c  kozani
c     ts_x=-0.402
c     ts_y=0.083
c     ts_z=-0.912
c chili
      ts_x=-0.36
      ts_y=0.09
      ts_z=-0.93
c
c landers 
      ts_x=-0.401
      ts_y=0.0636
      ts_z=-0.9135
c
c      frint=(0.0284/3)*5
       frint=(0.0284)
c      frint=(0.0284/3)
c
c     fichiers temporaires
c
      ficht1='./defx'
      ficht2='./defy'
      ficht3='./defz'
c
c-----------------------------------------------
c      valeurs de depart et increment des ecritures
c   dans param.txt
c
       xxdeb=-71.1
       yydeb=-22.1
       yyinc=-.1
c
c-----------------------------------------------
c-----------------------------------------------
C     CONSTANTES
c     transformation de m en km 
       FAC=1000.d0
       ISOW=1
c..........ISOW=0 sortie parametres de faille classique
c              =1 sortie sous forme coord. des pts milieux et extremites
      PI= 3.14159265359
      CC= 6399593.6257585d0
      PID=PI/2.d0
      AMIL=1000.
      AMOT=0.0001d0
      ELAS= 0.5
      I=0
      A=0.
      B=0.
      C=0.
c
c     ilo (unite logique fichier de sortie
      ilo=44
c
C -------- Parametres de failles:
c
      print*,'type de fichier parametre :'
      print*,' 1 ------> type jcr'
      print*,' 2 ------> type jbc'
      read(*,*) ityp
 2    print*,'entrees en : 1 --> UTM'
      print*,'             2 --> Geographique'
      read(*,*)isys
      if(isys.eq.1) then
         geo=.false.
      elseif(isys.eq.2) then
         geo=.true.
         print*,'fichier de projection'
         read(*,'(a)')fich4
         open(10,file=fich4)
         read(10,'(a)')script
         read(10,*)cc,eprim
         e2=eprim*eprim
         read(10,*)ifu
         close(10)
      else
         goto 2
      ENDIF
c
   10 write(*,*) 'nom du cas (sans extension, <9 caractere)'
      WRITE(*,*) ' le fichier parametre de faille doit s"appeler :'
      write(*,*)'nom_du_cas.fault'
      read(*,'(a)')nom
      call tail (nom,itai)
      if(itai.gt.9) goto 10
      print*,itai,nom
      write (form,'(a2,i1,a4)')'(a',itai,',a6)'
      print*,form
      write(fich1,form)nom,'.fault'
       print*,fich1
c      READ(*,'(A)',err=10) FICH1
      OPEN(31,STATUS='OLD',FILE=FICH1,err=10)
      print*,'ouverture ',fich1
C
c  11 WRITE(*,*) ' Nom du fichier edition - .ed :'
      write (form,'(a2,i1,a4)')'(a',itai,',a3)'
      write(fich2,form)nom,'.ed'
      OPEN(42,FILE=FICH2)
      print*,'ouverture ',fich2
c
C -------- Titre :
      if(ityp.eq.1) then
      READ(31,'(A)') TEXT1
      WRITE(*,'(A)') TEXT1
      WRITE(42,'(A)') TEXT1
      WRITE(42,98) FICH1
      endif
   98 format(' Nom du fichier parametres failles: ',A12)
      write(42,*)'  --------------------------------------------'
c
c
      open (66,file='trfault.xy')
      open (67,file='hafault.xy')
      open (77,file='slipd.txt')
      open (88,file='slips.txt')
      open(92,file='slip.txt')
      open (94,file='slip.xy')
      open(91,file='param.txt')
      write(91,*)'param : x,y,z,strike,dip,h1,h2,slip,rake,long'
      open(93,file='moment.txt')
      write(66,'(a1)')'>'
      write(67,'(a1)')'>'
      write(94,*)'x y z slip slips slipd '
C
C -------- Entree des 5ailles :
C              NFA : nb de failles
c              IFA : no de la faille; fin par ifa=0 
      write(*,*) ' Unite des valeurs angulaires : grad =1., deg=0.9 :'
      if (ityp.eq.1) then
        read(31,*) VANG
      else
        read (31,*)nff
        vang=.9
      endif
c
  101 I=I+1
      if(ityp.eq.2.and.i.gt.nff)goto 99
      READ(31,*) IFA
      IF(IFA.LT.1) GOTO 99
C
C------------XS,YS Coord. du milieu du bord supÇrieur de la faille
c            XQ,YQ coord. du milieu de la trace en surface
c            PHI, Azimut de la trace de la faille, plong. a droite
c            HA, profondeur de ce point - WA, largeur de la faille
c            RA, demi-longueur , TETA plongement
c            
      if(ityp.eq.1) then
        READ (31,*) XS,YS,PHI
        READ (31,*) USA,UDA,UTA
        READ (31,*) HA,WA,RA,TETA
        xxx=xs
        yyy=ys
        hb=wa
        if(geo) then
          if (ys.lt.0.d0) then
             kh=-1
          else
             kh=1
          endif
          call geoutm(XS,YS,XS,YS)
          xs=xs/amil
          ys=ys/amil
        endif        
c        call sr2usd(USA,UDA)
c        usa=usa/100
c        uda=uda/100
c        uta=uta/100
      else
        read (31,*) xs,wz,ys,wz
        read (31,*) phi,wz,teta,wz
        read (31,*) ha,wz,wa,wz,ra,wz
        read (31,*) usa,wz,uda,wz,uta,wz
        xxx=xs
        yyy=ys   
        hb=wa
        if(geo) then
          if (ys.lt.0.d0) then
             kh=-1
          else
             kh=1
          endif
          call geoutm(XS,YS,XS,YS) 
          xs=xs/amil 
          ys=ys/amil 
        endif
        call sr2usd(USA,UDA)
        usa=usa/100
        uda=uda/100
        uta=uta/100
        call new2jcr(xs,ys,ha,wa,phi,teta)
      endif
      if(usa.eq.0.d0.and.uda.eq.0.d0) then
        usa=0.00000001d0
        uda=0.00000001d0
      endif
      print*,XS,YS,PHI,USA,UDA,UTA,HA,WA,RA,TETA
c
      FANG=PI/200./VANG
      TFI=TETA*FANG
      TF(I)=TFI
      FFI=PHI*FANG
      FF(I)=FFI
      DA1=HA/DSIN(TF(I))
      DA2=DA1+WA
c
c..... Q = milieu trace en surface,  M = centre de la faille
c  XQ et YQ  FAUX!!!!!
      XQ=XS-DA1*dcos(FFI)
      YQ=YS+DA1*dsin(FFI)
      WM=WA/2.*dcos(TFI)
      XM=XS+WM*dcos(FFI)
      YM=YS-WM*dsin(FFI)
      ZM=HA+WA/2.*dsin(TFI)
c..........S1, S2 extremites superieures
      ZS=HA
      XS1=XS-RA*dsin(FFI)
      YS1=YS-RA*dcos(FFI)
      XS2=XS+RA*dsin(FFI)
      YS2=YS+RA*dcos(FFI)
c ......... XI,YI,ZI coord. point inferieur centre
      DHSI=WA*dcos(TFI)
      DZSI=WA*dsin(TFI)
      XI= XS+DHSI*dsin(FFI+PID)
      YI= YS+DHSI*dcos(FFI+PID)
      ZI= ZS+DZSI
c..........I1, I2 extremites inferieures
      XI1=XI-RA*dsin(FFI)
      YI1=YI-RA*dcos(FFI)
      XI2=XI+RA*dsin(FFI)
      YI2=YI+RA*dcos(FFI)
c
c................xq1 yq1 xq2 yq2 trace faille en surface
      xq1=xs1
      yq1=ys1
      wq1=wa
      call jcr2new(xq1,yq1,ha,wq1,phi,teta)
      xq2=xs2
      yq2=ys2 
      wq2=wa 
      call jcr2new(xq2,yq2,ha,wq2,phi,teta)
c
c -------------------------
c
      XF(I)=XS*AMIL
      YF(I)=YS*AMIL
      US(I)=USA
      UD(I)=UDA
      UT(I)=UTA
c
      H(I)=HA*AMIL
      W(I)=WA*AMIL
      R(I)= RA* AMIL 
c
      PHID= FFI*180./PI
      TETD= TFI*180./PI
      TETG= TFI*200./PI
      PHIG= FFI*200./PI
c
c--- Surface : S = WA * RT km2, mu = 3.3 Exp 10 Pa, du= dÇplacement (m)
c    moment = AMO (N-m), Magnitude: AMW par Log Mo = 9 + 1.5 x Mw
c
      RT=RA+RA
      SK = WA*RT
      SM = SK * 1.D6
      DU2= USA*USA+UDA*UDA
      DU = DSQRT(DU2)
c
      XCIS= 3.3d10
      AMO = SM*XCIS*DU
      ALMO= DLOG10(AMO) 
      AMW = (ALMO-9.)/1.5
c
c--sortie des parametres de failles :
c
  111 FORMAT(' Faille no ',I3)
  112 FORMAT (10x,5(3X,F10.3)) 
c
      write(42,*) '                    Moment sismique - Magnitude :'
      write(42,*) '                    ---(avec mu=3.3)------------ '
c
  105 format(30x,'Moment :',1x,D12.4,2X, ' Mw =',F6.3)
c
c ------ sortie des coordonnees de tous les points de la faille
c
  102 WRITE(*,111) IFA
      WRITE(*,*)' Coordonnees des points extremite et milieu :'
      write(*,116)
      write(*,113) XS1,YS1, XS,YS, XS2,YS2, ZS
      write(*,114) XM,YM, ZM
      write(*,115) XI1,YI1, XI,YI, XI2,YI2, ZI
  113 format(' Bord superieur:',3(F8.3,1X,F8.3,2x),1x,F6.3) 
  114 format(' Milieu faille :',19x,F8.3,1X,F8.3,22x,F6.3) 
  115 format(' Bord inferieur:',3(F8.3,1X,F8.3,2x),1x,F6.3) 
  116 format(19x,'x1, y1 '12x,' xm, ym ',12x,' x2, y2 ',5x,'     z ')
      write(*,*) ' Glissement (m), Azimut, pendage, rake (deg.):'
      rake=- datan2(UDA,USA)*180./PI
      slip= dsqrt(USA**2+UDA**2)   
      write(*,112) slip, PHID, TETD,rake
      write(*,*) ' Dimensions faille (km), Surface (km2) :'
      write(*,112) WA, RT, SK
      write(*,*) ' Moment sismique - Magnitude :'
      write(*,105) AMO, AMW
c
c ---
c    
      print*,XS1,YS1, XS,YS, XS2,YS2, ZS
      print*,XM,YM, ZM
       print*,XI1,YI1, XI,YI, XI2,YI2, ZI
      WRITE(42,111) IFA
      WRITE(42,*)' Coordonnees des points extremite et milieu :'
      write(42,116)
      write(42,113) XS1,YS1, XS,YS, XS2,YS2, ZS
      write(42,114) XM,YM, ZM
      write(42,115) XI1,YI1, XI,YI, XI2,YI2, ZI
      write(42,*) ' -------------------------- '
      write(42,*) ' Glissement (m), Azimut, pendage, rake (deg.):'
      write(42,112) slip, PHID, TETD,rake
      write(42,*) ' Dimensions faille (km), Surface (km2) :'
      write(42,112) WA, RT, SK
      write(42,*) ' Moment sismique - Magnitude :'
      write(42,105) AMO, AMW
      AMOT=AMOT+AMO
c
 334  format(3(1x,a20))
      if(geo) then
        xs1=xs1*fac
        ys1=ys1*fac
        xs2=xs2*fac
        ys2=ys2*fac
        xi1=xi1*fac
        yi1=yi1*fac
        xi2=xi2*fac
        yi2=yi2*fac
        xm=xm*fac
        ym=ym*fac
        call utmgeo(xs1,ys1,xs1,ys1)
        call utmgeo(xs2,ys2,xs2,ys2)
        call utmgeo(xi1,yi1,xi1,yi1)
        call utmgeo(xi2,yi2,xi2,yi2)
        call utmgeo(xq2*fac,yq2*fac,xq2,yq2)
        call utmgeo(xq1*fac,yq1*fac,xq1,yq1)
        call utmgeo(xm,ym,xm,ym)
      endif
      call cv20(xs1,xg)
      call cv20(ys1,yg)
      call cv20(-zs,vg)
      write(67,*)xq1,yq1,zq 
      write(67,*)xq2,yq2,zq  
      write(66,*)xq1,yq1,zq
      write(66,*)xq2,yq2,zq
      write(66,'(a1)')'>'
c     write(66,334)xg,yg,vg
      write(66,*)real(xs1),real(ys1),real(-zs)
      call cv20(xs2,xg)
      call cv20(ys2,yg)
c     write(66,334)xg,yg,vg
      write(66,*)real(xs2),real(ys2),real(-zs) 
      call cv20(xi2,xg)
      call cv20(yi2,yg)
      call cv20(-zi,vg)
c     write(66,334)xg,yg,vg
      write(66,*)real(xi2),real(yi2),real(-zi) 
      call cv20(xi1,xg)
      call cv20(yi1,yg)
c     write(66,334)xg,yg,vg
      write(66,*)real(xi1),real(yi1),real(-zi) 
      call cv20(xs1,xg)
      call cv20(ys1,yg)
      call cv20(-zs,vg)
c     write(66,334)xg,yg,vg
      write(66,*)real(xs1),real(ys1),real(-zs) 
      write(66,'(a1)')'>'
      write(67,'(a1)')'>'
c
      if(geo) then
        write(42,*)'coordonnees de la faille en degres'
        write(42,400)xs1,ys1,xs2,ys2,zs
        write(42,401)xm,ym,zm
        write(42,400)xi1,yi1,xi2,yi2,zi
      endif
 400  format(2(f11.5),2x,2(f11.5),2x,f6.3)
 401  format(10x,2(f11.5),8x,f6.3)
      call cv20(xm,v20)
      call cv20(ym,w20)
      write(77,335)v20,w20,uda
      write(88,335)v20,w20,usa
      write(92,335)v20,w20,slip
      write(94,336)v20,w20,zm,slip,usa,uda
 336  format(a20,1x,a20,1x,f9.3,3(1x,f6.2))
 335  format(a20,1x,a20,' 12 0. 1 4 ',f6.2)
      yzz=yydeb+(i-1)*yyinc*3
      xzz=xxdeb
c
      write(91,391)float(xm),float(ym),float(-zm),float(phid),float
     &(TETD),float(-zi),float(-zs),float(slip),float(rake),float(rt)
 391   format(2(E15.8,2X),f8.3,2x,f5.0,2x,f4.0,2x,
     & 2(f7.2,2x),f6.2,2x,f5.0,1x,f7.2)
      GOTO 101
C
  99  NFA=I-1
      CLOSE(31)
      close(66)
      close(67)
      close (77)
      close (88)
      close(92)
      close(91)
      close(94)
c
      ALMO= dlog10(AMOT) 
      AMW = (ALMO-9.)/1.5
c
      WRITE(*,*) ' Moment sismique Total - Magnitude Mw associee :'
      write(*,*) ' ---------------------------------------------- '
      WRITE(42,*) ' Moment sismique Total - Magnitude Mw associee :'
      write(42,*) ' ---------------------------------------------- '
      write(42,105) AMOT, AMW
      write(*,105) AMOT, AMW
      write(93,393)AMOT, AMW
 393  format(' 1.0 1.0 15 0.0 1 4 Moment='D12.4,2X,' Mw=',F6.3)
      write(42,*) ' ---------------------------------------------- '
c
      close(93)
c
c---------------------------------------------------------------------
C -------- Parametres des points de calcul :
c
 100  WRITE(*,*) '  DÇfinition d'' une GRILLE              --  : '
      IDE1=0
C
c -------------DÐfinition d'une grille :
C
      WRITE(*,*) '***- GRILLE NPX x NPY (=<1900); npx,npy :-***'
      READ(*,*) NPX,NPY
c
c      WRITE(*,*) ' Coord. du pt Haut Ö gauche XC1 , YC1 (km) :'
c      READ(*,*) XC1,YC1
c
      WRITE(*,*) ' Coordonnees du point en bas a gauche du dessin'
      if(.not.geo) then
      write(*,*) '        UTM     ---> XC1, YC1 (km)'
      else
      write(*,*) '        Degres dec-> LAT, LONG    '
      endif
      READ(*,*) XC1,YC1
      WRITE(*,*) xc1, yc1
      WRITE(42,*) xc1, yc1 
c
      if(.not.geo) then
        WRITE(*,*) 'PAS (en km):'
        READ(*,*) PAS
      else
        write(*,*) 'pas (en degres)'
        read(*,*)pas
      endif
      eps=0.000001d0
      K=0
c
      write(*,*) ' Calcul en cours.....'
c
  150 FORMAT(2(F10.1,2X),3(F8.4,3X))
  151 FORMAT(2(F10.1,2X),F8.4,3X,F8.2)
c
c-------------------------------
C     ---- CALCUL DEFORMATIONS :
c
C-------------------------------
c
      if(geo) then
        xdeb=xc1
        ydeb=yc1
        xfin=xdeb+(npx-1)*pas
        yfin=ydeb+(npy-1)*pas
        write (42,*) 'coordonnees des 4 coins de la grille'
        write(42,*) 'en coordonnees geographiques'
        write(42,403)xdeb,yfin,xfin,yfin
        write(42,403)xdeb,ydeb,xfin,ydeb
        call geoutm(xdeb,yfin,xdu,ydu)
        call geoutm(xfin,yfin,xfu,yfu)
        xdu=xdu/fac
        xfu=xfu/fac
        ydu=ydu/fac
        yfu=yfu/fac
        write (42,*) 'coordonnees des 4 coins de la grille'
        write(42,*) 'en coordonnees UTM'
        write(42,403)xdu,ydu,xfu,yfu
        call geoutm(xdeb,ydeb,xdu,ydu)
        call geoutm(xfin,ydeb,xfu,yfu)
        xdu=xdu/fac 
        xfu=xfu/fac 
        ydu=ydu/fac 
        yfu=yfu/fac 
        write(42,403)xdu,ydu,xfu,yfu
 403    format (2(f11.5),3x,2(f11.5))
      else
        write (42,*) 'coordonnees des 4 coins de la grille' 
        write(42,*) 'en coordonnees UTM' 
        write(42,403)xdeb,yfin,xfin,yfin 
        write(42,403)xdeb,ydeb,xfin,ydeb
      endif
c
c
      xmin=1.0e20
      xmax=-1.0e20
      ymin=1.0e20
      ymax=-1.0e20
      zmin=1.0e20
      zmax=-1.0e20
      hmin=1.0e20
      hmax=-1.0e20
      tmin=1.0e20
      tmax=-1.0e20
      smin=1.0e20
      smax=-1.0e20
      irec=npx*4
      open(97,file=ficht1,access='direct',recl=irec)
      open(98,file=ficht2,access='direct',recl=irec)
      open(99,file=ficht3,access='direct',recl=irec)
    1 DO 20 J=1,NPY
       Yw=YC1+PAS*(J-1)
       write(*,*) J
       DO 21 I=1,NPX
         Xw=XC1+PAS*(I-1)
         if(geo) then
c           print*,xw,yw,eee
            callgeoutm(xw,yw,xk,yk)
c           print*,xk,yk
            xk=xk/fac
            yk=yk/fac
         else
            xk=xw
            yk=yw
         endif
         xk=xk+eps
         yk=yk+eps
c         
         XX=XK*FAC
         YY=YK*FAC
c
         CALL DEFOR(NFA,XX,YY,DFX,DFY,DFZ)
c
c        print*,'&&&&&&&&&&&&&&&&&&&&&'
c        print*,xx,yy,dfx,dfy,dfz
c        print*,'&&&&&&&&&&&&&&&&&&&&&' 
         DX(I)=sngl(DFX*defac)
         DY(I)=sngl(DFY*defac)
         DZ(I)=sngl(DFZ*defac)
c
         xmin=dmin1(xmin,dble(dx(i)))
         xmax=dmax1(xmax,dble(dx(i)))
         ymin=dmin1(ymin,dble(dy(i)))
         ymax=dmax1(ymax,dble(dy(i)))
         zmin=dmin1(zmin,dble(dz(i)))
         zmax=dmax1(zmax,dble(dz(i)))
         ttt=dsqrt(dble(dx(i))**2+dble(dy(i))**2)
         hmin=dmin1(hmin,ttt)
         hmax=dmax1(hmax,ttt)
         ttt=dsqrt(dble(dx(i))**2+dble(dy(i))**2+dble(dz(i))**2)
         tmin=dmin1(tmin,ttt)
         tmax=dmax1(tmax,ttt) 
         ttt=dble(dx(i))*ts_x+dble(dy(i))*ts_y+dble(dz(i))*ts_z
         smin=dmin1(smin,ttt)
         smax=dmax1(smax,ttt) 
   21 continue
         write (97,rec=j)(dx(i),i=1,npx)
         write (98,rec=j)(dy(i),i=1,npx)
         write (99,rec=j)(dz(i),i=1,npx)
        write (*,*)(dx(i),i=1,npx)
        write (*,*)(dy(i),i=1,npx)
        write (*,*)(dz(i),i=1,npx)
   20 continue
c
c---------------------------------------------------------
c     
      nbts=0
      nbfr=0
      ik=1
 61   continue
      print*,'choix de la grandeur a representer'
      print*,'    1--------> deplacement vertical'
      print*,'    2--------> deplacement horizontal en vectoriel (GMT)'
      print*,'    3--------> deplacement proj. sur dir. TS'
      print*,'    4--------> frange interferometriques'
      print*,'    5--------> deplacement y'
      print*,'    6--------> deplacement x'
      print*,'    7--------> deplacement horizontal (amplitude)'
      print*,'    8--------> deplacement total'
c     print*,'    9-------> deplacement dans une direction particuliere'
      read(*,*)idep(ik)
      if(idep(ik).lt.1.or.idep(ik).gt.8) goto 61
c 
      if(idep(ik).ne.2)then
         print*,'format de sortie'
         print*,'1 -----> fichier code sur 1 octet'
         print*,'2 -----> fichier code sur 2 octets'
         print*,'3 -----> fichier real4 pour GMT'
         read(*,*)isor(ik)
      else
         print*,'echantillonage (1=1/1, 2=1/2, 3=1/3 ..)'
         read(*,*)iech
      endif
      print*,'autre type de grandeur a sortir 1/0'
      read(*,*) irep
      if(irep.eq.1) then
        ik=ik+1
        goto 61
      endif
      nk=ik
c
c     do 45 ik=1,nk
c     if(idep(ik).eq.4.or.idep(ik).eq.3) then
c       print*,'fichier de vecteur projection'
c       read(*,'(a)')ficvec
c       open(45,file=ficvec)
c        read(45,'(a1)')nul
c        read(45,*)npp
c        do iqp=1,npp
c          read(45,*)px(iqp),py(iqp),tsx(iqp),tsy(iqp),tsz(iqp)
c        enddo
c        call planvec(ta,tb,tc)
c        do iqp=1,3
c           print*,ta(iqp),tb(iqp),tc(iqp)
c        enddo
c        goto 46
c     endif
c45   continue
c46   continue
c
c
c
      do 300 ik=1,nk
        write(42,*)''
        write(42,*)'***************************************'
        write(42,*)'sortie numero ',idep(ik)
        print*,'sortie numero ',idep(ik)
        goto (311,312,313,314,315,316,317,318,319)idep(ik)
 311   continue
      write(42,*)'deplacement vertical'
      print*,'deplacement vertical'
      write (form,'(a2,i1,a4)')'(a',itai,',a6)'
      write(nomout,form) nom,'_depzz'
      itaio=itai+6
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif
      goto (3117,3118,3119)isor(ik)
 3117 write(fich4,form) nomout,'.1oc'
      nrec=npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3116
 3118 write(fich4,form) nomout,'.2oc'
      nrec=2*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3116
 3119 write(fich4,form) nomout,'.bin'
      nrec=4*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3116
 3116 continue
      vmax=zmax*100
      vmin=zmin*100
      kn=0
      do j=npy,1,-1
         kn=kn+1
         read  (99,rec=j)(dz(i),i=1,npx)
c        write (*,*)(dz(i),i=1,npx)
         do i=1,npx
           v(i)=dble(dz(i))*100.
         enddo
         goto (3111,3112,3113)isor(ik)
 3111    call wri1(v,ilo,kn)
         goto 3110
 3112    call wri2(v,ilo,kn)
         goto 3110
 3113    call wri4(v,ilo,kn)
         goto 3110
 3110    continue
      enddo
      close(ilo)
      write(42,form)fich4
      goto 310
c
 312   continue
      write(42,*)'deplacement horizontal en vectoriel'
      print*,'deplacement horizontal en vectoriel'
      write (form,'(a2,i1,a4)')'(a',itai,',a6)' 
      write(fich4,form) nom,'_vecxy' 
      vmax=hmax *100.
      vmin=hmin *100.
      open(ilo,file=fich4)
      do j=npy,1,-iech
         read  (97,rec=j)(dx(i),i=1,npx)
         read  (98,rec=j)(dy(i),i=1,npx)
         yk=yc1+pas*(j-1)
         do i=1,npx,iech
           xk=xc1+pas*(i-1)
           v(i)=dble(dx(i))*100.
           u(i)=dble(dy(i))*100.
           call cv20(xk,xg)
           call cv20(yk,yg)
           call cv20(v(i),vg)
	   call cv20(u(i),wg)
	   write(ilo,333) xg,yg,vg,wg,ww,ww,ww
         enddo
      enddo
 333  format(4(1x,a20),3(1x,f2.1))
      close (ilo)
      write(42,form)fich4
      goto 310
c
 313   continue
      write(42,*)'deplacement projete sur la direction terre-satellite'
      print*,'deplacement projete sur la direction terre-satellite'
       nbts=nbts+1
c
        print*,'fichier de vecteur projection'
        read(*,'(a)')ficvec
        open(45,file=ficvec)
         read(45,'(a1)')nul
         read(45,*)npp
         do iqp=1,npp
           read(45,*)px(iqp),py(iqp),tsx(iqp),tsy(iqp),tsz(iqp)
         enddo
         call planvec(ta,tb,tc)
         do iqp=1,3
            print*,ta(iqp),tb(iqp),tc(iqp)
         enddo
         close(45)
c
c     write (form,'(a2,i1,a4)')'(a',itai,',a6)' 
c     write(nomout,form) nom,'_depts' 
c     itaio=itai+6 
c     if(itaio.gt.9) then 
c       write (form,'(a2,i2,a4)')'(a',itaio,',a4)' 
c     else 
c       write (form,'(a2,i1,a4)')'(a',itaio,',a4)' 
c     endif 
c
      write (form,'(a2,i1,a7)')'(a',itai,',a6,i1)'
      write(nomout,form) nom,'_depts',nbts
      itaio=itai+7
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif

      goto (3137,3138,3139)isor(ik) 
 3137 write(fich4,form) nomout,'.1oc' 
      nrec=npx 
      open(ilo,file=fich4,access='direct',recl=nrec) 
      goto 3136  
 3138 write(fich4,form) nomout,'.2oc' 
      nrec=2*npx 
      open(ilo,file=fich4,access='direct',recl=nrec) 
      goto 3136  
 3139 write(fich4,form) nomout,'.bin' 
      nrec=4*npx 
      open(ilo,file=fich4,access='direct',recl=nrec) 
      goto 3136  
 3136 continue 
      vmax=-1.0e20
      vmin=1.0e20
      kn=0 
c     print*,npx,npy
c     print*,ta(1),tb(1),tc(1)
c     print*,ta(2),tb(2),tc(2) 
c     print*,ta(3),tb(3),tc(3) 
      do j=npy,1,-1
         yw=yc1+pas*(j-1)
         kn=kn+1
         read  (97,rec=j)(dx(i),i=1,npx)
         read  (98,rec=j)(dy(i),i=1,npx)
         read  (99,rec=j)(dz(i),i=1,npx)
        do i=1,npx
         xw=xc1+pas*(i-1)
           ts_x=ta(1)*xw+tb(1)*yw+tc(1)
           ts_y=ta(2)*xw+tb(2)*yw+tc(2)
           ts_z=ta(3)*xw+tb(3)*yw+tc(3)
c          print*,xw,yw
c          print*,dx(i),dy(i),dz(i)
c          print*,ts_x,ts_y,ts_z
           v(i)=dble(dx(i))*ts_x+dble(dy(i))*ts_y+dble(dz(i))*ts_z
c          print*,v(i)
           v(i)=v(i)*100.
           vmax=dmax1(v(i),vmax)
           vmin=dmin1(v(i),vmin)
        enddo
         goto (3131,3132,3133)isor(ik)
 3131    call wri1(v,ilo,kn)
         goto 3130
 3132    call wri2(v,ilo,kn)
         goto 3130
 3133    call wri4(v,ilo,kn) 
         goto 3130
 3130    continue
c
      enddo
      close (ilo)
      write(42,form)fich4
      print*,vmin,vmax
      goto 310
c
 314   continue
      write(42,*)'determination des franges interferometriques'
      print*,'determination des franges interferometriques'
      nbfr=nbfr+1
c
        print*,'fichier de vecteur projection'
        read(*,'(a)')ficvec
        open(45,file=ficvec)
         read(45,'(a1)')nul
         read(45,*)npp
         do iqp=1,npp
           read(45,*)px(iqp),py(iqp),tsx(iqp),tsy(iqp),tsz(iqp)
         enddo
         call planvec(ta,tb,tc)
         do iqp=1,3
            print*,ta(iqp),tb(iqp),tc(iqp)
         enddo
         close(45)
c
c     write (form,'(a2,i1,a4)')'(a',itai,',a6)' 
c     write(nomout,form) nom,'_inter' 
c     itaio=itai+6 
c     if(itaio.gt.9) then 
c       write (form,'(a2,i2,a4)')'(a',itaio,',a4)' 
c     else 
c       write (form,'(a2,i1,a4)')'(a',itaio,',a4)' 
c     endif 
c
      write (form,'(a2,i1,a7)')'(a',itai,',a6,i1)'
      write(nomout,form) nom,'_depfr',nbfr
      itaio=itai+7
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif
c
      goto (3147,3148,3149)isor(ik) 
 3147 write(fich4,form) nomout,'.1oc' 
      nrec=npx 
      open(ilo,file=fich4,access='direct',recl=nrec) 
      goto 3146  
 3148 write(fich4,form) nomout,'.2oc' 
      nrec=2*npx 
      open(ilo,file=fich4,access='direct',recl=nrec) 
      goto 3146  
 3149 write(fich4,form) nomout,'.bin' 
      nrec=4*npx 
      open(ilo,file=fich4,access='direct',recl=nrec) 
      goto 3146  
 3146 continue 
c
      vmin=0
      vmax=255
      kn=0 
      do j=npy,1,-1
         yw=yc1+pas*(j-1)  
         kn=kn+1
         read  (97,rec=j)(dx(i),i=1,npx)
         read  (98,rec=j)(dy(i),i=1,npx)
         read  (99,rec=j)(dz(i),i=1,npx)
        do i=1,npx
         xw=xc1+pas*(i-1) 
           ts_x=ta(1)*xw+tb(1)*yw+tc(1)  
           ts_y=ta(2)*xw+tb(2)*yw+tc(2)  
           ts_z=ta(3)*xw+tb(3)*yw+tc(3)
           amp=dble(dx(i))*ts_x+dble(dy(i))*ts_y+dble(dz(i))*ts_z
           v(i)=sngl((amp/frint-dint(amp/frint))*255)
           if(v(i).lt.0.)v(i)=v(i)+255
         enddo
  
         goto (3141,3142,3143)isor(ik)   
 3141    call wri1(v,ilo,kn)
         goto 3140
 3142    call wri2(v,ilo,kn)
         goto 3140
 3143    call wri4(v,ilo,kn) 
         goto 3140
 3140    continue
      enddo   
      close (ilo)
      write(42,form)fich4
      goto 310 
c
 315  continue
      write(42,*)'deplacement dans la direction y'
      print*,'deplacement dans la direction y'
      write (form,'(a2,i1,a4)')'(a',itai,',a6)' 
      write(nomout,form) nom,'_depyy' 
      itaio=itai+6
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif
      goto (3157,3158,3159)isor(ik)
 3157 write(fich4,form) nomout,'.1oc'
      nrec=npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3156
 3158 write(fich4,form) nomout,'.2oc'
      nrec=2*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3156
 3159 write(fich4,form) nomout,'.bin'
      nrec=4*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3156
 3156 continue
      vmin=ymin*100.
      vmax=ymax*100.
      kn=0 
      do j=npy,1,-1
         kn=kn+1
         read  (98,rec=j)(dy(i),i=1,npx)
        do i=1,npx
          v(i)=dble(dy(i))*100.
        enddo
         goto (3151,3152,3153)isor(ik)    
 3151    call wri1(v,ilo,kn)
         goto 3150
 3152    call wri2(v,ilo,kn)
         goto 3150
 3153    call wri4(v,ilo,kn) 
         goto 3150
 3150    continue
c
      enddo      
      close (ilo)
      write(42,form)fich4
      goto 310  
c 
 316  continue
      write(42,*)'deplacement dans la direction x'
      print*,'deplacement dans la direction x'
      write (form,'(a2,i1,a4)')'(a',itai,',a6)' 
      write(nomout,form) nom,'_depxx' 
      itaio=itai+6
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif
      goto (3167,3168,3169)isor(ik)
 3167 write(fich4,form) nomout,'.1oc'
      nrec=npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3166
 3168 write(fich4,form) nomout,'.2oc'
      nrec=2*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3166
 3169 write(fich4,form) nomout,'.bin'
      nrec=4*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3166
 3166 continue
      vmax=xmax *100.
      vmin=xmin *100.
      kn=0 
      do j=npy,1,-1
         kn=kn+1
         read  (97,rec=j)(dx(i),i=1,npx)
         do i=1,npx
          v(i)=dble(dx(i))*100.
        enddo
         goto (3161,3162,3163)isor(ik)    
 3161    call wri1(v,ilo,kn)
         goto 3160
 3162    call wri2(v,ilo,kn)
         goto 3160
 3163    call wri4(v,ilo,kn) 
         goto 3160
 3160    continue
c
      enddo      
      close (ilo)
      write(42,form)fich4
      goto 310 
c
 317  continue
      write(42,*)'deplacement horizontal (amplitude)'
      print*,'deplacement horizontal (amplitude)'
      write (form,'(a2,i1,a4)')'(a',itai,',a6)'
      write(nomout,form) nom,'_depxy'
      itaio=itai+6
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif
      goto (3177,3178,3179)isor(ik)
 3177 write(fich4,form) nomout,'.1oc'
      nrec=npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3176
 3178 write(fich4,form) nomout,'.2oc'
      nrec=2*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3176
 3179 write(fich4,form) nomout,'.bin'
      nrec=4*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3176
 3176 continue
      vmax=hmax *100.
      vmin=hmin *100.
      kn=0 
      do j=npy,1,-1
         kn=kn+1
         read  (97,rec=j)(dx(i),i=1,npx)
         read  (98,rec=j)(dy(i),i=1,npx)
         do i=1,npx
           v(i)=dsqrt(dble(dx(i))**2+dble(dy(i))**2)*100.
         enddo
         goto (3171,3172,3173)isor(ik)
 3171    call wri1(v,ilo,kn)
         goto 3170
 3172    call wri2(v,ilo,kn)
         goto 3170
 3173    call wri4(v,ilo,kn) 
         goto 3170
 3170    continue
      enddo
      close (ilo)
      write(42,form)fich4
      goto 310
c
 318  continue
      write(42,*)'deplacement total (amplitude)'
      print*,'deplacement total (amplitude)'
      write (form,'(a2,i1,a4)')'(a',itai,',a6)' 
      write(nomout,form) nom,'_deptt' 
      itaio=itai+6
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif
      goto (3187,3188,3189)isor(ik)
 3187 write(fich4,form) nomout,'.1oc'
      nrec=npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3186
 3188 write(fich4,form) nomout,'.2oc'
      nrec=2*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3186
 3189 write(fich4,form) nomout,'.bin'
      nrec=4*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3186
 3186 continue
      vmax=tmax *100.
      vmin=tmin *100.
      kn=0 
      do j=npy,1,-1
         kn=kn+1
         read  (97,rec=j)(dx(i),i=1,npx)
         read  (98,rec=j)(dy(i),i=1,npx)
         read  (99,rec=j)(dz(i),i=1,npx)
        do i=1,npx
          v(i)=dsqrt(dble(dx(i))**2+dble(dy(i))**2+dble(dz(i))**2)*100.
        enddo
         goto (3181,3182,3183)isor(ik)
 3181    call wri1(v,ilo,kn)
         goto 3180
 3182    call wri2(v,ilo,kn)
         goto 3180
 3183    call wri4(v,ilo,kn) 
         goto 3180
 3180    continue
      enddo
      close (ilo)
      write(42,form)fich4
      goto 310  
c
 319   continue
      write(42,*)'deplacement projete sur une direction particuliere'
      print*,'deplacement projete sur une direction particuliere'   
      print*,'vecteur x y z'
      read(*,*)ts_x,ts_y,ts_z
      write (form,'(a2,i1,a4)')'(a',itai,',a6)'
      write(nomout,form) nom,'_parts'
      itaio=itai+6
      if(itaio.gt.9) then
        write (form,'(a2,i2,a4)')'(a',itaio,',a4)'
      else
        write (form,'(a2,i1,a4)')'(a',itaio,',a4)'
      endif
      goto (3197,3198,3199)isor(ik)
 3197 write(fich4,form) nomout,'.1oc'
      nrec=npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3196
 3198 write(fich4,form) nomout,'.2oc'
      nrec=2*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3196
 3199 write(fich4,form) nomout,'.bin'
      nrec=4*npx
      open(ilo,file=fich4,access='direct',recl=nrec)
      goto 3196
 3196 continue
      vmax=-1.0e20
      vmin=1.0e20
      kn=0
      do j=npy,1,-1
         kn=kn+1
         read  (97,rec=j)(dx(i),i=1,npx)
         read  (98,rec=j)(dy(i),i=1,npx)
         read  (99,rec=j)(dz(i),i=1,npx)
        do i=1,npx
           v(i)=dble(dx(i))*ts_x+dble(dy(i))*ts_y+dble(dz(i))*ts_z
           v(i)=-v(i)*100.
           vmax=dmax1(v(i),vmax)
           vmin=dmin1(v(i),vmin)
c          print*,dx(i),dy(i),dz(i),v(i)
        enddo
         goto (3191,3192,3193)isor(ik)
 3191    call wri1(v,ilo,kn)
         goto 3190
 3192    call wri2(v,ilo,kn)
         goto 3190
 3193    call wri4(v,ilo,kn)
         goto 3190
 3190    continue
c
      enddo
      close (ilo)
      write(42,form)fich4
      goto 310
c
 310  continue
      goto 300
      if(idep(ik).eq.4)then
        call cregmt1(nomout,form,fich4)
      elseif(idep(ik).eq.2) then
      else
        call cregmt(nomout,form,fich4)
      endif
 300  continue
      close (97)
      close (98)
      close (99) 
      open(97,file=ficht1)
      open(98,file=ficht2)
      open(99,file=ficht3)
      write(97,'(a1)')nul
      write(98,'(a1)')nul
      write(99,'(a1)')nul
      close (97)
      close (98)
      close (99) 
      print*,'apgmt'
c
c 
      write(42,'(A)') 'FIN' 
      print*,'**** FIN ****'
      CLOSE (42) 
c     CLOSE(6) 
c 
      STOP 
      END 
c-----------------Fin prog. principal --------------- 
c--------------------------------------------------------------------------------------------------------------------------------------------------------------------------
c-----------------Fin prog. principal ---------------
      subroutine cregmt1(nomout,form,fich4)
c
      parameter (nbc=256)
      IMPLICIT REAL*8(A-H,L,O-Z)
      character*80 nomout,form*20,fichgmt,fich4,fichgrd,fichcpt
      character v10*10,w10*10,script*300,ch20*20
      dimension ir(nbc),ig(nbc),ib(nbc),tmin(nbc),tmax(nbc)
      common/exit/npx,npy,vmin,vmax,xc1,yc1,pas
c
      write(fichcpt,form) nomout,'.cpt'
c
        xdeb=xc1
        ydeb=yc1
        xfin=xdeb+(npx-1)*pas
        yfin=ydeb+(npy-1)*pas
      call palet(ir,ig,ib,vmin,vmax,tmin,tmax)
      open(56,file=fichcpt)
      do i=1,256
         call cv10(tmin(i),v10)
         callcv10(tmax(i),w10)
         write(56,57)v10,ir(i),ig(i),ib(i),w10,ir(i),ig(i),ib(i)
 57    format (2(1x,a10,1x,i4,i4,i4))
       enddo
       close(56)
c
      write(fichgmt,form) nomout,'.gmt'
      write(fichgrd,form) nomout,'.grd'
c
      write(*,'(a)') form,fichgmt
c     ll=lnblnk(fichgmt)
      open(55,file=fichgmt)
c
      call tail (fich4,il)
      script(1:8)='xyz2grd '
      nd=8+1
      nf=nd+il-1
      script(nd:nf)=fich4(1:il)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -G'
      nd=nf+1
      nf=nd+il-1
      script(nd:nf)=fichgrd(1:il)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -I'
      nd=nf+1
      nf=nd+20-1
      call cv20(pas,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(pas,ch20)
      script(nd:nf)=ch20
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -R'
      nd=nf+1
      nf=nd+20-1
      call cv20(xdeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(xfin,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(ydeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(yfin,ch20)
      script(nd:nf)=ch20
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -b'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -Z'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -V'
      do i=nf+1,300
        script(i:i)=' '
      enddo
c
      write(55,'(a)')script(1:nf)
c
c  dessin image
c
      script(1:9)='grdimage '
      nd=9+1
      nf=nd+il-1
      script(nd:nf)=fichgrd(1:il)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -C'
      call tail (fichcpt,it)
      nd=nf+1
      nf=nd+it-1
      script(nd:nf)=fichcpt(1:it)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -R'
      nd=nf+1
      nf=nd+20-1
      call cv20(xdeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(xfin,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(ydeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(yfin,ch20)
      script(nd:nf)=ch20
      ex=xfin-xdeb
      ey=yfin-ydeb
      if(ey.gt.ex) then
        ech1=18/ex
        ech2=24/ey
        ech=dmin1(ech1,ech2)
        nd=nf+1  
        nf=nd+2
        script(nd:nf)=' -P'
      else
        ech1=24/ex
        ech2=18/ey
        ech=dmin1(ech1,ech2)
      endif
      nd=nf+1
      nf=nd+3
      call cv20(ech,ch20)
      script(nd:nf)=' -Jx'
      nd=nf+1
      nf=nd+20-1
      script(nd:nf)=ch20
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -K'
      nd=nf+1
      nf=nd+4
      script(nd:nf)='>i.ps'
      do i=nf+1,300
        script(i:i)=' '
      enddo
      write(55,'(a)')script(1:nf)
c
c dessin trace faille
c
      script(1:5)='psxy '
      nd=5+1
      nf=nd+9
      script(nd:nf)='trfault.xy'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -R'
      nd=nf+1
      nf=nd+3
      script(nd:nf)=' -Jx'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -O'
      nd=nf+1
      nf=nd+5
      script(nd:nf)=' -F255'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -M'
      nd=nf+1
      nf=nd+5
      script(nd:nf)='>>i.ps'
      write(55,'(a)')script(1:nf)
      close(55)
c
      return
      end

c

c
      subroutine cregmt(nomout,form,fich4)
c
      parameter (nbc=256)
      IMPLICIT REAL*8(A-H,L,O-Z)
      character*80 nomout,form*20,fichgmt,fich4,fichgrd,fichcpt
      character v10*10,w10*10,script*300,ch20*20
      dimension ir(nbc),ig(nbc),ib(nbc),tmin(nbc),tmax(nbc)
      common/exit/npx,npy,vmin,vmax,xc1,yc1,pas
c
      write(fichcpt,form) nomout,'.cpt'
c
        xdeb=xc1
        ydeb=yc1
        xfin=xdeb+(npx-1)*pas
        yfin=ydeb+(npy-1)*pas
      call palet(ir,ig,ib,vmin,vmax,tmin,tmax)
      open(56,file=fichcpt)
      do i=1,256
         call cv10(tmin(i),v10)
         callcv10(tmax(i),w10)
         write(56,57)v10,ir(i),ig(i),ib(i),w10,ir(i),ig(i),ib(i)
 57    format (2(1x,a10,1x,i4,i4,i4))
       enddo
       close(56)
c
      write(fichgmt,form) nomout,'.gmt'
      write(fichgrd,form) nomout,'.grd'
c
      write(*,'(a)') form,fichgmt
c     ll=lnblnk(fichgmt)
      open(55,file=fichgmt)
c
      call tail (fich4,il)
      script(1:8)='xyz2grd '
      nd=8+1
      nf=nd+il-1
      script(nd:nf)=fich4(1:il)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -G'
      nd=nf+1
      nf=nd+il-1
      script(nd:nf)=fichgrd(1:il)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -I'
      nd=nf+1
      nf=nd+20-1
      call cv20(pas,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(pas,ch20)
      script(nd:nf)=ch20
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -R'
      nd=nf+1
      nf=nd+20-1
      call cv20(xdeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(xfin,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(ydeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(yfin,ch20)
      script(nd:nf)=ch20
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -b'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -Z'
      nd=nf+1 
      nf=nd+2
      script(nd:nf)=' -V'
      do i=nf+1,300 
        script(i:i)=' ' 
      enddo 
c
      write(55,'(a)')script(1:nf)
c
c  dessin image
c
      script(1:9)='grdimage '
      nd=9+1
      nf=nd+il-1
      script(nd:nf)=fichgrd(1:il)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -C'
      call tail (fichcpt,it)
      nd=nf+1
      nf=nd+it-1
      script(nd:nf)=fichcpt(1:it)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -R'
      nd=nf+1
      nf=nd+20-1
      call cv20(xdeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(xfin,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(ydeb,ch20)
      script(nd:nf)=ch20
      script(nf+1:nf+1)='/'
      nd=nf+2
      nf=nd+20-1
      call cv20(yfin,ch20)
      script(nd:nf)=ch20
      script(nd:nf)=ch20 
      ex=xfin-xdeb 
      ey=yfin-ydeb
      if(ey.gt.ex) then 
        ech1=18/ex 
        ech2=24/ey 
        ech=dmin1(ech1,ech2) 
        nd=nf+1   
        nf=nd+2 
        script(nd:nf)=' -P' 
      else 
        ech1=24/ex 
        ech2=18/ey 
        ech=dmin1(ech1,ech2) 
      endif 
      nd=nf+1
      nf=nd+3
      call cv20(ech,ch20)
      script(nd:nf)=' -Jx'
      nd=nf+1
      nf=nd+20-1
      script(nd:nf)=ch20 
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -K'
      nd=nf+1
      nf=nd+4
      script(nd:nf)='>i.ps'
      do i=nf+1,300 
        script(i:i)=' ' 
      enddo 
      write(55,'(a)')script(1:nf)
c
c   contour
c
      print*,"contour"
      eca=vmax-vmin
      ee=1
      if(eca.eq.0.d0) eca=0.01
 58   continue
      if(eca.lt.1) then
        eca=eca*10
        ee=ee/10
        goto 58
      elseif (eca.gt.10) then
        eca=eca/10
        ee=ee*10
        goto 58
      else
        continue
        ee=ee/2
      endif
      script(1:11)='grdcontour '
      print*,script
      nd=11+1
      nf=nd+il-1
      script(nd:nf)=fichgrd(1:il)
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -R'
      nd=nf+1
      nf=nd+3
      script(nd:nf)=' -Jx'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -O'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -K'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -A'
      if (idep.eq.4)then
         ee=254
      endif
      call cv10(ee,w10)
      nd=nf+1
      nf=nd+10-1
      script(nd:nf)=w10
      nd=nf+1
      nf=nd+6
      script(nd:nf)=' >>i.ps'
      print*,script
      write(55,'(a)')script(1:nf)
c           
c dessin trace faille
c
      script(1:5)='psxy '
      nd=5+1
      nf=nd+9
      script(nd:nf)='trfault.xy'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -R'
      nd=nf+1
      nf=nd+3
      script(nd:nf)=' -Jx'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -O'
      nd=nf+1
      nf=nd+5
      script(nd:nf)=' -F255'
      nd=nf+1
      nf=nd+2
      script(nd:nf)=' -M'
      nd=nf+1
      nf=nd+5
      script(nd:nf)='>>i.ps'
      write(55,'(a)')script(1:nf)
      close(55)
c
      return
      end
c
c-----------------------------------------
      SUBROUTINE UTMGEO(XX,YY,LON,LAT)
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,L,O-Z)
      COMMON/COMOU/CC,E2,PI,IFU,AK,KH
c
      trg=200.d0/180.d0 
      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)
      L=P*KH
      S=E2*CO*CO
      R=CC/DSQRT(1.d0+S)
      PP=P*KH
      YR=YR-ARCME(L)
      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=L+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(L)*E2*(G-L)*DCOS(L)
      LON=AM+AM0
      LAT=L+V*(G-L)
      LON=(AM+AM0)/RD
      LAT=LAT/RD*KH
      lon=lon/trg
      lat=lat/trg
c
      RETURN
      END
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/COMOU/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 ------------------------------------------------
       FUNCTION ARCME(AL)
C--------------------------
C         Fonction Arc de mÇridien en fonction de la latitude :
c
      IMPLICIT REAL*8(A-H,L,O-Z)
      COMMON/COMOU/CC,E2,PI,IFU,AK,KH
c
      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--------------------------------------------------------------------------
      SUBROUTINE DEFOR(NFA,XX,YY,DFX,DFY,DFZ)
C-----------------------------------------
C           Calcul de la deformation due a un ensemble de failles
c
      parameter (nbf=200)
      IMPLICIT REAL*8(A-H,L,O-Z)
      COMMON/F1/XF(nbf),YF(nbf),FF(nbf),H(nbf),W(nbf),R(nbf),TF(nbf)
      COMMON/F2/US(nbf),UD(nbf),UT(nbf)
c
      A=0.
      B=0.
      C=0.
c
      DO 200 I=1,NFA
       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 ---------Calcul des deplacements pour chaque faille :
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)
       A=A+A1
       B=B+B1
       C=C+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)
       A=A+A1
       B=B+B1
       C=C+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)
       A=A+A1
       B=B+B1
       C=C+C1
   53 CONTINUE
c
  200 CONTINUE
      DFX=A
      DFY=B
      DFZ=C
c   
      RETURN
      END
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 X1,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----------------- FIN -----------------
c
      subroutine cv20(val,ch20)
c
      implicit real*8 (a-h,l,o-z)
      character*20 ch20
c
      write(ch20,'(e20.13)')val
      if(val.ge.0.d0) then
        ch20(1:1)='+'
      endif
      ch20(17:17)='e'
      return
      end
c
c
      subroutine cv10(val,ch10)
      character*10 ch10
      real*8 val 
c
      write(ch10,'(e10.3)')val
      if(val.ge.0.d0) then
        ch10(1:1)='+'
      endif
      ch10(7:7)='e'
      return
      end
c 
c
c------------------------------------------
c
      subroutine tail(fich,ii)
c
      character*80 fich
c
      ii=lnblnk(fich)
      return
      end
c
c
c------------------------------------------
c
      subroutine palet (ir,ig,ib,zmin,zmax,tab1,tab2)
c
      implicit real*8(a-h,l,o-z)
      parameter (nbc=256)
      dimension ir(nbc),ig(nbc),ib(nbc)
      dimension tab1(nbc),tab2(nbc)
c
         do i=1,103
            ir(i)=0
         enddo
         do i=104,154
           ir(i)=ir(i-1)+5
         enddo
         do i=155,256
            ir(i)=255
         enddo
c
         ig(1)=0
         do i=2,52
           ig(i)=ig(i-1)+5
         enddo
         do i=53,154
            ig(i)=255
         enddo
         do i=155,205
            ig(i)=ig(i-1)-5
         enddo
         do i=206,256
            ig(i)=0
         enddo
c
         do i=1,52
           ib(i)=255
         enddo
         do i=53,103
          ib(i)=ib(i-1)-5
         enddo
         do i=104,205
          ib(i)=0
         enddo
         do i=206,256
            ib(i)=ib(i-1)+5
         enddo
c
         ecart=(zmax-zmin)/256.d0
         do i=1,256
              tab1(i)=zmin+(ecart*(i-1))
              tab2(i)=zmin+(ecart*i)
         enddo
         return
         end
c
c--------------------------------
c
         subroutine wri1(v,ilo,j)
c
      parameter (nbpc=3000)
      IMPLICIT REAL*8 (A-H,L,O-Z)
      dimension v(nbpc)
      character*1 tab(nbpc)
      common/exit/npx,npy,vmin,vmax,xc1,yc1,pas
c
        do i=1,npx
           v(i)=(v(i)-vmin)*255/(vmax-vmin)
            tab(i)=char(nint(v(i)))
        enddo
        write(ilo,rec=j)(tab(i),i=1,npx)
      return
      end
c--------------------------------
c
         subroutine wri2(v,ilo,j)
c
      parameter (nbpc=3000)
      IMPLICIT REAL*8 (A-H,L,O-Z)
      dimension v(nbpc)
      integer*2 tab(nbpc)
      common/exit/npx,npy,vmin,vmax,xc1,yc1,pas
c
        do i=1,npx
           v(i)=(v(i)-vmin)*255/(vmax-vmin)
            tab(i)=nint(v(i))
        enddo
        write(ilo,rec=j)(tab(i),i=1,npx)
      return 
      end 
c
c-------------------------------- 
c 
         subroutine wri4(v,ilo,j)
c 
      parameter (nbpc=3000)
      IMPLICIT REAL*8 (A-H,L,O-Z) 
      dimension v(nbpc)
      real*4  tab(nbpc) 
      common/exit/npx,npy,vmin,vmax,xc1,yc1,pas
c 
        do i=1,npx 
           tab(i)=sngl(v(i))
        enddo 
        write(ilo,rec=j)(tab(i),i=1,npx) 
      return  
      end   
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 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-----------------------------------------
cc
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
        slip=dsqrt(us**2+ud**2)
        rake=-datan2d(ud,us)
c
        return
        end
c
c----------------------------------------------
        subroutine planvec(taba,tabb,tabc) 
c
        common/plan/px(4),py(4),tsx(4),tsy(4),tsz(4),npp
c
      implicit real*8(a-h,l,o-z)
        dimension taba(1),tabb(1),tabc(1)
c
      write(42,*)' '
      write(42,*)'+++++++++++++++++++++++++++++++++++++++++++++++'
      write(42,*)'vecteurs et plans  de projection'
      do i=1,npp
         write(42,*)px(i),py(i),tsx(i),tsy(i),tsz(i)
      enddo
        sax=0.   
        sbx=0.
        scx=0.
        say=0.
        sby=0.
        scy=0.
        saz=0.
        sbz=0.
        scz=0.
        nc=0
        do i=1,npp
          do j=i+1,npp
            do k=j+1,npp
        write(42,*)i,j,k
c x
      bx=((tsx(i)-tsx(k))-(((tsx(j)-tsx(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=((tsx(i)-tsx(k))-bx*(py(i)-py(k)))/(px(i)-px(k))
          cx=tsx(i)-ax*px(i)-bx*py(i)
          sax=sax+ax 
          sbx=sbx+bx 
          scx=scx+cx 
c y
      by=((tsy(i)-tsy(k))-(((tsy(j)-tsy(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=((tsy(i)-tsy(k))-by*(py(i)-py(k)))/(px(i)-px(k))
          cy=tsy(i)-ay*px(i)-by*py(i)
          say=say+ay
          sby=sby+by 
          scy=scy+cy
c z
      bz=((tsz(i)-tsz(k))-(((tsz(j)-tsz(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=((tsz(i)-tsz(k))-bz*(py(i)-py(k)))/(px(i)-px(k))
          cz=tsz(i)-az*px(i)-bz*py(i)
          saz=saz+az 
          sbz=sbz+bz 
          scz=scz+cz
          write(42,*)ax,bx,cx
          write(42,*)ay,by,cy
          write(42,*)az,bz,cz
          nc=nc+1
            enddo
          enddo
        enddo
          taba(1)=sax/npp
          tabb(1)=sbx/npp
          tabc(1)=scx/npp
          taba(2)=say/npp
          tabb(2)=sby/npp
          tabc(2)=scy/npp
          taba(3)=saz/npp
          tabb(3)=sbz/npp
          tabc(3)=scz/npp

c
        write (42,*)'plan x ',taba(1),tabb(1),tabc(1)
        write (42,*)'plan y ',taba(2),tabb(2),tabc(2)
        write (42,*)'plan z ',taba(3),tabb(3),tabc(3)
c
        print*,'TEST'
        do k=1,npp
           tx=taba(1)*px(k)+tabb(1)*py(k)+tabc(1)
           ty=taba(2)*px(k)+tabb(2)*py(k)+tabc(2) 
           tz=taba(3)*px(k)+tabb(3)*py(k)+tabc(3) 
           write(42,*)k
           write(42,*)tsx(k),tx
           write(42,*)tsy(k),ty
           write(42,*)tsz(k),tz
        enddo
      write(42,*)'+++++++++++++++++++++++++++++++++++++++++++++++' 
      write(42,*)''
c
        return
        end
