$debug
C-----------------------------------------------------------------------
C       T R E L I 2 .for
C-------------------------------------------------------------------------
C
C      GEODESIE - Passage coordonnÜes UTM vers  wgs84
C                 
C                 connaissant les paramÊtres de la transformation
C-------------------------------------------------------------------------
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*12 FICH2,FICH3,FICH4,FICH5,FICH6
      CHARACTER*4 CHSTA
c
C
      COMMON/COELI/CC,EP2,PI,RAD,IFU
      COMMON/COMOI/MPAR,NDAT,NK
      DIMENSION XU(50),YU(50),ZU(50),NSTU(50)
      DIMENSION P(7),PO(7)
      DIMENSION TEX(20)
      DIMENSION CHSTA(50)
C     
c -- initialisation des parametres :
c
      MPAR=7
      TX=100.
      TY= 50.
      TZ=200.
      DELTA=50.D-6
      RX=1.d-8
      RY=1.d-8
      RZ=1.d-8
c
   91  PI=4.d0*DATAN(1.d0)
c ..............3.141592653589793D0
       RAD=PI/200.d0
c
C----------
      write(*,*) ' Programme TRELI 2 .for :'
      write(*,*) ' -----------------------'
c
      WRITE(*,*) ' Nom du fichier EDITION :'
      READ(*,'(A)') FICH4
      OPEN(4,FILE=FICH4,status='NEW',err=91)
      write(4,*) ' Programme TRELI 2.for :'
      write(4,*) ' ---- utm-->-wgs84 -----'
c
   11 WRITE(*,*) ' Nom fichier projection UTM : '
      READ(*,'(A)') FICH2  
      OPEN(2,FILE=FICH2,STATUS='OLD',err=11)
      READ(2,99) (TEX(I),I=1,20)
   99 FORMAT(20A4)
      WRITE(*,99) (TEX(I),I=1,20)
      READ(2,*) CC1,EPR1
      READ(2,*) IFU
      EP1C=EPR1*EPR1
      write(*,*) CC1, EPR1, EP1C
      write(4,'(A)')  FICH2
      WRITE(4,99) (TEX(I),I=1,20)
      write(4,*) CC1, EPR1, EP1C
       AA1C=CC1*CC1/(1.d0+EP1C)
       AA1=dsqrt(AA1C)
       BB1=AA1C/CC1
       FF1=(AA1-BB1)/AA1     
c      write(4,*) '     - demi-axes et applatissement :'
c      write(4,*) AA1,BB1,FF1
      write(4,*) '-----------------------------------------------------'
      close(2)
c 
c -- paramÊtres wgs84 :
c
      CC2=6399593.6257585d0
      EPR2=0.08209443794984d0
      EP2C=EPR2*EPR2
c      WRITE(4,*) ' ParamÊtres projection WGS84 : '
c      write(4,*) CC2, EPR2, EP2C
       AA2C=CC2*CC2/(1.d0+EP2C)
       AA2=dsqrt(AA2C)
       BB2=AA2C/CC2
       FF2=(AA2-BB2)/AA2     
c      write(4,*) '     - demi-axes et applatissement :'
c      write(4,*) AA2,BB2,FF2
c      write(4,*) '-----------------------------------------------------'
C
      CALL GETDAT(IAN,MOIS,JOUR)
c      write(4,114) IAN,MOIS,JOUR
C
c -- coordonnees UTM a transformer :
c
   13 WRITE(*,*) ' Nom du fichier des coordonnees UTM :'
c      WRITE(4,*) '---------------- Coordonnees UTM -------- :'
      READ(*,102) FICH6
      OPEN(6,FILE=FICH6,STATUS='OLD',err=13) 
c
      CC=CC1
      EP2=EP1C
c
c      write(4,*) '------ Coordonnees UTM -------- : '
      READ(6,99) (TEX(I),I=1,20)
      WRITE(4,99) (TEX(I),I=1,20)
c ---
c   10 write(4,*) ' Nom fichier paramÊtres transformation:'
c      READ(*,'(A)') FICH3  
c      OPEN(3,FILE=FICH3,STATUS='OLD',err=10)
c      READ(3,99) (TEX(I),I=1,20)
c      WRITE(*,99) (TEX(I),I=1,20)
c      READ(3,*) TX,TY,TZ
c      READ(3,*) DELTA
c      READ(3,*) RX,RY,RZ
c      write(4,*) ' -------Transformation UTM - wgs84 -----------------'
c      write(4,'(A)')  FICH3
c      WRITE(4,99) (TEX(I),I=1,20)
c ---
      READ(6,*) TX,TY,TZ,DELTA,RX,RY,RZ
c
      write(4,150) TX,TY,TZ
      write(4,151) DELTA
      write(4,152) RX,RY,RZ      
      write(4,*) '----------------------------------------------------'
      close(3)
c
c -- entree des coordonnees UTM :
c
      J=0
  131 J=J+1
      READ(6,154,end=132) IC,NSTU(J),XU(J),YU(J),ZU(J),ER,CHSTA(J)
      write(4,133) J,NSTU(J),XU(J),YU(J),ZU(J)
c
      goto 131
  132 KMU=J-2
  133 format(I3,2x,I4,2x,3(F12.3,1x))
c
c -- Initialisation des paramŠtres
c
      P(1)=TX
      P(2)=TY
      P(3)=TZ
      P(4)=DELTA
      P(5)=RX
      P(6)=RY
      P(7)=RZ
c
 102  format(A12)
 114  format('  Date : ',I4,2(1X,I2))
 115  format('  Fichier Ellipsoide UTM : ',A12,3X,
     *' Fichier Ellipsoide wgs84 : ',A12)
 116  format('  Fichier coord. wgs84 (data) : ',A12)
 117  format('  Fichier coordonnees U T M   : ',A12)
c
c --- calcul coordonnees wgs de tous les pts UTM :
      write(4,*) ' Coordonnees WGS84 et UTM:'
      write(4,*) ' --------------------------'
c
      IC=1
      ER=1.
      DO 40 J=1,KMU
c
      CALL UTMGEO(XU(J),YU(J),AL,AF)
      CALL GEOCENTR(AL,AF,ZU(J),X,Y,Z)
c
      XT=X+P(1)+P(4)*X+P(5)*Z-P(6)*Y
      YT=Y+P(2)+P(4)*Y+P(6)*X-P(7)*Z
      ZT=Z+P(3)+P(4)*Z+P(7)*Y-P(5)*X
c
      write(4,155) IC,NSTU(J),XT,YT,ZT,ER,CHSTA(J)
   40 continue
c
      write(4,*) '0    0   0   0   0   0   0   0'
c
  150 format(10x,' TX =',F12.4,3x,' TY =',F12.4,3x,' TZ =',F12.4)
  151 format(10x,' Delta=',F12.8,3x,'sig.:',F12.9)
  152 format(10x,' RX =',F12.9,3x,' RY =',F12.9,3x,' RZ =',F12.9)
  154 FORMAT (I1,1X,I4,1X,2(F11.3,1X),F8.3,1X,F8.3,A4)
  155 FORMAT (I1,1X,I4,1X,3(F12.3,1X),F8.3,1X,A4)
  156 format(I4,1x,3(F11.3,1x),2x,2(F10.2,1x),1x,F7.2)
c
   90 write(*,*) ' -- FIN --'
      write(4,*)'-----------------------------------------------------'
      write(4,114) IAN,MOIS,JOUR
      write(4,117) FICH6
      write(4,*)'-----------------------------------------------------'
      close (4)
      CLOSE(6)
      CLOSE(7,err=93)
  93  STOP
      END
C -----------------------------------------------------------
       SUBROUTINE GEOCENTR(LON,LAT,H,X,Y,Z)
C-----------------------
C---Passage coordonnees geographiques --> geocentriques
c         unitÜs : grades et mÊtres
c
       IMPLICIT REAL*8(A-H,L,O-Z)
       COMMON/COELI/CC,EP2,PI,RAD,IFU
c
       LATR=LAT*PI/200.D0
       LONR=LON*PI/200.D0
       CB=DCOS(LATR)
       SB=DSIN(LATR)
       CL=DCOS(LONR)
       SL=DSIN(LONR)
C---Grande normale GN
c   Note: CC rayon de courbure au pole - EP2 = EPRIM **2
c         E = seconde excentricite - e2 =(a2-b2)/b2
       GN=CC/DSQRT(1.D0+EP2*CB*CB)
c
c---Coordonnees tridimensionnelles
c
       X=(GN+H)*CB*CL
       Y=(GN+H)*CB*SL
       Z=(GN/(1.d0+EP2)+H)*SB
c
       RETURN
       END
C-------------------------------------------------------------
      SUBROUTINE UTMGEO(XX,YY,LON,LAT)
C --------------------------------------
C         TRANSFORMATION  UTM - GEOGRAPHIQUE
C           ENTREE :  COORD. UTM  EN  METRES
C           SORTIE :  COORD. GEOGRAPH. EN GRADES
C
C           HEMISPHERE NORD : KH=1
C                      SUD  : KH=-1
C
      IMPLICIT REAL*8 (A-H,L,O-Z)
      COMMON/COELI/CC,EP2,PI,RAD,IFU
c
      E2=EP2      
      KH=1
      AK=0.9996d0
      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
c
      RETURN
      END
C ------------------------------------------------
       FUNCTION ARCME(AL)
C--------------------------
C         Fonction Arc de m‚ridien en fonction de la latitude :

      IMPLICIT REAL*8(A-H,L,O-Z)
      COMMON/COELI/CC,EP2,PI,RAD,IFU
c
      E2=EP2
      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--------------------------FIN-----------------------------------