$large
$debug
C-----------------------------------------------------------------------
C        REPET2.FOR
C-------------------------------------------------------------------------
C
C      Calcul de la repetabilite de mesures GPS
C
C	- lecture d'un fichier "coordonnees" issu d'un ajustement
c
c	- lecture d'un fichier donnees au format AG3D
c
c  	- classement par ordre (du fichier.cor)
c	- calcul du residu O-C pour chaque composante dX,dY,dZ
c                                                     dE,dN,dH et distance
C-------------------------------------------------------------------------
C        J.C. Ruegg -  Juin 93
c--------------------------------------------------------------------------
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*12 FICH0,FICH1,FICH2,FICH3,FICH4
      CHARACTER*8 FIC,FFIC
      CHARACTER*4 CHSTA, CHVIS, CH1,CH2
c
      COMMON/COMOB/CC,EP2,PI,RAD,IFU
c
c -- coordonnees des pts
c
      DIMENSION XC(50),YC(50),ZC(50),NST(50)
c
c -- dimensionnement donnees lues
c
      DIMENSION LSTA(400),LVIS(400)
      DIMENSION CHSTA(400),CHVIS(400)
      DIMENSION VX(400),VY(400),VZ(400),EX(400),EY(400),EZ(400)
c                           
c -- dimensionnement des ecarts o-c --> classement 1 par ordre des stations
c
      DIMENSION EEX(400),EEY(400),EEZ(400),EEE(400),EEN(400),EEH(400)
      DIMENSION EED(400),DDK(400)
c                           
c -- dimensionnement des ecarts o-c --> classement 2 par
c                           ordre des ecarts en longueur
c
      DIMENSION FEX(400),FEY(400),FEZ(400),FEE(400),FEN(400),FEH(400)
      DIMENSION FED(400),FDK(400)
c
      DIMENSION JOUR(400),IDAT(400),JDAT(400),FFIC(400)
      DIMENSION TEX(20),TEX1(20)
c
       PI=4.d0*DATAN(1.d0)
c ..............3.141592653589793D0
       RAD=PI/200.d0
       AK=6.371d6*RAD
c
C----------
      write(*,*) ' Programme REPET2.for'
      write(*,*) ' --------------------'
c
   1  WRITE(*,*) ' Nom du fichier EDITION :'
      READ(*,'(A)') FICH4
      OPEN(4,FILE=FICH4,status='NEW',err=1)
      write(4,*) ' Programme REPET2.for :'
      write(4,*) ' ---------------------'
c
      CALL GETDAT(IAN,MOIS,JOU)
c      write(4,114) IAN,MOIS,JOU
c
c------- Fichier projection = paramŠtres de l'ellipsoide CC, Eprim
c           cc=rayon courbure au pole - Eprim=seconde excentricite
c
   2  WRITE(*,*) ' Nom fichier projection : PRJ : '
      READ(*,'(A)') FICH0  
      OPEN(3,FILE=FICH0,STATUS='OLD',err=2)
      READ(3,99) (TEX(I),I=1,20)
   99 FORMAT(20A4)
      WRITE(*,99) (TEX(I),I=1,20)
      READ(3,*) CC,EPRIM
      READ(3,*) IFU
      EP2=EPRIM*EPRIM
      write(*,*) CC, EPRIM, EP2
      write(4,*) ' ---------- Ellipsoide : ---------------------------'
      write(4,'(A)')  FICH0
      WRITE(4,99) (TEX(I),I=1,20)
      write(4,*) CC, EPRIM, EP2
       AA2=CC*CC/(1.d0+EP2 )
       AA=dsqrt(AA2)
       BB=AA2/CC
       FF=(AA-BB)/AA     
      write(4,*) '     - demi-axes et applatissement :'
      write(4,*) AA,BB,FF
      write(4,*) '-----------------------------------------------------'
      close(3)
c
c---------- Entree des coordonnees de stations :
c
   3  WRITE(*,*) ' Nom du fichier stations : STA :'
      WRITE(4,*) '---------- CoordonnÜes a priori ------------------- :'
      READ(*,'(A)') FICH1
      OPEN(5,FILE=FICH1,STATUS='OLD',err=3) 
c
      READ(5,99) (TEX(I),I=1,20)
      WRITE(*,99)(TEX(I),I=1,20)
      WRITE(4,99)(TEX(I),I=1,20)
c
   4  WRITE(*,*) ' ----------------------------------'
      write(*,*) ' Type de coordonnÜes a priori :'
      write(*,*) '                - geocentriques:  0'             
      write(*,*) '                - grades     :    1'             
      write(*,*) '                - degres dec.:    2'
      write(*,*) '                - U.T.M.     :    3'
      write(*,*) '                - coord. locales: 4'
      write(*,*) ' ----------------------------------'
c
      READ(*,102,err=4) ISYS
      FAC=1.d0
c
      WRITE(*,105)
  105 FORMAT(/,5X,'COORDONNEES INITIALES :')
      WRITE(*,106)
  106 FORMAT(7X,'N0',6X,'XI',15X,'YI',16X,'Z',10X,'ERP',4X,'ERPZ',/)
      write(4,*) ' '
c
       K=0
   10  K=K+1
       if(ISYS.NE.0) goto 11
       READ(5,*,end=19) ICOD,NSTA,XX,YY,ZZ
       WRITE(*,100) ICOD,NSTA,XX,YY,ZZ
       WRITE(4,100) ICOD,NSTA,XX,YY,ZZ
  100 FORMAT (I1,1X,I4,1X,3(F12.3,3X),2x,3(1x,F8.4))
  103 FORMAT (12x,2(F13.8,2x),F8.3)
      goto 15
c
   11 if(ISYS.EQ.2) FAC=0.9d0
      if(ISYS.GT.2) goto 12 
       READ(5,*,end=19) ICOD,NSTA,XLON,YLAT,HE
       XLON=XLON/FAC
       YLAT=YLAT/FAC
      CALL GEOCENTR(XLON,YLAT,HE,XX,YY,ZZ)
       WRITE(*,101) ICOD,NSTA,XLON,YLAT,HE
       WRITE(4,101) ICOD,NSTA,XLON,YLAT,HE
       WRITE(4,100) ICOD,NSTA,XX,YY,ZZ
  101  FORMAT (I1,1X,I4,5X,2(F13.9),F8.3,1X,3(1x,F8.4))
      goto 15
c
   12  READ(5,*,end=19) ICOD,NSTA,XUTM,YUTM,HE
       CALL UTMGEO(XUTM,YUTM,XLON,YLAT)
       CALL GEOCENTR(XLON,YLAT,HE,XX,YY,ZZ)
       WRITE(*,102) ICOD,NSTA,XUTM,YUTM,HE
       WRITE(4,102) ICOD,NSTA,XUTM,YUTM,HE
  102 FORMAT (I1,1X,I4,1X,2(F11.3,1X),F8.3,1X,3(1x,F8.4))
       WRITE(4,100) ICOD,NSTA,XX,YY,ZZ
       goto 15
c
   13  READ(5,*,end=19) ICOD,NSTA,XL,YL,ZL
       if(ICOD.NE.9) goto 14
       CALL CENTRGEO(XL,YL,ZL,XLON0,YLAT0,HEL0)
       NSTA0=NSTA
      write(4,*) ' Les coordonnees entrees sont des coordonnees locales'
      write(4,*) ' Coordonnees geocentriques et ellipsoidiques :'
      write(4,*) '                                 du point principal :'
      write(4,143) NSTA0, XL, YL, ZL
      write(4,144) XLON0,YLAT0,HEL0
  143 format(5x,I4,2X,' Xg = ',F13.4,'  Yg = ',F13.4,'  Zg = ',F13.4) 
  144 format(11x,' Xe = ',F13.8,' Ye = ',F13.8,' He = ', F8.3)
       goto 13
c
   14 if(NSTA.NE.NSTA0) goto 141
      XX0=XX
      YY0=YY
      ZZ0=ZZ
      XLON=XLON0
      YLAT=YLAT0
      HE=HEL0
       CFI=DCOS(YLAT*RAD)
       RN0=CC/dsqrt(1.d0+EP2*CFI*CFI)
       RO=RN0/(1.d0+EP2*CFI*CFI)
       FAL=RN0/(RN0+HEL0)
       WRITE(*,101) ICOD,NSTA,XLON,YLAT,HE
       WRITE(*,104) XX,YY,HE
       WRITE(4,101) ICOD,NSTA,XLON,YLAT,HE
       WRITE(4,104) XX,YY,HE
  104  FORMAT (11X,2(F13.4),F8.3)
      goto 15
c
  141 DLAM=(XX-XX0)*FAL/RN0/CFI
      DFI =(YY-YY0)*FAL/RO
      XLON=XLON0+DLAM/RAD
      YLAT=YLAT0+DFI/RAD
      HE=HEL0+ZZ-ZZ0
c
      CALL GEOUTM(XLON,YLAT,XUTM,YUTM)
c
       WRITE(*,101) ICOD,NSTA,XLON,YLAT,HE
       WRITE(*,104) XX,YY,HE
       WRITE(4,101) ICOD,NSTA,XLON,YLAT,HE
       WRITE(4,104) XUTM,YUTM
       WRITE(4,104) XX,YY,HE
c
   15 continue
c
      XC(K)=XX
      YC(K)=YY
      ZC(K)=ZZ      
      NST(K)=NSTA
c
      GOTO 10
c
   19 MSTA=K-1
c
C---------- Entree des donnees :
c
    6 WRITE(*,*) 'Nom du fichier donn‚es : DAT :'
      READ(*,'(A)') FICH2
      OPEN (6,FILE=FICH2,STATUS='OLD',err=6)
c
      write(4,*) ' '
      WRITE(4,*) '----------- Donnees ----------------------- :'
      write(4,*) ' '
c
      READ(6,99) (TEX1(I),I=1,20)
      WRITE(*,99)(TEX1(I),I=1,20)
      WRITE(4,99)(TEX1(I),I=1,20)
c
      JJ=0
   60 continue
c
      READ(6,108,end=70) IDON,NSTA,NVIS,DATI,ERD,IBID1,IBID2,CH1,CH2,JR
     *,FIC
      write(*,108) IDON,NSTA,NVIS,DATI,ERD,IBID1,IBID2,CH1,CH2,JR
  107 format(I5,2I6,F12.3,F8.3,2F6.0,A4,1x,A4,2x,I6)
  108 format(I2,2x,I3,2x,I3,1x,F12.4,1x,F8.4,2x,I2,2x,I2,2x,A4,2x,A4,I6,
     *7x,A8)
c
      IF(IDON.EQ.12) goto 62
      IF(IDON.EQ.13) goto 63
      IF(IDON.EQ.14) goto 64
      goto 60
c
   62 JS=NSTA
      JV=NVIS
      JJ=JJ+1
      LSTA(JJ)=NSTA
      LVIS(JJ)=NVIS
      VX(JJ)=DATI
      EX(JJ)=ERD
      CHSTA(JJ)=CH1
      CHVIS(JJ)=CH2
      JOUR(JJ)=JR
      FFIC(JJ)=FIC
      goto 60
c
   63 if(NSTA.NE.JS) goto 60
      if(NVIS.NE.JV) goto 60
      VY(JJ)=DATI
      EY(JJ)=ERD
      goto 60
c
   64 if(NSTA.NE.JS) goto 60
      if(NVIS.NE.JV) goto 60
      VZ(JJ)=DATI
      EZ(JJ)=ERD
      goto 60
c
   70 NDAT=JJ     
      write(*,*) ' Entree des donnees terminees !'
c
c -- sortie des donnees lues :
c
      DO 20 J=1,NDAT
       DD=DSQRT(VX(J)**2+VY(J)**2+VZ(J)**2)
      write(4,201) J,LSTA(J),LVIS(J),VX(J),VY(J),VZ(J),DD,EX(J),EY(J),
     *EZ(J)
   20 continue
c
      write(4,*) '-----------------------------------------------------'
  201 format(I3,2I4,3(1x,F12.4),1x,F10.2,1x,3(F6.4))
c
c -- sortie des donnees
c
      WRITE(*,*) ' ------------------------------------------'
      write(*,*) ' Sortie : ecarts sur dX, dY, dZ, dD : 1'
      write(*,*) '          ecarts sur dE, dN, dH, dD : 2'
      write(*,*) '                '             
      write(*,*) ' ------------------------------------------'
c
      READ(*,102) ISOR
c
c -- classement et calcul des ecarts O-C
c          en fonction de l'ordre des stations
c
      write(4,*) ' ----------------------------------------------------'
      write(4,*) ' Ecart / classement des stations :'
      if(ISOR.EQ.2) goto 21
      write(4,*) '     Dist.(km), dX, dY, dZ, dD, NSTA, NVIS --' 
      goto 211
   21 write(4,*) '     Dist.(km)  dE, dN, dH, dD, NSTA, NVIS --'
  211 write(4,*) ' ------'
c
       JK=0
c
c -- notations : JK : indice du classement en fonction de l'ordre des stations
c                K1,K2,MSTA : ordre des stations et nb de stations
c                J, NDAT : ordre des donnees lues et nb de donnees lues
c
      DO 30 K1=1,MSTA	
       NSTA=NST(K1)
       KK=K1+1
      DO 30 K2=KK,MSTA
       NVIS=NST(K2)
c
      DO 32 J=1,NDAT
       if(LSTA(J).NE.NSTA) goto 33
       if(LVIS(J).NE.NVIS) goto 33
c
c -- calcul de dx,dy,dz, dE,dN,dH d'apres les coordonnees
c
      XAC=XC(K1)
      YAC=YC(K1)
      ZAC=ZC(K1) 
      XBC=XC(K2)
      YBC=YC(K2)
      ZBC=ZC(K2)
      DXC=XBC-XAC
      DYC=YBC-YAC
      DZC=ZBC-ZAC
      CALL TRLOC (XBC,YBC,ZBC,XAC,YAC,ZAC,DEC,DNC,DHC)
      DISTC=DSQRT(DXC**2+DYC**2+DZC**2)
       goto 35
   33 if(LSTA(J).NE.NVIS) goto 32
      if(LVIS(J).NE.NSTA) goto 32
      XAC=XC(K2)
      YAC=YC(K2)
      ZAC=ZC(K2) 
      XBC=XC(K1)
      YBC=YC(K1)
      ZBC=ZC(K1)
      DXC=XBC-XAC
      DYC=YBC-YAC
      DZC=ZBC-ZAC
      CALL TRLOC (XBC,YBC,ZBC,XAC,YAC,ZAC,DEC,DNC,DHC)
      DISTC=DSQRT(DXC**2+DYC**2+DZC**2)
c
c -- calcul de dE,DN,DH observe, correspondant a la donnee dX,dY,dZ
c
   35 JK=JK+1
      XA=XAC
      YA=YAC
      ZA=ZAC
      XB=XA+VX(J)
      YB=YA+VY(J)
      ZB=ZA+VZ(J)
      CALL TRLOC (XB,YB,ZB,XA,YA,ZA,DEO,DNO,DHO)
      DISTO=DSQRT(VX(J)**2+VY(J)**2+VZ(J)**2)
c
c -- calcul des ecarts O-C
c
      ECX=VX(J)-DXC
      ECY=VY(J)-DYC
      ECZ=VZ(J)-DZC
      ECE=DEO-DEC
      ECN=DNO-DNC
      ECH=DHO-DHC
      ECDD=DISTO-DISTC
      DISTK=DISTC*.001d0
c
c
c	write(4,205) J,VX(J),DXC,ECX      
c	write(4,205) J,VY(J),DYC,ECY      
c	write(4,205) J,VZ(J),DZC,ECZ      
c	write(4,205) J,DISTO,DISTC,ECDD
  205 format('     - ',I4,3(1x,F12.3))      
c
c -- mise en memoire des O-C
c
      EEX(JK)=ECX
      EEY(JK)=ECY
      EEZ(JK)=ECZ
      EEE(JK)=ECE
      EEN(JK)=ECN
      EEH(JK)=ECH
      EED(JK)=ECDD
      DDK(JK)=DISTK
      IDAT(JK)=J
c
c -- sortie
c
      if (ISOR.EQ.2) goto 36
      write(4,202) DISTK,ECX,ECY,ECZ,ECDD,NSTA,NVIS
      goto 32
   36 write(4,202) DISTK,ECE,ECN,ECH,ECDD,NSTA,NVIS
c
  202 format(F8.3,3(1x,F7.4),2X,F7.4,3x,2I5)
c

   32 continue
   30 continue
      NJK=JK
c
      write(4,206) MSTA,NDAT,NJK
  206 format('   MSTA=',I4,'   NDAT='I4,'   NJK=',I4)
c
c --classement par ordre croissant de l'erreur sur D
c           NJK= nb de donnees classees
c           JL = indice du nouveau classement
c
      write(4,*) ' ----------------------------------------------------'
      write(4,*) ' Ecart / classement par ordre croissant erreur sur D:'
      if(ISOR.EQ.2) goto 31
      write(4,*) '     Dist.(km), dX, dY, dZ, dD, NSTA, NVIS --' 
      goto 311
   31 write(4,*) '     Dist.(km)  dE, dN, dH, dD, NSTA, NVIS --'
  311 write(4,*) ' ------'
c
      JL=0
c
c  recherche du maximum de EED
c    
   41 JL=JL+1
       VMAX=DABS(EED(1))
      DO 40 K=1,NJK
       VAL=DABS(EED(K))
       if(VAL.LE.VMAX) goto 40
       VMAX=VAL
       KVAL=K
   40 continue
c
c      write(4,401) JL,KVAL,VMAX,IDAT(KVAL)
  401 format(I5,I5,F10.3,I5)
c
      FEX(JL)=EEX(KVAL)
      FEY(JL)=EEY(KVAL)
      FEZ(JL)=EEZ(KVAL)
      FEE(JL)=EEE(KVAL)
      FEN(JL)=EEN(KVAL)
      FEH(JL)=EEH(KVAL)
      FED(JL)=EED(KVAL)
      FDK(JL)=DDK(KVAL)
      JDAT(JL)=IDAT(KVAL)
       EED(KVAL)=0.d0
       KVAL=1
c
      if (JL.LT.NJK) goto 41
      NJL=JL
c
      DO 42 L=1,NJL 
c
c -- recherche des sigles et no de la donnee originale
c
      J=JDAT(L)
      NSTA=LSTA(J)
      NVIS=LVIS(J)
      CH1=CHSTA(J)
      CH2=CHVIS(J)
      JR=JOUR(J)
      FIC=FFIC(J)
c
      ECX=FEX(L)
      ECY=FEY(L)
      ECZ=FEZ(L)
      ECE=FEE(L)
      ECN=FEN(L)
      ECH=FEH(L)
      ECD=FED(L)
      ECX=FEX(L)
      DISTK=FDK(L)
c
      if (ISOR.EQ.2) goto 421
      write(4,204) DISTK,ECX,ECY,ECZ,ECD,NSTA,NVIS,CH1,CH2,JR,FIC
      goto 422
  421 write(4,204) DISTK,ECE,ECN,ECH,ECD,NSTA,NVIS,CH1,CH2,JR,FIC
  422 continue
   42 continue
c
      write(4,*)'-----------------------------------------------------'
      write(4,114) IAN,MOIS,JOU
      WRITE(4,112) FICH0
      write(4,115) FICH1,FICH2
      write(4,*)'-----------------------------------------------------'
 112  format('  Nom fichier projection : PRJ : ', A12)
 114  format('  Date : ',I4,2(1X,I2))
 115  format('  Fichier entr‚e  STA : ',A12,3X,
     *' Fichier donn‚es DAT : ',A12)
 116  format('  Fichier sortie  COR : ',A12)
 204  format(F8.3,1x,4(F9.4,1x),1x,2I5,2x,A4,1x,A4,1x,I5,1x,A8)
      close (4)
      CLOSE(5)
      CLOSE(6)
      STOP
      END
c --------------------------------------------------
      SUBROUTINE TRLOC (X,Y,Z,X0,Y0,Z0,U,V,W)
c --------------------------------------------------
c	Transformation de coordonn‚es geocentriques en coordonnees locales
c	sur l'ellipsoide : origine au pt  X0,Y0,Z0
c	entree : X,Y,Z - sortie:  U,V,W
c       angles de rotation en grades: AL = Lambda0, longitude de Mo (xo,yo,zo)
c   	                 FI = Phi0 = latitude du pt Mo (xo,yo,zo)
c	Rotation de (AL+100gr) autour de Oz, et (100gr-FI) autour de Ox
c
      IMPLICIT REAL*8(A-H,O-Z)
      PI=4.d0*DATAN(1.d0)
       RAD=PI/200.d0
c
c -- calcul des composantes des rotations:
c
       CALL CENTRGEO(X0,Y0,Z0,AL0,FI0,H0)
c
       CF0=dcos(FI0*RAD)
       SF0=dsin(FI0*RAD)
       CL0=dcos(AL0*RAD)
       SL0=dsin(AL0*RAD)
c
c ------  translation :
       X1=X-X0
       Y1=Y-Y0
       Z1=Z-Z0
c -------- matrice rotation z(lambda+100), x(100-phi):
       r11=-SL0
       r12=CL0
       r13= 0.d0
       r21=-CL0*SF0
       r22=-SL0*SF0
       r23= CF0
       r31= CL0*CF0
       r32= SL0*CF0
       r33= SF0
c ------- nouvelles coordonnees:
c
       U = X1*r11+Y1*r12+Z1*r13
       V = X1*r21+Y1*r22+Z1*r23
       W = X1*r31+Y1*r32+Z1*r33
c
       RETURN
       END
c---------------------------------------------------------------------
c=============================================
      SUBROUTINE CNORM(X,Y,Z,DN1,DN2,DN3)    
C---------------------------------------------
c----Calcul des composantes du vecteur normal Û l'ellipsoide
c        au pt Po 
c
       IMPLICIT REAL*8 (A-H,L,O-Z)
      COMMON/COMOB/CC,EP2,PI,RAD,IFU
C
       A=CC/DSQRT(1.D0+EP2)
       B=A*A/CC
       P= DSQRT(X**2+Y**2)
       PB=B*P
       ZA=Z*A
       TETA= DATAN2(ZA,PB)
       LBD= DATAN2(Y,X)
       ST=DSIN(TETA)
       CT=DCOS(TETA)
C
       FX=Z+EP2*B*ST*ST*ST
       FY=P-EP2*A*CT*CT*CT/(1.d0+EP2)
       FI= DATAN2(FX,FY)
c
       DN1=dcos(FI)*dcos(LBD)
       DN2=dcos(FI)*dsin(LBD)
       DN3=dsin(FI)
c
       RETURN
       END
C-------------------------------------------------
       SUBROUTINE CENTRGEO(X,Y,Z,LON,LAT,H)
c---------------------------
c
c----Transformation de coordonnees geocentriques --> ellipsoidiques
c        unitÜs : mÊtres et grades
c
       IMPLICIT REAL*8 (A-H,L,O-Z)
      COMMON/COMOB/CC,EP2,PI,RAD,IFU
C
       A=CC/DSQRT(1.D0+EP2)
       B=A*A/CC
       P= DSQRT(X**2+Y**2)
       PB=B*P
       ZA=Z*A
       TETA= DATAN2(ZA,PB)
       LBD= DATAN2(Y,X)
       ST=DSIN(TETA)
       CT=DCOS(TETA)
C
       FX=Z+EP2*B*ST*ST*ST
       FY=P-EP2*A*CT*CT*CT/(1.d0+EP2)
       FI= DATAN2(FX,FY)
       CFI=DCOS(FI)
       AN=CC/DSQRT(1.D0+EP2*CFI*CFI)
c
       LAT=FI*200.D0/PI
       LON=LBD*200.D0/PI
       H= P/CFI-AN
       H1=Z/DSIN(FI)-AN/(1.d0+EP2)
       RETURN
       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/COMOB/CC,EP2,PI,RAD,IFU
       LATR=LAT*RAD
       LONR=LON*RAD
       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 GEOUTM(LONG,LAT,XX,YY)
C---------------------------------------------
C          TRANSFORMATION GEOGRAPHIQUE - U.T.M.
C            ENTREE : LONGITUDE ET LATITUDE EN GRADES
C            SORTIE :  COORD. UTM  EN METRES
C
      IMPLICIT REAL*8 (A-H,L,O-Z)
      COMMON/COMOB/CC,EP2,PI,RAD,IFU
c
      KH=1
      AK=0.9996d0
      if(IFU.LT.0) KH=-1
c
      MI=DABS(IFU)*6-183
      AMI=DBLE(MI)/0.9D0
      RD=PI/200.D0
      AM=LONG*RD
      AL=LAT*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+EP2*CO*CO)
      CX=EP2*(CO*XI)**2
      DX=GN*XI*(1.D0+CX/6.D0)
      DY=GN*TA*(1.D0+CX/2.D0)
      XX=500000.D0+AK*DX
      YY=AK*(DY+ARCME(AL))
      if(LAT) 1,2,2
    1 YY=(YY+10000000.d0)
    2 RETURN
      END
C-------------------------------------------------------------
      SUBROUTINE UTMGEO(XX,YY,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/COMOB/CC,EP2,PI,RAD,IFU
c
      KH=1
      AK=0.9996d0
      if(IFU.LT.0) KH=-1
c
      MI=DABS(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=EP2*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)*EP2*(G-L)*DCOS(L)
      LON=AM+AM0
      LAT=L+V*(G-L)
      LON=(AM+AM0)/RD
      LAT=LAT/RD*KH
       write(4,*) IFU,KH
       write(4,*) LON,LAT
c
      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/COMOB/CC,EP2,PI,RAD,IFU
c
      AQ=EP2*3.d0/4.d0
      BQ=EP2*EP2*15.d0/16.d0
      CQ=(EP2**3)*35.d0/64.d0
      DQ=(EP2**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--------------------------------------------------------------------