The subroutine converts UTM north- and east-coordinates to latitude and longitude
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | iutm |
UTM coordinate system indicator |
||
integer, | intent(in) | :: | isone_in |
UTM zone |
||
real, | intent(in) | :: | utmn_in |
UTM north-coordinate (X) (meter from equator) |
||
real, | intent(in) | :: | utme |
UTM east-coordinate (Y) (meter from west border) |
||
real, | intent(out) | :: | lat |
Latitude in decimal degrees |
||
real, | intent(out) | :: | lon |
Longitude in decimal degrees |
subroutine utm2ll(iutm, isone_in, utmn_in, utme, lat, lon) !! The subroutine converts UTM north- and east-coordinates to latitude !! and longitude integer, intent(in) :: iutm !! UTM coordinate system indicator integer, intent(in) :: isone_in !! UTM zone real, intent(in) :: utmn_in !! UTM north-coordinate (X) (meter from equator) real, intent(in) :: utme !! UTM east-coordinate (Y) (meter from west border) real, intent(out) :: lat !! Latitude in decimal degrees real, intent(out) :: lon !! Longitude in decimal degrees ! Local variables and parameters integer :: isone real :: utmn ! UTM north-coordinate real(dp), parameter :: deast = 500000.0 real(dp), parameter :: scale = 0.9996 real(dp) :: a ! Big semiaxis real(dp) :: f ! Flattening real(dp) :: la0 ! Tangeringsmeridian real(dp) :: x ! Scaled north-coordinate real(dp) :: y ! Scaled east-coordinate real(dp) :: bb0, e, e2, fi, m, n ! Intermediate value select case(iutm) case(1) ! UTM WGS84 EUREF 89 (AirQUIS) a = 6378137.0 f = 1.0/298.257222101 case(2) ! UTM WGS84 OLD a = 6378137.0 f = 1.0/298.257223563 case(3) ! UTM ED50 a = 6378388.0 f = 1.0/297.0 case default print *, "ERROR: Unknown UTM coordinate system indictor, iutm: ", iutm stop 1 end select ! Adjust for southern hemisphere, specified by negative isone_in if (isone_in < 0) then utmn = utmn_in - 10000000.0 else utmn = utmn_in end if isone = abs(isone_in) ! Scale coordinates x = utmn/scale y = (utme - deast)/scale la0 = real(isone - 30)*6.0 - 3.0 ! Calculate some intermediate quantities e2 = f*(2.0 - f) bb0 = (1.0 - f/2.0 + f*f/16.0 + f*f*f/32.0)*a fi = x/bb0 & + (3.0*f/4.0 + 3.0*f*f/8.0 + 21.0*f*f*f/256.0)*dsin(2.0*x/bb0) & + (21.0*f*f/64.0 + 21.0*f*f*f/64.0)*dsin(4.0*x/bb0) & + (151.0*f*f*f/768.0)*dsin(6.0*x/bb0) n = a/dsqrt(1.0 - e2*dsin(fi)*dsin(fi)) e = dsqrt(e2*dcos(fi)*dcos(fi)/(1.0 - e2)) m = n/(1.0 + e*e) ! Calculate latitude and longitude in radians lat = fi - (y*y*dtan(fi))/(2.0*m*n) & + (y*y*y*y*dtan(fi))/(24.0*m*n*n*n) & * (5.0 + 3.0*dtan(fi)*dtan(fi) + e*e & - 9.0*e*e*dtan(fi)*dtan(fi) - 4.0*e*e*e*e) & - (y*y*y*y*y*y*dtan(fi))/(720.0*m*n*n*n*n*n) & * (61.0 + 90.0*dtan(fi)*dtan(fi) & + 45.0*dtan(fi)*dtan(fi)*dtan(fi)*dtan(fi)) lon = y/(n*dcos(fi)) & - (y*y*y*(1.0 + 2.0*dtan(fi)*dtan(fi) + e*e)) & / (6.0*n*n*n*dcos(fi)) & + (y*y*y*y*y*(5.0 + 28.0*dtan(fi)*dtan(fi) & + 24.0*dtan(fi)*dtan(fi)*dtan(fi)*dtan(fi))) & / (120.0*n*n*n*n*n*dcos(fi)) + la0*pi/180.0 ! Convert from radians to degrees lat = lat*180.0/pi lon = lon*180.0/pi end subroutine utm2ll