conv_lambert_to_degdec.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\conv_lambert_to_degdec.f
00002 !
00055                         SUBROUTINE CONV_LAMBERT_TO_DEGDEC
00056 !                       *********************************
00057 !
00058      &(NTAB,XTAB,YTAB,LAMBDATAB,PHITAB,NUMZONE)
00059 !
00060 !***********************************************************************
00061 ! TELEMAC2D   V6P2                                   25/06/2010
00062 !***********************************************************************
00063 !
00064 !
00065 !
00066 !
00067 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00068 !| LAMBDATAB      |<--| LONGITUDE
00069 !| NUMZONE        |-->| NUMBER OF LAMBERT ZONE
00070 !| NTAB           |-->| NUMBER OF COORDINATES
00071 !| PHITAB         |<--| LATITUDE
00072 !| XTAB           |-->| METRIC COORDINATES (LAMBERT)
00073 !| YTAB           |-->| METRIC COORDINATES (LAMBERT)
00074 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00075 !
00076       USE INTERFACE_TELEMAC2D, EX_CONV_LAMBERT_TO_DEGDEC
00077      &                         => CONV_LAMBERT_TO_DEGDEC
00078 !
00079       IMPLICIT NONE
00080       INTEGER LNG,LU
00081       COMMON/INFO/LNG,LU
00082 !
00083 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00084 !
00085       INTEGER, INTENT(IN) :: NTAB,NUMZONE
00086       DOUBLE PRECISION, INTENT(IN)  :: XTAB(NTAB),YTAB(NTAB)
00087       DOUBLE PRECISION, INTENT(OUT) :: LAMBDATAB(NTAB),PHITAB(NTAB)
00088 !
00089 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00090 !
00091       DOUBLE PRECISION PI,DTR,RTD
00092       DOUBLE PRECISION X,Y,LAMBDA,PHI
00093       DOUBLE PRECISION LAMBDAC,EEE,ES2,NNN,CCC,XS,YS,RRR,GAMMA,LATISO
00094       DOUBLE PRECISION HE,AAA,GNORM,Z,TX,TY,TZ,PHIM,EPSI,CORRPHI
00095 !
00096       INTEGER I,J
00097 !
00098 !-----------------------------------------------------------------------
00099 !
00100       PI  = 4.D0*ATAN(1.D0)
00101       DTR = PI/180.D0
00102       RTD = 180.D0/PI
00103 !
00104       EPSI = 1.D-11
00105 !
00106 !  LAMBDAC : PARIS MERIDIAN / GREENWICH MERIDIAN
00107 !
00108       LAMBDAC = (2.D0+20.D0/60.D0+14.025D0/3600.D0)*DTR
00109 !
00110 !  TRANSF PROJECTION CC LAMBERT X, Y --> COORD GEO LAMBERT LAMBDA, PHI
00111 !  (IGN: ALG0004)
00112 !
00113       EEE = 0.08248325676D0
00114 !
00115       ES2 = EEE/2.D0
00116 !
00117 !  NNN, CCC, XS, YS IN CONSTANTES DE PROJECTION, PROJECTION LAMBERT FRANCE
00118 !
00119 !  LAMBERT 1
00120       IF(NUMZONE.EQ.1) THEN
00121 !       NNN = 0.760405966D0
00122 !       CCC = 11603796.9767D0
00123         NNN = 0.7604059656D0
00124         CCC = 11603796.98D0
00125         XS  = 600000.D0
00126         YS  = 5657616.674D0
00127 !  LAMBERT 2
00128       ELSEIF(NUMZONE.EQ.2) THEN
00129         NNN = 0.7289686274D0
00130         CCC = 11745793.39D0
00131         XS  = 600000.D0
00132         YS  = 6199695.768D0
00133 !  LAMBERT 3
00134       ELSEIF(NUMZONE.EQ.3) THEN
00135         NNN = 0.6959127966D0
00136         CCC = 11947992.52D0
00137         XS  = 600000.D0
00138         YS  = 6791905.085D0
00139 !  LAMBERT 4
00140       ELSEIF(NUMZONE.EQ.4) THEN
00141         NNN = 0.6712679322D0
00142         CCC = 12136281.99D0
00143         XS  = 234.358D0
00144         YS  = 7239161.542D0
00145 !  LAMBERT 2 ETENDU
00146       ELSEIF(NUMZONE.EQ.22) THEN
00147         NNN = 0.7289686274D0
00148         CCC = 11745793.39D0
00149         XS  = 600000.D0
00150         YS  = 8199695.768D0
00151 !  LAMBERT 93
00152       ELSEIF(NUMZONE.EQ.93) THEN
00153         IF(LNG.EQ.1) THEN
00154           WRITE(LU,*) 'CONVERSION DE LAMBERT 93 PAS ENCORE IMPLEMENTEE'
00155         ENDIF
00156         IF(LNG.EQ.2) THEN
00157           WRITE(LU,*) 'CONVERSION FROM LAMBERT 93 NOT YET IMPLEMENTED'
00158         ENDIF
00159         CALL PLANTE(1)
00160         STOP
00161 !  DEFAULT VALUE
00162       ELSEIF(NUMZONE.EQ.-1) THEN
00163         IF(LNG.EQ.1) THEN
00164           WRITE(LU,*) 'VALEUR PAR DEFAUT INCORRECTE POUR LA PROJECTION'
00165           WRITE(LU,*) 'LAMBERT. CHOISIR PARMI LES CHOIX POSSIBLES :'
00166           WRITE(LU,*) '  -1  : LAMBERT 1 NORD ;'
00167           WRITE(LU,*) '  -2  : LAMBERT 2 CENTRE ;'
00168           WRITE(LU,*) '  -3  : LAMBERT 3 SUD ;'
00169           WRITE(LU,*) '  -4  : LAMBERT 4 CORSE ;'
00170           WRITE(LU,*) '  -22 : LAMBERT 2 ETENDU.'
00171         ENDIF
00172         IF(LNG.EQ.2) THEN
00173           WRITE(LU,*) 'INCORRECT DEFAULT VALUE FOR LAMBERT PROJECTION.'
00174           WRITE(LU,*) 'TO BE CHOSEN AMONG THE POSSIBLE CHOICES:'
00175           WRITE(LU,*) '  -1 : LAMBERT 1 NORTH ;'
00176           WRITE(LU,*) '  -2 : LAMBERT 2 CENTER ;'
00177           WRITE(LU,*) '  -3 : LAMBERT 3 SOUTH ;'
00178           WRITE(LU,*) '  -4 : LAMBERT 4 CORSICA ;'
00179           WRITE(LU,*) '  -22: LAMBERT 2 EXTENDED.'
00180         ENDIF
00181         CALL PLANTE(1)
00182         STOP
00183       ELSE
00184         IF(LNG.EQ.1) THEN
00185           WRITE(LU,*) 'NUMERO DE PROJECTION LAMBERT',NUMZONE,'INCONNUE.'
00186           WRITE(LU,*) 'CHOISIR PARMI LES CHOIX POSSIBLES :'
00187           WRITE(LU,*) '  -1  : LAMBERT 1 NORD ;'
00188           WRITE(LU,*) '  -2  : LAMBERT 2 CENTRE ;'
00189           WRITE(LU,*) '  -3  : LAMBERT 3 SUD ;'
00190           WRITE(LU,*) '  -4  : LAMBERT 4 CORSE ;'
00191           WRITE(LU,*) '  -22 : LAMBERT 2 ETENDU.'
00192         ENDIF
00193         IF(LNG.EQ.2) THEN
00194           WRITE(LU,*) 'UNKNOWN NUMBER OF LAMBERT PROJECTION',NUMZONE
00195           WRITE(LU,*) 'TO BE CHOSEN AMONG THE POSSIBLE CHOICES:'
00196           WRITE(LU,*) '  -1 : LAMBERT 1 NORTH ;'
00197           WRITE(LU,*) '  -2 : LAMBERT 2 CENTER ;'
00198           WRITE(LU,*) '  -3 : LAMBERT 3 SOUTH ;'
00199           WRITE(LU,*) '  -4 : LAMBERT 4 CORSICA ;'
00200           WRITE(LU,*) '  -22: LAMBERT 2 EXTENDED.'
00201         ENDIF
00202         CALL PLANTE(1)
00203         STOP
00204       ENDIF
00205 !
00206       DO J=1,NTAB
00207         X = XTAB(J)
00208         Y = YTAB(J)
00209         RRR = SQRT((X-XS)**2+(Y-YS)**2)
00210         GAMMA = ATAN((X-XS)/(YS-Y))
00211 !
00212         LAMBDA = GAMMA/NNN+LAMBDAC
00213 !
00214         LATISO = -1.D0/NNN*LOG(ABS(RRR/CCC))
00215 !
00216 !  CALCULATION OF LATITUDE PHI FROM ISOMETRIC LATITUDE LATISO
00217 !  (IGN: ALG0002)
00218 !
00219 !  I = 0
00220         PHIM = 2.D0*ATAN(EXP(LATISO))-PI/2.D0
00221 !  I = 1
00222         PHI  = 2.D0*ATAN(EXP(LATISO)*( (1.D0+EEE*SIN(PHIM))
00223      &                                /(1.D0-EEE*SIN(PHIM)))**ES2)
00224      &        -PI/2.D0
00225 !
00226         I = 1
00227 !
00228         DO WHILE (ABS(PHI-PHIM).GE.EPSI)
00229           PHIM = PHI
00230           PHI  = 2.D0*ATAN(EXP(LATISO)*( (1.D0+EEE*SIN(PHIM))
00231      &                                  /(1.D0-EEE*SIN(PHIM)))**ES2)
00232      &          -PI/2.D0
00233           I = I + 1
00234         ENDDO
00235 !
00236 !  TRANSFORMATION GEOGRAPHIC COORDINATES LAMBDA, PHI --> CARTESIAN COORDINATES X, Y, Z
00237 !  (IGN: ALG0009)
00238 !
00239 !  COMPUTATION OF THE BIG NORMAL (GRANDE NORMALE) GNORM
00240 !  LAT PHI --> GNORM (IGN: ALG0021)
00241 !
00242         HE = 0.D0               ! FOR THE TEST, NOT USEFUL FINALLY
00243 !
00244         AAA = 6378249.2D0
00245         EEE = 0.08248325679D0
00246 !
00247         GNORM = AAA/SQRT(1.D0-EEE**2*SIN(PHI)**2)
00248 !
00249         X = (GNORM+HE)*COS(PHI)*COS(LAMBDA)
00250         Y = (GNORM+HE)*COS(PHI)*SIN(LAMBDA)
00251         Z = (GNORM*(1.D0-EEE**2)+HE)*SIN(PHI)
00252 !
00253 !  COORDINATE TRANSFORMATION AT 7 PARAMETERS BETWEEN 2 GEODESIC SYSTEMS
00254 !  SIMPLIFIED TRANSFORMATION, JUST A TRANSLATION WITH 3 PARAMETERS (IGN: ALG0013 SIMPLIFIED)
00255 !  NTF LAMBERT --> WGS84 (IGN)
00256 !
00257         TX = -168.D0
00258         TY = -60.D0
00259         TZ =  320.D0
00260 !
00261         X = X + TX
00262         Y = Y + TY
00263         Z = Z + TZ
00264 !
00265 !  TRANSFORMATION CARTESIAN COORDINATES X, Y, Z --> GEOGRAPHIC COORDINATES LAMBDA, PHI
00266 !  HEISKANEN-MORITZ-BOUCHER METHOD (IGN: ALG0012)
00267 !
00268         LAMBDA = ATAN(Y/X)
00269 !  I = 0
00270         PHIM   = ATAN(Z/(SQRT(X**2+Y**2)*(1.D0-(AAA*EEE**2)
00271      &                                        /(SQRT(X**2+Y**2+Z**2)))))
00272 !  I = 1
00273         PHI = ATAN(Z/(SQRT(X**2+Y**2)*(1.D0-(AAA*EEE**2*COS(PHIM))
00274      &              /(SQRT(X**2+Y**2)*SQRT(1.D0-EEE**2*SIN(PHIM)**2)))))
00275 !
00276         I = 1
00277 !
00278         DO WHILE (ABS(PHI-PHIM).GE.EPSI)
00279           PHIM = PHI
00280           PHI = ATAN(Z/(SQRT(X**2+Y**2)*(1.D0-(AAA*EEE**2*COS(PHIM))
00281      &              /(SQRT(X**2+Y**2)*SQRT(1.D0-EEE**2*SIN(PHIM)**2)))))
00282           I = I + 1
00283         ENDDO
00284 !
00285 !  ANGLE CORRECTION FOR PHI (RE-CALIBRATION)
00286 !
00287         CORRPHI = -5.439609D-5
00288         PHI = PHI + CORRPHI
00289 !
00290 !  CONVERSION INTO DECIMAL DEGREES
00291 !
00292         LAMBDA = LAMBDA*RTD
00293         PHI    = PHI*RTD
00294 !
00295         LAMBDATAB(J) = LAMBDA
00296         PHITAB(J)    = PHI
00297 !
00298       ENDDO
00299 !
00300 !-----------------------------------------------------------------------
00301 !
00302       RETURN
00303       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0