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