conv_mercator_to_degdec.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\conv_mercator_to_degdec.f
00002 !
00063                       SUBROUTINE CONV_MERCATOR_TO_DEGDEC
00064 !                     **********************************
00065 !
00066      &(NTAB,XTAB,YTAB,LAMBDATAB,PHITAB,GEOSYST,NUMZONE,LONG0,LAT0)
00067 !
00068 !***********************************************************************
00069 ! TELEMAC2D   V6P2                                   22/04/2011
00070 !***********************************************************************
00071 !
00072 !
00073 !
00074 !
00075 !
00076 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00077 !| LAMBDATAB      |<--| LONGITUDE (DECIMAL DEGREES)
00078 !| LAT0           |-->| LATITUDE OF ORIGIN POINT (KEYWORD, IN DEGREES)
00079 !|                |   | BEWARE!!! IN TELEMAC PHI0 IS LONGITUDE!!!
00080 !| LONG0          |-->| LONGITUDE OF ORIGIN POINT (KEYWORD, IN DEGREES)
00081 !|                |   | BEWARE!!! IN TELEMAC LAMBD0 IS LATITUDE!!!
00082 !| NTAB           |-->| NUMBER OF COORDINATES
00083 !| GEOSYST        |-->| TYPE OF GEOGRAPHIC SYSTEM
00084 !|                |   | 2: UTM NORTH     3: UTM SOUTH
00085 !|                |   | 5: MERCATOR FOR TELEMAC
00086 !| NUMZONE        |-->| NUMBER OF UTM ZONE
00087 !| PHITAB         |<--| LATITUDE (DECIMAL DEGREES)
00088 !| XTAB           |-->| METRIC COORDINATES (WGS84 UTM)
00089 !| YTAB           |-->| METRIC COORDINATES (WGS84 UTM)
00090 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00091 !
00092       USE INTERFACE_TELEMAC2D, EX_CONV_MERCATOR_TO_DEGDEC
00093      &                         => CONV_MERCATOR_TO_DEGDEC
00094 !
00095       IMPLICIT NONE
00096       INTEGER LNG,LU
00097       COMMON/INFO/LNG,LU
00098 !
00099 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00100 !
00101       INTEGER, INTENT(IN)           :: NTAB,GEOSYST,NUMZONE
00102       DOUBLE PRECISION, INTENT(IN)  :: LAT0,LONG0
00103       DOUBLE PRECISION, INTENT(IN)  :: XTAB(NTAB),YTAB(NTAB)
00104       DOUBLE PRECISION, INTENT(OUT) :: LAMBDATAB(NTAB),PHITAB(NTAB)
00105 !
00106 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00107 !
00108       DOUBLE PRECISION PI,DTR,RTD,CONST,AAA,FFF,EEE,EEE2,EEE4,EEE6,EEE8
00109       DOUBLE PRECISION LAMBDAC,NNN,XS,YS,PHIM,LATISO,LATISOS,ES2,EPSI
00110       DOUBLE PRECISION X,Y,LAMBDA,PHI
00111       DOUBLE PRECISION CITM(5)
00112       DOUBLE PRECISION, PARAMETER :: RADIUS = 6371000.D0
00113       COMPLEX(KIND(1.D0)) ZPRIME,ZZZ
00114 !
00115       INTEGER I,J,K
00116 !
00117 !-----------------------------------------------------------------------
00118 !
00119       PI  = 4.D0*ATAN(1.D0)
00120       DTR = PI/180.D0
00121       RTD = 180.D0/PI
00122       CONST = TAN(0.5D0*LAT0*DTR+0.25D0*PI)
00123 !
00124       EPSI = 1.D-11
00125 !
00126       AAA = 6378137.D0
00127       FFF = 1.D0/298.257223563D0
00128       EEE = SQRT(2.D0*FFF-FFF**2)
00129 !     EEE = 0.081991889980000D0
00130 !
00131       IF(GEOSYST.EQ.2.AND.NUMZONE.EQ.-1) THEN
00132 !     NORTHERN UTM
00133         IF(LNG.EQ.1) THEN
00134           WRITE(LU,*) 'VALEUR PAR DEFAUT INCORRECTE POUR LE NUMERO'
00135           WRITE(LU,*) 'DE ZONE UTM. A CHOISIR ENTRE 1 ET 60.'
00136           WRITE(LU,*) 'PAR EXEMPLE ENTRE 30 ET 32 EN FRANCE.'
00137         ENDIF
00138         IF(LNG.EQ.2) THEN
00139           WRITE(LU,*) 'INCORRECT DEFAULT VALUE FOR UTM ZONE.'
00140           WRITE(LU,*) 'TO BE CHOSEN BETWEEN 1 AND 60.'
00141           WRITE(LU,*) 'E.G. BETWEEN 30 AND 32 FOR FRANCE.'
00142         ENDIF
00143         CALL PLANTE(1)
00144         STOP
00145       ELSEIF(GEOSYST.EQ.3.AND.NUMZONE.EQ.-1) THEN
00146 !     SOUTHERN UTM
00147         IF(LNG.EQ.1) THEN
00148           WRITE(LU,*) 'VALEUR PAR DEFAUT INCORRECTE POUR LE NUMERO'
00149           WRITE(LU,*) 'DE ZONE UTM. A CHOISIR ENTRE 1 ET 60.'
00150         ENDIF
00151         IF(LNG.EQ.2) THEN
00152           WRITE(LU,*) 'INCORRECT DEFAULT VALUE FOR UTM ZONE.'
00153           WRITE(LU,*) 'TO BE CHOSEN BETWEEN 1 AND 60.'
00154         ENDIF
00155         CALL PLANTE(1)
00156         STOP
00157       ENDIF
00158 !
00159       NNN = 0.9996D0 * AAA
00160       LAMBDAC = (6.D0*REAL(NUMZONE)-183.D0)*DTR ! RADIANS
00161       XS = 500000.D0
00162 !
00163 !     NORTHERN UTM
00164 !
00165       IF(GEOSYST.EQ.2) THEN
00166         YS = 0.D0
00167 !       SOUTHERN UTM
00168       ELSEIF(GEOSYST.EQ.3) THEN
00169         YS = 10000000.D0
00170       ENDIF
00171 !
00172       ES2 = EEE/2.D0
00173 !
00174       EEE2 = EEE**2
00175       EEE4 = EEE2**2
00176       EEE6 = EEE2*EEE4
00177       EEE8 = EEE4**2
00178 !
00179 !  PROJECTION COEFFICIENTS
00180 !  MERCATOR TRANSVERSE PROJECTION "BACKWARD"
00181 !  (IGN: ALG0029)
00182 !
00183       CITM(1) = 1.D0 - 1.D0/4.D0*EEE2 - 3.D0/64.D0*EEE4
00184      &               - 5.D0/256.D0*EEE6 - 175.D0/16384.D0*EEE8
00185       CITM(2) = 1.D0/8.D0*EEE2 + 1.D0/48.D0*EEE4 + 7.D0/2048.D0*EEE6
00186      &                         + 1.D0/61440.D0*EEE8
00187       CITM(3) = 1.D0/768.D0*EEE4 + 3.D0/1280.D0*EEE6
00188      &                           + 559.D0/368640.D0*EEE8
00189       CITM(4) = 17.D0/30720.D0*EEE6 + 283.D0/430080.D0*EEE8
00190       CITM(5) = 4397.D0/41287680.D0*EEE8
00191 !
00192       IF(GEOSYST.EQ.2.OR.GEOSYST.EQ.3) THEN
00193 !
00194 !  BEGINNING OF LOOP ON POINTS
00195 !
00196       DO J=1,NTAB
00197         X = XTAB(J)
00198         Y = YTAB(J)
00199 !       ZPRIME = DCMPLX((Y-YS), (X-XS))/NNN/CITM(1)
00200         ZPRIME = CMPLX(Y-YS,X-XS,KIND(1.D0))/NNN/CITM(1)
00201 
00202         ZZZ = ZPRIME
00203 
00204         DO K=1,4
00205           ZZZ = ZZZ - CITM(K+1)*SIN(2.D0*REAL(K)*ZPRIME)
00206         ENDDO
00207 
00208         LATISO  =  REAL(ZZZ)
00209         LATISOS = AIMAG(ZZZ)
00210 !
00211         LAMBDA = LAMBDAC + ATAN(SINH(LATISOS)/COS(LATISO))
00212         PHI    =           ASIN(SIN(LATISO)/COSH(LATISOS))
00213 !
00214 !  COMPUTATION OF LATITUDE PHI FROM ISOMETRIC LATITUDE LATISO
00215 !  (ELLIPSOID OF 1ST EXCENTRICITY EEE AT POINT OF LATITUDE PHI)
00216 !  (IGN: ALG0001)
00217 !
00218         LATISO = LOG(TAN(PI/4.D0+PHI/2.D0)) ! EEE = 0.D0 HERE
00219 !
00220 !  COMPUTATION OF LATITUDE PHI FROM ISOMETRIC LATITUDE LATISO
00221 !  (IGN: ALG0002)
00222 !
00223 !  I = 0
00224         PHIM = 2.D0*ATAN(EXP(LATISO))-PI/2.D0
00225 !  I = 1
00226         PHI  = 2.D0*ATAN(EXP(LATISO)*( (1.D0+EEE*SIN(PHIM))
00227      &                                /(1.D0-EEE*SIN(PHIM)))**ES2)
00228      &        -PI/2.D0
00229 !
00230         I = 1
00231 !
00232         DO WHILE (ABS(PHI-PHIM).GE.EPSI)
00233           PHIM = PHI
00234           PHI  = 2.D0*ATAN(EXP(LATISO)*( (1.D0+EEE*SIN(PHIM))
00235      &                                  /(1.D0-EEE*SIN(PHIM)))**ES2)
00236      &          -PI/2.D0
00237           I = I + 1
00238         ENDDO
00239 !
00240 !  CONVERSION INTO DECIMAL DEGREES
00241 !
00242         LAMBDA = LAMBDA*RTD
00243         PHI    = PHI*RTD
00244 !
00245         LAMBDATAB(J) = LAMBDA
00246         PHITAB(J)    = PHI
00247       ENDDO
00248 !
00249       ELSEIF(GEOSYST.EQ.5) THEN
00250         DO J=1,NTAB
00251           LAMBDATAB(J) = XTAB(J)/RADIUS*RTD + LONG0
00252           PHITAB(J)    = ( 2.D0*ATAN(CONST*EXP(YTAB(J)/RADIUS))
00253      &                    -0.5D0*PI)*RTD
00254         ENDDO
00255       ENDIF
00256 !
00257 !-----------------------------------------------------------------------
00258 !
00259       RETURN
00260       END

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