c c programme inversion par gradient c avec pregradient integre et ponderation des fichiers de donnnees c c janv2000 - essaie de fixer des patchs a des valeurs fixees par modele a priori c mars 2000 - augmentation de delta tous les 10 iterations c 31 mars 2000 modif pour fichier de sortie lorsque nrec!=1 (*alea(i,j)) c implicit real*8 (a-h,o-z) parameter (mxd=50000,mxs=200,maxrec=10,mxf=9) cc mxs=maxi de sources, mxf=maxi de fichiers de donnees, mxd=maxi de donnees character*100 fires,form,formf,formo,fida,fimo character*100 nom1,nom2,fisou,fault,nom3,nom4,nom5 c dimension a(mxd,mxs),grad(mxs) dimension bs(mxd),x(mxs) dimension alea(maxrec,mxs),io(maxrec) dimension bests(maxrec),bestx(maxrec,mxs) dimension at(mxs,mxd),tapar(mxs,10),libre(mxs) dimension xdata(mxd),ydata(mxd),ddata(mxd) dimension weight(mxf),edata(mxd) dimension nv(mxf,2),irv(mxf) dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2) dimension ixxx(mxs,2),iyyy(mxs,2) dimension xf(mxs),yf(mxs),ff(mxs),h(mxs) dimension w(mxs),r(mxs),tf(mxs) dimension us(mxs),ud(mxs),ut(mxs) dimension t_momu(mxs) real*4, rand_unif c c common/param/tapar,numsour,libre common/data1/xdata,ydata,ddata,edata,weight common/data2/nv,tsa,tsb,tsc,irv,dref common/data3/nfic common/defo/a COMMON/F1/XF,YF,FF,H,W,R,TF COMMON/F2/US,UD,UT COMMON/COMOB/Q,ST,CT,ELAS,U0 COMMON/COMOG/CC,E2,PI,IFU,AK,KH common/grid_mom/ixxx,iyyy c pi=4*datan(1.d0) smu = 3.0d10 elas=0.5d0 cc=6399593.6258000d0 AK=0.9996d0 eprim=0.08209443779496d0 EP2=EPRIM*EPRIM E2=EP2 c idimax=100 idim=0 c open(99,file='GRAD.LOG') call pregradient (ndat,nstri,ndip,igeo) print*,'FIN CALCUL DIRECT' c c numx=nstri numy=ndip npar=numsour c do 50 jj=1,nfic do 51 kk=nv(jj,1),nv(jj,2) xx=xdata(kk)/1000. yy=ydata(kk)/1000. if (igeo.eq.2) then call utmgeo (xx,yy,xx,yy) endif do ll=1,npar x(ll)=1.d0 enddo 51 enddo 50 enddo c c call mvmult(mxd,mxs,ndat,npar,a,x,bs) do k=1,npar t_momu(k)=tapar(k,8) enddo print*,ndat,' donnees ; ',npar,' parametres' write(99,*)ndat,' donnees ; ',npar,' parametres' c write(99,*)' SOURCES' print*,'SOURCES' do kd=1,numy k1=(kd-1)*numx+1 k2=k1+numx-1 write(99,*)(real(tapar(j,8)),j=k1,k2) write(*,*)(real(tapar(j,8)),j=k1,k2) enddo c print*,'nb de recherches aleatoires <',maxrec read(*,*)nrec write(99,*)'nb de recherches aleatoires <',maxrec nrec=nrec+1 c print*,' glissement min, max, increment (cm)' read(*,*)smin,smax,sinc write(99,*)' glissement min, max, increment (cm)' write(99,*)sngl(smin),sngl(smax),sngl(sinc) c print*,'nb d"iterations pour gradient' read(*,*)niter write(99,*)' nb d"iterations pour gradient',niter c k=1 do j=1,npar alea(k,j)=1.d0 enddo c c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ c irec=0 c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ BOUCLE SUR les REC 300 continue irec=irec+1 write(99,*)'' write(99,*)'' write(99,*)'+++++++++++++++++++' write(99,*)' IREC ',irec idim=0 delta=1 c c print*,'ALEA' do 13 j=1,npar c print*,alea(irec,j) x(j)=alea(irec,j) bestx(irec,j)=x(j) 13 enddo c c calcul de l'erreur quadratique c call mvmult(mxd,mxs,ndat,npar,a,x,bs) ssr=0.d0 sw=0.d0 c call calcssr(ndat,ddata,bs,ssr) write(99,*)'' write(99,*)' CALCUL DE L ERREUR INITIALE' write(99,*)'no fich, rms fich, rms tout' do if=1,nfic n1=nv(if,1) n2=nv(if,2) call calcssrt(ddata,bs,n1,n2,ssrs) ssr=ssr+ssrs*weight(if) sw=sw+weight(if)*(n2-n1+1) ww=weight(if) print*,if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/(sw)) write(99,*)if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/(sw)) enddo print*,'RMS TOT',dsqrt(ssr/sw) write(99,*)'RMS TOT ',dsqrt(ssr/sw) c bests(irec)=ssr c write(99,*) ' ' print*,'DELTA INI ',delta write(99,*)'DELTA INI ',real(delta) write(*,'(1000f8.1)')(x(j)*tapar(j,8),j=1,npar) write(99,'(1000f8.1)')(real(x(j)*tapar(j,8)),j=1,npar) print*,'PAR INI' do j=1,npar print*,x(j) enddo c c call mtxtrp(mxd,mxs,ndat,npar,a,at) c 12 continue c increase_delta=0 c c****************** c*** ITER ******* c****************** do 100 iter=1,niter increase_delta=increase_delta+1 if (increase_delta.gt.10) then print*,'INCREASE DELTA ', delta*2. delta=delta*2. increase_delta=0 endif c print*,'' print*,'ITER',iter write(99,*)'' write(99,*)'ITER',iter call findgrad(mxd,mxs,ndat,npar,a,at,ddata,x,grad) c do 200 j=1,npar if (libre(j).ne.0) then x(j)=x(j)+delta*grad(j) endif slipr=x(j)*tapar(j,8) if(slipr.lt.0.0) x(j)=0.d0 200 continue call mvmult(mxd,mxs,ndat,npar,a,x,bs) ssr=0.d0 sw=0.d0 write(99,*)'no fich, rms fich, rms tout' do if=1,nfic n1=nv(if,1) n2=nv(if,2) call calcssrt(ddata,bs,n1,n2,ssrs) ssr=ssr+ssrs*weight(if) sw=sw+weight(if)*(n2-n1+1) ww=weight(if) c print*,if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/sw) write(99,*)if,dsqrt(ssrs/(n2-n1+1)),dsqrt(ssr/sw) enddo print*,'RMS TOT',dsqrt(ssr/sw) write(99,*)'RMS TOT ',dsqrt(ssr/sw) c c write(99,*)' SOURCES' c do kd=1,numy c k1=(kd-1)*numx+1 c k2=k1+numx-1 c write(99,*)(real(x(j)*tapar(j,8)),j=k1,k2) c write(99,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2) c write(*,*)(real(x(j)*tapar(j,8)),j=k1,k2) c write(*,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2) c enddo c if (ssr.lt.bests(irec)) then brms=dsqrt(ssr/sw) bests(irec)=ssr do j=1,npar bestx(irec,j)=x(j) enddo else idim=idim+1 print*,'idim DIMIN DELTA',idim,delta*.5 write(99,*)'idim,DIMIN DELTA', idim,delta*.5 increase_delta=0 ssr=bests(irec) delta=delta*.5 do j=1,npar x(j)=bestx(irec,j) enddo iter=iter-1 endif c c if(idim.gt.idimax) then write(99,*)' MAXIMUM DIMINUTION' write(99,*)' idimax idim',idimax,idim print*,' MAXIMUM DIMINUTION' print*,'idimax ',idim goto 101 endif 100 continue 101 continue print*,iter, ' iterations' write (99,*)'' write (99,*)' NB ITERATIONS EFFECTUEES', iter c cccc c998 continue c brms=bests(irec)/sw c print*,' essais autres valeurs de slip' c do j=1,npar c print*,real(bestx(irec,j)*tapar(j,8)) c print*,' slip ?' c read(*,*)slip c x(j)=slip/tapar(j,8) c enddo c call mvmult(mxd,mxs,ndat,npar,a,x,bs) c do if=1,nfic c n1=nv(if,1) c n2=nv(if,2) c call calcssrt(ddata,bs,n1,n2,ssrs) c ssr=ssr+ssrs*weight(if) c sw=sw+weight(if)*(n2-n1+1) c ww=weight(if) c enddo c print*,'slip RMS ',slip, dsqrt(ssr/sw) c print*,'BEST', (bestx(irec,1)*tapar(1,8)),brms c if(ssr.lt.bests(irec)) then c bests(irec)=ssr c brms=dsqrt(ssr/sw) c do j=1,npar c bestx(irec,j)=x(j) c enddo c endif c goto 998 c c c----------------------meilleur modele c print*,'MEILLEUR MODELE' print*,bests(irec),dsqrt(bests(irec)/sw) write(*,'(1000f8.1)')(bestx(irec,j)*tapar(j,8),j=1,npar) c write(99,*)'' write(99,*)'meilleur modele' write(99,*)'RMS ',real(dsqrt(bests(irec)/sw)) write(99,*)' SOURCES' print*,' SOURCES' do kd=1,numy k1=(kd-1)*numx+1 k2=k1+numx-1 write(99,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2) write(*,*)(real(bestx(irec,j)*tapar(j,8)),j=k1,k2) enddo c c c-----------------------recherche modele aleatoire print*,'' if(irec.ge.nrec) then goto 400 endif do k=1,npar tmax=smax/tapar(k,8) tmin=smin/tapar(k,8) tdif=tmax-tmin tinc=sinc/tapar(k,8) ns=idnint(tdif/tinc)+1 alea(irec+1,k)=tmin+dble(nint(ns*rand_unif(0)))*tinc enddo goto 300 c 400 continue print*,'nrec ',nrec c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ BOUCLE SUR les REC c c------------------------------------------------- c-------------tri et ecriture des modeles aleatoires recherches c do k=1,nrec io(k)=k enddo c call tri (bests,io,nrec) c c open(51,file='resu') nw=npar*4 open(52,file='resu.bin',access='direct',recl=nw) c c write(51,*)'nb modeles sources ',nrec,npar do k=1,nrec write(51,*)bests(k),dsqrt(bests(irec)/ndat) write(51,'(1000f8.1)')(bestx(k,j)*tapar(j,8),j=1,npar) c write(*,'(1000f8.1)')(bestx(k,j)*tapar(j,8),j=1,npar) write(52,rec=k)(sngl(bestx(k,j)*tapar(j,8)),j=1,npar) enddo close(51) close(52) c nom1='resi' nom2='sour' nom3='out' nom4='data' nom5='modl' c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c write(99,*)'' write(99,*)'*********** SORTIE ****************' c do 60 k=1,nrec do j=1,npar tapar(j,8)=t_momu(j)*bestx(k,j)*alea(k,j) enddo c call direct do j=1,npar x(j)=1.0 enddo call mvmult(mxd,mxs,ndat,npar,a,x,bs) c if(k.lt.10) then form='(a4,i1,a4)' formf='(a3,i1,a6)' formo='(a4,i1,i1,a4)' elseif(k.lt.100) then form='(a4,i2,a4)' formf='(a3,i2,a6)' formo='(a4,i2,i1,a4)' else form='(a4,i3,a4)' formf='(a3,i3,a6)' formo='(a4,i2,i1,a4)' endif c c rmst=0.d0 rmsp=0.d0 sw=0.d0 do 61 jj=1,nfic write(99,*)'' print*,'FICHIER',jj write(99,*)'FICHIER',jj write(fida,formo)nom4,k,jj,'.out' write(fimo,formo)nom5,k,jj,'.out' write(fires,formo)nom1,k,jj,'.out' open(56,file=fida) open(57,file=fimo) open(53,file=fires) c rms=0.0 nndd=0 write(99,*)' reference donnees modl',dref(jj,1),dref(jj,2) c do 62 kk=nv(jj,1),nv(jj,2) nndd=nndd+1 if(igeo.eq.2) then call utmgeo (xdata(kk)/1000.,ydata(kk)/1000.,xx,yy) else xx=xdata(kk)/1000. yy=ydata(kk)/1000. endif deff=bs(kk)+dref(jj,2) ddaa=ddata(kk)+dref(jj,2) adiff=ddaa-deff write(56,*)xx,yy,float(ddaa) write(57,*)xx,yy,float(deff) write(53,*)xx,yy,float(adiff) rms=adiff**2+rms 62 enddo c c rmsp=rms*weight(jj)+rmsp sw=sw+weight(jj)*(nv(jj,2)-nv(jj,1)+1) rmst=rms+rmst c c rms=dsqrt(rms/dble(nndd)) write(53,*)' RMS= ',rms write(99,*)'ndata, RMS, poids ',nndd,rms,weight(jj) print*,'ndata RMS= ',nndd,rms c c close(53) close(56) close(57) 61 enddo c c rmst=dsqrt(rmst/dble(ndat)) rmsp=dsqrt(rmsp/sw) print*,'rms tot',rmst write(99,*)'' write(99,*)' ndat RMS TOTAL ',ndat,rmst write(99,*)' RMS PONDERE ',rmsp c c------- fichier sour_k_.out (source en binaire) write(fisou,form)nom2,k,'.out' print*,fisou nw=4*numy open(54,file=fisou,access='direct',recl=nw) l=0 do 73 j=numx,1,-1 l=l+1 c write(*,'(1000f7.1)')(*tapar(i,8), c & i=j,j+numx*(numy-1),numx) write(*,*)(real(tapar(i,8)) & ,i=j,j+numx*(numy-1),numx) write(54,rec=l)(real(tapar(i,8)) & ,i=j,j+numx*(numy-1),numx) 73 enddo close(54) 60 enddo c c----------------------------Ecriture sources---- c do 72 i=1,npar call jcr2new(tapar(i,1),tapar(i,2),tapar(i,5), & tapar(i,6),tapar(i,3),tapar(i,4)) if(igeo.eq.2) then call utmgeo(tapar(i,1),tapar(i,2),tapar(i,1),tapar(i,2)) endif 72 enddo c do 71 k=1,nrec c print*,"RECHERCHE ",K c do 74 i=1,npar tapar(i,8)=t_momu(i)*bestx(k,i)*alea(k,i) c print*,t_momu(i),bestx(k,i),tapar(i,8) c print*,tapar(i,1),tapar(i,2) 74 enddo c write(fault,formf)nom3,k,'.fault' print*,fault open(55,file=fault) qq=0.d0 write(55,*)npar,' 5 0 0 ' do i=1,npar c write(55,*)i write(55,*)tapar(i,1),qq,tapar(i,2),qq write(55,*)tapar(i,3),qq,tapar(i,4),qq write(55,*)tapar(i,5),qq,tapar(i,6),qq,tapar(i,7),qq write(55,*)tapar(i,8)+0.001,qq,tapar(i,9),qq,tapar(i,10),qq enddo close(55) c c------- fichier sour_k_.out (source en binaire) c c write(fisou,form)nom2,k,'.out' c print*,fisou c nw=4*numy c open(54,file=fisou,access='direct',recl=nw) c l=0 c do 73 j=numx,1,-1 c l=l+1 c write(*,'(1000f7.1)')(tapar(i,8), c & i=j,j+numx*(numy-1),numx) c write(*,*)(tapar(i,8) c & ,i=j,j+numx*(numy-1),numx) c write(54,rec=l)(real(tapar(i,8)) c & ,i=j,j+numx*(numy-1),numx) c73 enddo c close(54) c 71 enddo c c----------------------- stop end c c--------------------------------- subroutine mtxtrp(mxd,mxs,ndat,npar,a,at) c implicit real*8 (a-h,o-z) dimension at(mxs,mxd) dimension a(mxd,mxs) c do i=1,ndat do j=1,npar at(j,i)=a(i,j) enddo enddo return end c------------------------------------ c subroutine findgrad(mxd,mxs,ndat, & npar,a,at,ddata,x,grad) c implicit real*8 (a-h,o-z) dimension a(mxd,mxs),at(mxs,mxd) dimension x(mxs),ddata(mxd) dimension grad(mxs) dimension passa(100000) c do 10 i=1,ndat sum=0.d0 do 20 j=1,npar sum=sum+a(i,j)*x(j) 20 continue passa(i)=ddata(i)-sum c print*,'passa ',passa(i) 10 continue call mvmult(mxs,mxd,npar,ndat,at,passa,grad) return end c------------------------------------ c subroutine mvmult(max1,max2,n1,n2,a,v,r) c implicit real*8 (a-h,o-z) parameter (mxd=50000,mxs=200) c dimension a(max1,max2),v(1),r(1) c do 10 i=1,n1 sum=0.d0 do 20 j=1,n2 sum=sum+a(i,j)*v(j) 20 continue r(i)=sum 10 continue return end c c----------------------------------------- c subroutine calcssr(ndat,b,bs,ssr) c implicit real*8 (a-h,o-z) dimension b(1),bs(1) c ssr=0.d0 do 10 i=1,ndat dif=bs(i)-b(i) ssr=ssr+dif**2 10 continue return end c c----------------------------------------- c subroutine calcssrt(b,bs,n1,n2,ssr) c implicit real*8 (a-h,o-z) c dimension bs(1),b(1) c c ssr=0.d0 do 10 i=n1,n2 dif=bs(i)-b(i) c PRINT*,bs(i),b(i) ssr=ssr+dif**2 10 continue return end c c------------------------------------------------ c c subroutine ecr_par(igeo,numx,numy) c c implicit real*8 (a-h,o-z) c c parameter (mxs=200) c parameter (mxf=9) c c character*100 fpar c dimension tapar(mxs,10) c c common/param/tapar,numsour c c print*,'nom du fichier sortie parametre de faille' c read(*,'(a)')fpar c open (40,file=fpar) c c nnuummss=0 c write(40,*)numsour,' 5 0 0 ' c do k=1,numsour c nnuummss=nnuummss+1 c write(40,*)nnuummss c write(40,*)tapar(k,1),zz,tapar(k,2),zz c write(40,*)tapar(k,3),zz,tapar(k,4),zz c write(40,*)tapar(k,5),zz,tapar(k,6),zz,tapar(k,7),zz c write(40,*)tapar(k,8),zz,tapar(k,9),zz,tapar(k,10),zz c enddo c c c close(40) c return c end c c-------------------------------------------------------- c subroutine maxv (mxs,npar,grad,vmin,vmax) c c implicit real*8 (a-h,o-z) dimension grad(mxs) c vmin=1e50 vmax=-1e50 do j=1,npar vmin=dmin1(vmin,grad(j)) vmax=dmax1(vmax,grad(j)) enddo return end c c C**************************************************************************** c SUBROUTINE RAND_INIT c integer time, isec c real*4 dtime,tarray(2),adum c INTEGER*4 I,J c COMMON /RANDCOM/ INITFLAG,RANTAB,INDEX1,INDEX2,INDEX3 LOGICAL INITFLAG INTEGER*4 RANTAB(0:54) INTEGER*4 INDEX1,INDEX2,INDEX3 DATA INITFLAG /.FALSE./ c c get the year,month,day,hour,min,sec from time() and calculate a seed c based on this value [secs] (relative to 1970) plus the time elapsed for c executing a 'waiting' loop - to make sure that two subsequent calls to c RAND_INIT don't yield the same seed c isec = int(dtime(tarray)) isec = time() adum = 0.01 do j = 1,1000000 adum = atan(atan(atan(atan(atan(adum))))) enddo isec = isec + int(dtime(tarray)*1.e6) RANTAB(0) = isec c C C Use a linear congruential psuedo random number generator to C initialize the table for the RAND_$UNIF function. Note that C Knuth's rules for 'good' generators requires a multiplier of C the form ending in ...x21 where 'x' is an even digit. The C choice of 31415821 is taken from one of his examples. C DO 100 I = 1,54 RANTAB(I) = RANTAB(I-1)*31415821+1 100 CONTINUE c INITFLAG = .TRUE. INDEX1 = 0 INDEX2 = 23 INDEX3 = 54 RETURN END C**************************************************************************** **************************************************************************** C C Uniform Random Number Generator. C Generates a psuedo-random number in the range 0.0 to 1.0 with C a uniform distribution. The first time RAND_$UNIF is called the C argument is used to seed the random number generator if the C random number generator has not already been seeded by a call C to the RAND_$INIT routine. If the seed passed on the first call C is greater than 0, then the seed is used to initialize the C random number generator, otherwise a built in seed number is C used to set up the initial table. Since IEEE standard floating C point numbers have a 24 bit mantissa (one bit being the hidden bit) C we generate random integers in the range 0 to 2**24-1 and then C divide them by 2**24-1 to get a real number in the range 0 to C 1.0, inclusive. C C**************************************************************************** REAL*4 FUNCTION RAND_UNIF (ISEED) COMMON /RANDCOM/ INITFLAG,RANTAB,INDEX1,INDEX2,INDEX3 LOGICAL INITFLAG INTEGER*4 RANTAB(0:54) INTEGER*4 INDEX1,INDEX2,INDEX3 INTEGER*4 ISEED C C If random number generator has not already been seeded C by calling the RAND_$INIT routine, then set up the initial C table for the additive congruential method used here. If C the seed number supplied by the caller is 0, then use our C own seed number. C IF (INITFLAG) GOTO 100 IF (ISEED.GT.0) THEN RANTAB(0) = ISEED ELSE RANTAB(0) = 31415926 ENDIF C C Use a linear congruential psuedo random number generator to C initialize the table for the RAND_$UNIF function. Note that C Knuth's rules for 'good' generators requires a multiplier of C the form ending in ...x21 where 'x' is an even digit. The C choice of 31415821 is taken from one of his examples. C DO 10 I = 1,54 RANTAB(I) = RANTAB(I-1)*31415821+1 10 CONTINUE INITFLAG = .TRUE. INDEX1 = 0 INDEX2 = 23 INDEX3 = 54 C C Generate the next random integer in the range 0 to 2**24-1 C and then convert it into a real number in the range 0.0 to 1.0 C 100 INDEX1 = MOD(INDEX1+1,55) INDEX2 = MOD(INDEX2+1,55) INDEX3 = MOD(INDEX3+1,55) RANTAB(INDEX1) = RANTAB(INDEX2)+RANTAB(INDEX3) RAND_UNIF = FLOAT(IAND(RANTAB(INDEX1),2**24-1))/FLOAT(2**24-1) RETURN END c c------------------------------------ subroutine tri(tab,io,n) c implicit real*8(a-h,o-z) logical drap dimension tab(1),io(1) c 200 drap=.true. c do 201 i=1,n if(i+1.gt.n)goto 201 if(tab(i).gt.tab(i+1)) then c bidu=tab(i) tab(i)=tab(i+1) tab(i+1)=bidu c ibi=io(i) io(i)=io(i+1) io(i+1)=ibi c drap=.false. endif 201 continue if(.not.drap) then goto 200 endif c return end c----------------------------------------------------------- subroutine pregradient(ndat,numx,numy,igeo) c implicit real*8 (a-h,o-z) c parameter (mxs=200,mxd=50000) parameter (mxf=9) c c mxs=maxi de sources, mxf=maxi de fichiers de donnees, mxd=maxi de donnees c c dimension xdata(mxd),ydata(mxd),ddata(mxd) dimension weight(mxf),edata(mxd) dimension nv(mxf,2),irv(mxf) dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2) dimension a(mxd,mxs) dimension tapar(mxs,10),libre(mxs) dimension ixxx(mxs,2),iyyy(mxs,2) dimension xf(mxs),yf(mxs),ff(mxs),h(mxs) dimension w(mxs),r(mxs),tf(mxs) dimension us(mxs),ud(mxs),ut(mxs) c common/param/tapar,numsour,libre common/data1/xdata,ydata,ddata,edata,weight common/data2/nv,tsa,tsb,tsc,irv,dref common/data3/nfic common/defo/a COMMON/F1/XF,YF,FF,H,W,R,TF COMMON/F2/US,UD,UT COMMON/COMOB/Q,ST,CT,ELAS,U0 COMMON/COMOG/CC,E2,PI,IFU,AK,KH common/grid_mom/ixxx,iyyy c c c c ouverture fichier LOG c call lec_par (igeo,numx,numy) c c print*,'nombre de fichiers de donnees' read(*,*)nfic write(99,*)nfic,' fichiers de donnees' c call lec_data (igeo,numx,numy) c ndat =nv(nfic,2) c print*,'DIRECT' call direct c c sigma2=0.d0 do i=1,ndat sigma2=sigma2+(ddata(i))**2 enddo sigma2=dsqrt(sigma2) sigma2=1.0 print*,' COV= ',sigma2 c c cc...print write(99,*)'' print*,' ' print*,'DATA no fich , n1, n2, nn, reference rel, abs, poids' write(99,*)'DATA no , n1, n2, nn, reference, poids' do ii=1,nfic n1=nv(ii,1) n2=nv(ii,2) nn=n2-n1+1 if(irv(ii).ne.0) then nref=irv(ii) mref=nv(ii,1)+irv(ii)-1 else nref=0 mref=0 endif write(99,*)ii,n1,n2,nn,nref,mref,sngl(weight(ii)) print*,ii,n1,n2,nn,nref,mref,sngl(weight(ii)) enddo c return end c c---------------------------------------- subroutine lec_par (igeo,numx,numy) c------------------------------- c implicit real*8 (a-h,o-z) c parameter (mxs=200,mxd=50000) parameter (mxf=9) c character*100 fpar dimension tapar(mxs,10),libre(mxs) c common/param/tapar,numsour,libre COMMON/COMOG/CC,E2,PI,IFU,AK,KH c print*,'nom du fichier parametre de faille' read(*,'(a)')fpar write(99,*)fpar c print*,' 2 : geographique' print*,' 1 : UTM' read(*,*)igeo if(igeo.eq.2) then print*,' fuseau' read(*,*)ifu write(99,*)'geographique , fuseau ',ifu endif c print*,'nombre de failles along strike' read(*,*)numx print*,'nombre de failles along dip' read(*,*)numy write(99,*)'nb failles en X et Y ',numx,numy c print*,' voulez vous normaliser le slip 1/0' read(*,*)lnorm if(lnorm.eq.1) then print*,'valeur du slip (en cm)' read(*,*)slipnorm write(99,*)'slip normalise ', slipnorm else write(99,*)'slip non normalise' endif c c open (30,file=fpar) c read(30,*)numsour nnppaarr=0 do k=1,numsour read(30,*)nnuummss write(99,*)'' write(99,*)'source no ',nnuummss read(30,*)tapar(k,1),zz,tapar(k,2),zz read(30,*)tapar(k,3),zz,tapar(k,4),zz read(30,*)tapar(k,5),zz,tapar(k,6),zz,tapar(k,7),zz read(30,*)tapar(k,8),wli,tapar(k,9),zz,tapar(k,10),zz libre(k)=nint(wli) if(lnorm.eq.1.and.libre(k).ne.0) tapar(k,8)=slipnorm write(99,*)tapar(k,1),tapar(k,2) write(99,*)tapar(k,3),tapar(k,4) write(99,*)tapar(k,5),tapar(k,6),tapar(k,7) write(99,*)tapar(k,8),tapar(k,9),tapar(k,10) if(libre(k).ne.0) then write(99,*)' SLIP LIBRE',libre(k) nnppaarr=nnppaarr+1 else write(99,*)' SLIP FIXE',libre(k) endif c if(igeo.eq.2) then if(tapar(k,2).lt.0.0d0) then kh=-1 else kh=1 endif call geoutm(tapar(k,1),tapar(k,2),tapar(k,1),tapar(k,2)) endif call new2jcr(tapar(k,1),tapar(k,2),tapar(k,5), & tapar(k,6),tapar(k,3),tapar(k,4)) c write(*,*)(tapar(k,il),il=1,10) enddo close(30) write(99,*)'nb de parametres libre = ',nnppaarr write(*,*)'nb de parametres libre = ',nnppaarr if(nnppaarr.eq.0) then print*,' ARRET PROGRAMME' write(99,*)' ARRET PROGRAMME' stop endif c return end c c------------------------------------------------ c subroutine ecr_par(igeo,numx,numy) c implicit real*8 (a-h,o-z) c parameter (mxs=200,mxd=50000) parameter (mxf=9) c character*100 fpar dimension tapar(mxs,10),libre(mxs) c common/param/tapar,numsour,libre COMMON/COMOG/CC,E2,PI,IFU,AK,KH c print*,'nom du fichier sortie parametre de faille' read(*,'(a)')fpar write(99,*)fpar open (40,file=fpar) c nnuummss=0 write(40,*)numsour,' 5 0 0 ' do k=1,numsour nnuummss=nnuummss+1 write(40,*)nnuummss call jcr2new(tapar(k,1),tapar(k,2),tapar(k,5), & tapar(k,6),tapar(k,3),tapar(k,4)) if(igeo.eq.2) then call utmgeo(tapar(k,1),tapar(k,2),tapar(k,1),tapar(k,2)) endif write(99,*)'source no ',nnuummss write(40,*)tapar(k,1),zz,tapar(k,2),zz write(40,*)tapar(k,3),zz,tapar(k,4),zz write(40,*)tapar(k,5),zz,tapar(k,6),zz,tapar(k,7),zz write(40,*)tapar(k,8),zz,tapar(k,9),zz,tapar(k,10),zz enddo c c close(40) return end c c----------------------------------------------- subroutine lec_data(igeo,numx,numy) c implicit real*8 (a-h,o-z) c parameter (mxs=200,mxd=50000) parameter (mxf=9) c character*100 ficdata,fivec c dimension xdata(mxd),ydata(mxd),ddata(mxd) dimension weight(mxf),edata(mxd) dimension nv(mxf,2),irv(mxf) dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2) c common/data1/xdata,ydata,ddata,edata,weight common/data2/nv,tsa,tsb,tsc,irv,dref common/data3/nfic COMMON/COMOG/CC,E2,PI,IFU,AK,KH c amil=1000.d0 nt=1 do kk=1,nfic nv(kk,1)=nt print*,'' print*,'nom du fichier no ',kk write(99,*)'' write(99,*)'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' write(99,*)'nom du fichier no ',kk read(*,'(a)')ficdata write(99,*)ficdata c print*,'incertitude donnees (si lu dans fichier : 999)' read(*,*)erda write(99,*)'incertitude ',erda print*,'no du point de reference (=0 si absolu)' read(*,*)irv(kk) print*,'fichier de 4 vecteurs de projection' read(*,'(A)')fivec call planv(igeo,fivec,kk) print*,'poids sur les donnees' read(*,*) weight(kk) write(99,*)'poids ',weight(kk) c open(31,file=ficdata) if(erda.eq.999.d0) then 2 read(31,*,end=4) x,y,dam,erda if(igeo.eq.2) then call geoutm(x,y,x,y) endif xdata(nt)=x*amil ydata(nt)=y*amil ddata(nt)=dam edata(nt)=erda nt=nt+1 goto 2 else 3 read(31,*,end=4) x,y,dam if(igeo.eq.2) then call geoutm(x,y,x,y) endif xdata(nt)=x*amil ydata(nt)=y*amil ddata(nt)=dam edata(nt)=erda nt=nt+1 goto 3 endif 4 close(31) c nv(kk,2)=nt-1 print*,nv(kk,1),nv(kk,2) c c difference par rapport a une reference eventuelle c print*,'reference ?',irv(kk) c if(irv(kk).ne.0) then nref=nv(kk,1)+irv(kk)-1 dref(kk,1)=ddata(nref) print*,'point de reference' print*,'no, no, x , y, def ' write(99,*)'mesures relatives au point : no, x y def' write(99,*)irv(kk),sngl(xdata(nref)),sngl(ydata(nref)), *sngl(dref(kk,1)) call utmgeo(xdata(nref)/1000.,ydata(nref)/1000.,xdd,ydd) write(99,*)irv(kk),sngl(xdd),sngl(ydd), *sngl(dref(kk,1)) print*,irv(kk),nref,xdata(nref),ydata(nref),dref(kk,1) do jk=nv(kk,1),nv(kk,2) ddata(jk)=ddata(jk)-dref(kk,1) enddo else print*,'mesures absolues' write(99,*)'mesures absolues' dref(kk,1)=0.0 endif enddo c return end c c-------------------------------------------- subroutine planv(igeo,fivec,kk) c-------------------------------------------- c IMPLICIT REAL*8(A-H,O-Z) parameter (mxf=9) character*100 fivec,nul*1 dimension nv(mxf,2),irv(mxf) dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2) dimension px(10),py(10),tx(10),ty(10),tz(10) c common/data2/nv,tsa,tsb,tsc,irv,dref COMMON/COMOG/CC,E2,PI,IFU,AK,KH c amil=1000.d0 c open(35,file=fivec) write(99,*)fivec c read(35,'(a1)')nul read(35,*)nd do i=1,nd 1 read(35,*)px(i),py(i),tx(i),ty(i),tz(i) write(99,*)sngl(px(i)),sngl(py(i)),sngl(tx(i)), # sngl(ty(i)),sngl(tz(i)) if(igeo.eq.2) then call geoutm(px(i),py(i),px(i),py(i)) px(i)=px(i)*amil py(i)=py(i)*amil endif enddo close(35) c sax=0. sbx=0. scx=0. say=0. sby=0. scy=0. saz=0. sbz=0. scz=0. nc=0 do i=1,nd do j=i+1,nd do k=j+1,nd c x bx=((tx(i)-tx(k))-(((tx(j)-tx(k))*(px(i)-px(k))) &/(px(j)-px(k))))/(((py(i)-py(k))*(px(j)-px(k))-(py(j)-py(k))* &(px(i)-px(k)))/(px(j)-px(k))) ax=((tx(i)-tx(k))-bx*(py(i)-py(k)))/(px(i)-px(k)) cx=tx(i)-ax*px(i)-bx*py(i) sax=sax+ax sbx=sbx+bx scx=scx+cx c y by=((ty(i)-ty(k))-(((ty(j)-ty(k))*(px(i)-px(k))) &/(px(j)-px(k))))/(((py(i)-py(k))*(px(j)-px(k))-(py(j)-py(k))* &(px(i)-px(k)))/(px(j)-px(k))) ay=((ty(i)-ty(k))-by*(py(i)-py(k)))/(px(i)-px(k)) cy=ty(i)-ay*px(i)-by*py(i) say=say+ay sby=sby+by scy=scy+cy c z bz=((tz(i)-tz(k))-(((tz(j)-tz(k))*(px(i)-px(k))) &/(px(j)-px(k))))/(((py(i)-py(k))*(px(j)-px(k))-(py(j)-py(k))* &(px(i)-px(k)))/(px(j)-px(k))) az=((tz(i)-tz(k))-bz*(py(i)-py(k)))/(px(i)-px(k)) cz=tz(i)-az*px(i)-bz*py(i) saz=saz+az sbz=sbz+bz scz=scz+cz c print*,ax,bx,cx c print*,ay,by,cy c print*,az,bz,cz nc=nc+1 enddo enddo enddo tsa(kk,1)=sax/nc tsb(kk,1)=sbx/nc tsc(kk,1)=scx/nc tsa(kk,2)=say/nc tsb(kk,2)=sby/nc tsc(kk,2)=scy/nc tsa(kk,3)=saz/nc tsb(kk,3)=sbz/nc tsc(kk,3)=scz/nc write(99,*)'plan a b c ' do k=1,3 write(99,*)sngl(tsa(kk,k)),sngl(tsb(kk,k)),sngl(tsc(kk,k)) enddo return end c------------------------------------- subroutine direct c c IMPLICIT REAL*8(A-H,O-Z) c parameter (mxs=200,mxd=50000) parameter (mxf=9) c dimension xf(mxs),yf(mxs),ff(mxs),h(mxs) dimension w(mxs),r(mxs),tf(mxs) dimension us(mxs),ud(mxs),ut(mxs) dimension tapar(mxs,10) ,libre(mxs) dimension a(mxd,mxs) dimension xdata(mxd),ydata(mxd),ddata(mxd) dimension weight(mxf),edata(mxd) dimension nv(mxf,2),irv(mxf) dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2) c COMMON/F1/XF,YF,FF,H,W,R,TF COMMON/F2/US,UD,UT COMMON/COMOB/Q,ST,CT,ELAS,U0 COMMON/COMOG/CC,E2,PI,IFU,AK,KH common/param/tapar,numsour,libre common/defo/a common/data1/xdata,ydata,ddata,edata,weight common/data2/nv,tsa,tsb,tsc,irv,dref common/data3/nfic c trad=pi/180.d0 amil=1000.d0 smu = 3.3d10 c c do i=1,numsour xf(i)=tapar(i,1)*amil yf(i)=tapar(i,2)*amil ff(i)=tapar(i,3)*trad tf(i)=tapar(i,4)*trad h(i)=tapar(i,5)*amil w(i)=tapar(i,6)*amil r(i)=tapar(i,7)*amil slip=tapar(i,8) c usa=slip*dcosd(tapar(i,9)) uda=-slip*dsind(tapar(i,9)) c print*,usa,uda us(i)=usa ud(i)=uda ut(i)=0.0000001d0 c write(*,*)(tapar(i,k),k=1,10) c print*,'PARAMETRES',slip c print*,xf(i),yf(i),ff(i),tf(i),h(i),w(i),r(i),us(i),ud(i),ut(i) enddo call def_unit(numsour) return end c c------------------------------------- subroutine direct1 c c IMPLICIT REAL*8(A-H,O-Z) c parameter (mxs=200,mxd=50000) parameter (mxf=9) c dimension xf(mxs),yf(mxs),ff(mxs),h(mxs) dimension w(mxs),r(mxs),tf(mxs) dimension us(mxs),ud(mxs),ut(mxs) dimension tapar(mxs,10),libre(mxs) dimension a(mxd,mxs) dimension xdata(mxd),ydata(mxd),ddata(mxd) dimension weight(mxf),edata(mxd) dimension nv(mxf,2),irv(mxf) dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2) c COMMON/F1/XF,YF,FF,H,W,R,TF COMMON/F2/US,UD,UT COMMON/COMOB/Q,ST,CT,ELAS,U0 COMMON/COMOG/CC,E2,PI,IFU,AK,KH common/param/tapar,numsour,libre common/defo/a common/data1/xdata,ydata,ddata,edata,weight common/data2/nv,tsa,tsb,tsc,irv,dref common/data3/nfic c trad=pi/180.d0 amil=1000.d0 smu = 3.0d10 c c do i=1,numsour xf(i)=tapar(i,1)*amil yf(i)=tapar(i,2)*amil ff(i)=tapar(i,3)*trad tf(i)=tapar(i,4)*trad h(i)=tapar(i,5)*amil w(i)=tapar(i,6)*amil r(i)=tapar(i,7)*amil c slip=1.d20/(w(i)*r(i)*2*smu) usa=(tapar(i,8)) uda=(tapar(i,9)) us(i)=usa ud(i)=uda ut(i)=0.0000001d0 print*,'PAR',xf(i),yf(i),ff(i),tf(i),h(i),w(i),r(i) print*,usa,uda c write(*,*)(tapar(i,k),k=1,10) enddo call def_unit(numsour) return end C-------------------------------------------------------------------------- subroutine def_unit(numsour) C----------------------------------------- C Calcul de la deformation due a un ensemble de failles c IMPLICIT REAL*8(A-H,L,O-Z) c parameter (mxs=200,mxd=50000) parameter (mxf=9) dimension xf(mxs),yf(mxs),ff(mxs),h(mxs) dimension w(mxs),r(mxs),tf(mxs) dimension xdata(mxd),ydata(mxd),ddata(mxd) dimension weight(mxf),edata(mxd) dimension us(mxs),ud(mxs),ut(mxs) dimension nv(mxf,2),irv(mxf) dimension tsa(mxf,3),tsb(mxf,3),tsc(mxf,3),dref(mxf,2) dimension a(mxd,mxs) c COMMON/F1/XF,YF,FF,H,W,R,TF COMMON/F2/US,UD,UT common/defo/a common/data1/xdata,ydata,ddata,edata,weight common/data2/nv,tsa,tsb,tsc,irv,dref common/data3/nfic COMMON/COMOG/CC,E2,PI,IFU,AK,KH COMMON/COMOB/Q,ST,CT,ELAS,U0 c c cc...initialize c do jj=1,nfic dref(jj,2)=0.d0 enddo do kk = 1,nv(nfic,2) do i = 1,numsour a(kk,i) = 0. enddo enddo c c boucle sur les failles c DO 200 I=1,numsour FA=FF(I) TA=TF(I) HA=H(I) XPA=XF(I)-HA*DCOS(FA)/DTAN(TA) YPA=YF(I)+HA*DSIN(FA)/DTAN(TA) USA=US(I) UDA=UD(I) UTA=UT(I) WA = W(I) DA1=HA/DSIN(TA) DA2=DA1+WA RA=R(I) HD= DA2*DSIN(TA) C -------- boucle sur les points c do 201 jj=1,nfic do 202 kk=nv(jj,1),nv(jj,2) dxunit=0.d0 dyunit=0.d0 dzunit=0.d0 xx=xdata(kk) yy=ydata(kk) C Coulissage : IF(DABS(USA).LE.0.001) GOTO 51 UF=USA IM=2 CALL FAULT( XX,YY,XPA,YPA,FA,TA,RA,WA,HD,UF,IM,A1,B1,C1) dxunit = dxunit+A1 dyunit = dyunit+B1 dzunit = dzunit+C1 c Plongeante : 51 IF(DABS(UDA).LE.0.001) GOTO 52 IM=3 UF=-UDA CALL FAULT(XX,YY,XPA,YPA,FA,TA,RA,WA,HD,UF,IM,A1,B1,C1) dxunit = dxunit+A1 dyunit = dyunit+B1 dzunit = dzunit+C1 c Tensile : 52 IF(UTA.LE.0.001) GOTO 53 UF=UTA IM=1 CALL FAULT(XX,YY,XPA,YPA,FA,TA,RA,WA,HD,UF,IM,A1,B1,C1) dxunit = dxunit+A1 dyunit = dyunit+B1 dzunit = dzunit+C1 53 CONTINUE tsx=tsa(jj,1)*xx+tsb(jj,1)*yy+tsc(jj,1) tsy=tsa(jj,2)*xx+tsb(jj,2)*yy+tsc(jj,2) tsz=tsa(jj,3)*xx+tsb(jj,3)*yy+tsc(jj,3) a(kk,i)=dxunit*tsx+dyunit*tsy+dzunit*tsz c print*,a(kk,i),ddata(kk) c 202 continue c c difference a une reference eventuelle c if(irv(jj).ne.0) then c print*,'******** reference ********' nref=nv(jj,1)+irv(jj)-1 c print*,nref ref=a(nref,i) dref(jj,2)=ref+dref(jj,2) c print*,'ref sommeref',ref,dref(jj,2) do kk=nv(jj,1),nv(jj,2) a(kk,i)=a(kk,i)-ref enddo else dref(jj,2)=0.0 endif 201 continue c 200 CONTINUE write(99,*)'' write(99,*)'calcul direct' do jj=1,nfic write(99,*)'fichier no ',jj refd=(dref(jj,1)) refc=(dref(jj,2)) write(99,*)'no point,ref mesure, calcule, dif',irv(jj),sngl(refd), &sngl(refc),sngl(refc-refd) print*,'no point, ref mesure, calcule dif ',irv(jj),sngl(refd), &sngl(refc),sngl(refc-refd) enddo c RETURN END c cc---------------------------------------------------------------- C------------------------------------------------ SUBROUTINE FAULT(XX,YY,XP,YP,FI,TE,AL,W,D,U,IMOD,A,B,C) C----------------------------- C C Calcul du deplacement dans le repere geographique C J.C.Ruegg-1986 d'apres Okada(1985,BSSA) c c Faille en tension IMOD= 1 c Faille en coulissage IMOD= 2 c Faille plongeante IMOD= 3 C Coordonnes initiales XX,YY dans repere geographique UTM C Repere local de faille,Y1 c Repere Okada X2,Y2 c Sortie : Deplacements A,B,C surface libre c IMPLICIT REAL*8(A-H,O-Z) COMMON/COMOB/Q,ST,CT,ELAS,U0 c PI= 3.14159265359d0 SF=DSIN(FI) CF=DCOS(FI) CT=DCOS(TE) ST=DSIN(TE) C---------- Passage repere geographique au repere local X1,Y1 X1= (XX-XP)*SF + (YY-YP)*CF Y1=-(XX-XP)*CF + (YY-YP)*SF C---------- Passage au repere Okada X2= X1 Y2= Y1 + D*CT/ST C---------- Deformation dans le repere Okada QS1= X2+AL QS2= X2-AL P1= Y2*CT +D*ST P2= P1 - W Q= Y2*ST - D*CT U0 = U c C WRITE(*,*) U0,CT,ST,CF,SF c CALL DISLOK (QS1,P1,UX1,UY1,UZ1,IMOD) CALL DISLOK (QS1,P2,UX2,UY2,UZ2,IMOD) CALL DISLOK (QS2,P1,UX3,UY3,UZ3,IMOD) CALL DISLOK (QS2,P2,UX4,UY4,UZ4,IMOD) c A2 = UX1-UX2-UX3+UX4 B2 = UY1-UY2-UY3+UY4 C2 = UZ1-UZ2-UZ3+UZ4 c C--------- Deformation dans le repere geographique c A = A2*SF - B2*CF B = A2*CF + B2*SF C = C2 RETURN END C-------------------------------------------------- SUBROUTINE DISLOK(QSI,ETA,FUX,FUY,FUZ,IMOD) C-------------------------------------------------- C Calcul des deplacements ux,uy,uz C associes a une dislocation C FORMULES DE OKADA (1985, BSSA) C IMPLICIT REAL*8(A-H,O-Z) COMMON/COMOB/Q,ST,CT,ELAS,U0 c PI= 3.14159265359d0 PIP=PI+PI TT=ST/CT c BY=ETA*CT+Q*ST BD=ETA*ST-Q*CT R= DSQRT(QSI**2+ETA**2+Q**2) XX=DSQRT(QSI**2+Q**2) c RE= R+ETA RQ= R+QSI QRE= Q/R/RE QRQ= Q/R/RQ c ATQE=DATAN(QSI*ETA/Q/R) c AH= ETA*(XX+Q*CT)+XX*(R+XX)*ST BH= QSI*(R+XX)*CT TH= AH/BH c AJ1= -ELAS*QSI/(R+BD)/CT AJ2= -ELAS*DLOG(R+ETA) AJ3= ELAS*(BY/(R+BD)/CT-DLOG(R+ETA)) AJ4= ELAS*(DLOG(R+BD)-ST*DLOG(R+ETA))/CT AJ5=2.d0*ELAS*DATAN(TH)/CT c IF (IMOD.NE.1) GOTO 10 C----------Faille en tension : c FUX= Q*QRE - AJ3*ST*ST - AJ4*ST*ST*TT FUY=-BD*QRQ - ST*QRE*QSI + ST*ATQE - AJ1*ST*ST + AJ5*ST*ST*TT FUZ= BY*QRQ + CT*QSI*QRE - CT*ATQE - AJ5*ST*ST GOTO 20 c 10 IF (IMOD.NE.2) GOTO 11 C------------Faille de coulissage pur (strike-slip) : c FUX= -QSI*QRE - ATQE - AJ1*ST + AJ5*ST*TT FUY= -BY*QRE - Q*CT/RE - AJ2*ST + AJ3*ST + AJ4*ST*TT FUZ= -BD*QRE - Q*ST/RE - AJ4*ST GOTO 20 c 11 IF(IMOD.NE.3) GOTO 12 C----------Faille plongeante pure ( dip-slip) c FUX= -Q/R + AJ3*ST*CT + AJ4*ST*ST FUY= -BY*QRQ - CT*ATQE + AJ1*ST*CT - AJ5*ST*ST FUZ= -BD*QRQ - ST*ATQE + AJ5*ST*CT GOTO 20 c 20 FUX= FUX*U0/PIP FUY= FUY*U0/PIP FUZ= FUZ*U0/PIP c 12 RETURN END c c-------------------------------------------- C ------------------------------------------------ SUBROUTINE UTMGEO(XE,YE,xLON,yLAT) C -------------------------------------- C TRANSFORMATION UTM - GEOGRAPHIQUE C ENTREE : COORD. UTM EN KILOMETRES C SORTIE : COORD. GEOGRAPH. EN DEGRES C C HEMISPHERE NORD : KH=1 C SUD : KH=-1 C IMPLICIT REAL*8 (A-H,L,O-Z) COMMON/COMOG/CC,E2,PI,IFU,AK,KH c c print*,'utmgeo i ',xe,ye XX=XE*1000. YY=YE*1000. MI=IFU*6-183 AMI=FLOAT(MI)/0.9D0 RD=PI/200.D0 AM0=AMI*RD YR=YY if(KH.EQ.1) goto 11 YR=10.D6-YR 11 YR=YR/AK P=KH*YR*PI*0.5D-7 CO=DCOS(P) zL=P*KH S=E2*CO*CO R=CC/DSQRT(1.d0+S) PP=P*KH YR=YR-ARCME(zL) U=(XX-5.D5)/AK/R V=S*U*U Q=YR/R*(1.D0-V/2.D0) U=U*(1.D0-V/6.D0) P=zL+Q W=DEXP(U) AM=DATAN((W*W-1.D0)/(2.D0*W*DCOS(P))) V=DSIN(P)/DCOS(P) G=DATAN(V*COS(AM)) V=1.D0+S-1.5D0*DSIN(zL)*E2*(G-zL)*DCOS(zL) xLON=AM+AM0 yLAT=zL+V*(G-zL) xLON=(AM+AM0)/RD yLAT=yLAT/RD*KH xlon=xlon*0.9d0 ylat=ylat*0.9d0 c print*,'utmgeo o ',xlon,ylat c RETURN END C ------------------------------------------------ C ------------------------------------------------ FUNCTION ARCME(AL) C-------------------------- C Fonction Arc de m\002ridien en fonction de la latitude : IMPLICIT REAL*8(A-H,O-Z) COMMON/COMOG/CC,E2,PI,IFU,AK,KH AQ=E2*3.d0/4.d0 BQ=E2*E2*15.d0/16.d0 CQ=(E2**3)*35.d0/64.d0 DQ=(E2**4)*315.d0/512.d0 U=DSIN(AL)*DCOS(AL) V=DCOS(AL)**2 AJ2=AL+U AJ4=(AJ2*3.d0+U*V*2.d0)/4.d0 AJ6=(AJ4*5.d0+2.d0*U*V**2)/3.d0 AJ8=(AJ6*7.d0+4.d0*U*V**3)/8.d0 ARCME=CC*(AL-AQ*AJ2+BQ*AJ4-CQ*AJ6+DQ*AJ8) RETURN END c c----------------------------------------- C----------------------------------------------------------- SUBROUTINE GEOUTM(LON,LAT,XX,YY) C--------------------------------------------- C TRANSFORMATION GEOGRAPHIQUE - U.T.M. C ENTREE : LONGITUDE ET LATITUDE EN degres C SORTIE : COORD. UTM EN KILOMETRES C IMPLICIT REAL*8 (A-H,L,O-Z) COMMON/COMOG/CC,E2,PI,IFU,AK,KH c----------------------------------- c PI=4.d0*DATAN(1.d0) eps = 0.00000001 trg=200.d0/180.d0 c CC= 6399593.6257585d0 c EPRIM= 0.08209443794984d0 c E2=EPRIM**2 c IFU= 54 AK= 0.9996d0 c ---------------------------------- c print*,'geoutm i ',lon,lat lo=lon*trg la=lat*trg MI=IFU*6-183 AMI=DBLE(MI)/0.9D0 RD=PI/200.D0 AM=LO*RD AL=LA*RD AM0=AMI*RD CO=DCOS(AL) AMU=AM-AM0 XI=CO*DSIN(AMU) XI=0.5D0*DLOG((1.D0+XI)/(1.D0-XI)) c TA=DATAN(DSIN(AL)/(CO*DCOS(AMU)))-AL TA=DATAN2(DSIN(AL),CO*DCOS(AMU))-AL GN=CC/DSQRT(1.D0+E2*CO*CO) CX=E2*(CO*XI)**2 DX=GN*XI*(1.D0+CX/6.D0) DY=GN*TA*(1.D0+CX/2.D0) XX=500000.D0+AK*DX YY=AK*(DY+ARCME(AL)) if(LAT) 1,2,2 1 YY=(YY+10000000.d0) 2 CONTINUE xx=xx/1000. yy=yy/1000. c print*,'geoutm o ',xx,yy RETURN END c----------------------------------------- c-------------------------------------------- c subroutine de transformation de : c us,ud en slip, rake c subroutine usd2sr (slip,rake) c c us>0 senestre c us<0 dextre c ud>0 normal c ud<0 inverse c c -180 0 inverse c rake<0 normal c 90< |rake| <180 senestre c 0 < |rake| <90 dextre c IMPLICIT REAL*8(A-H,O-Z) c 1 continue us=slip ud=rake c slip=dsqrt(us**2+ud**2) rake=-datan2d(ud,us) c return end c c----------------------------------------- c c subroutine de transformation de : c slip,rake ---> en us ud c subroutine sr2usd (slip,rake) c---------------------------------------- c c us>0 senestre c us<0 dextre c ud>0 normal c ud<0 inverse c c -180 0 inverse c rake<0 normal c 90< |rake| <180 senestre c 0 < |rake| <90 dextre c IMPLICIT REAL*8(A-H,O-Z) c 1 continue c us=slip*dcosd(rake) ud=-slip*dsind(rake) c slip=us rake=ud c return end c c-------------------------------------------- c subroutine jcr2new (xs,ys,h,w,str,dip) c c transformation de xs ys h (prof) w (largeur) c en xq yq (en surface) h1 h2 (deux profondeurs de la faille) c IMPLICIT REAL*8(A-H,O-Z) c wh=h/dtand(dip) xq=xs-wh*dcosd(str) yq=ys+wh*dsind(str) h2=w*dsind(dip)+h c xs=xq ys=yq w=h2 c return end c c c----------------------------------------- c subroutine new2jcr (xq,yq,h1,h2,str,dip) c c transformation de xq yq (en surface) h1 h2 (deux profondeurs de la faille) c en xs ys h (prof) w (largeur) c IMPLICIT REAL*8(A-H,O-Z) c wh=h1/dtand(dip) xs=xq+wh*dcosd(str) ys=yq-wh*dsind(str) wa=(h2-h1)/dsind(dip) c xq=xs yq=ys h1=h1 h2=wa return end c c----------------------------------------