ll2ltm Subroutine

    subroutine ll2ltm(iutm, lon0, lat, lon, utmn, utme)
        !! Local lon version (of ll2utm) without zone, so just typical Transverse
        !! Mecantor (Local Transverse Mecantor)
        integer, intent(in) :: iutm !! UTM coordinate system indicator
        real, intent(in) :: lon0 !! Central meridian of UTM zone
        real, intent(in) :: lat !! Latitude in decimal degrees
        real, intent(in) :: lon !! Longitude in decimal degrees
        real, intent(out) :: utmn !! UTM north-coordinate (x) (meter from equator)
        real, intent(out) :: utme !! UTM  east-coordinate (y) (meter from west border)

        ! Local variables and paramters
        real(dp), parameter :: deast = 500000.0 ! East movement UTM
        real(dp), parameter :: scale = 0.9996 ! Scale UTM
        real(dp) :: a ! Big semiaxis
        real(dp) :: f ! Flattening
        real(dp) :: latv ! Scaled latitude
        real(dp) :: lonv ! Scaled longitude
        real(dp) :: b, e, e2, m, n ! Intermediate values

        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

        ! Scale coordinates
        latv = lat*pi/180.0
        lonv = (lon - lon0)*pi/180.0

        ! Calculate some intermediate quantities
        e2 = f*(2.0 - f)
        n = a/dsqrt(1.0 - e2*dsin(latv)*dsin(latv))
        e = dsqrt(e2*dcos(latv)*dcos(latv)/(1.0 - e2))
        m = n/(1.0 + e*e)
        b = ((1.0 - f/2.0 + f*f/16.0 + f*f*f/32.0)*latv &
            - (3.0*f/4.0 - 3.0*f*f*f/128.0)*dsin(2.0*latv) &
            + (15.0*f*f/64.0 + 15.0*f*f*f/128.0)*dsin(4.0*latv) &
            - (35.0*f*f*f/384.0)*dsin(6.0*latv))*a

        ! Calculate UTM north coordinate

        utmn = (b + lonv*lonv*n*dsin(latv)*dcos(latv)/2.0 &
            + lonv*lonv*lonv*lonv*n*dsin(latv) &
            * dcos(latv)*dcos(latv)*dcos(latv) &
            * (5.0 - dtan(latv)*dtan(latv) + 9.0*e*e + 4.0*e*e*e*e) &
            / 24.0 + lonv*lonv*lonv*lonv*lonv*lonv*n*dsin(latv) &
            * dcos(latv)*dcos(latv)*dcos(latv)*dcos(latv)*dcos(latv) &
            * (61.0 - 58.0*dtan(latv)*dtan(latv) &
            + dtan(latv)*dtan(latv)*dtan(latv)*dtan(latv)) &
            / 720.0)*scale
      
        if (lat < 0.0) then
            utmn = utmn + 10000000.0
        end if

        ! Calculate UTM east coordinate
        utme = (lonv*n*dcos(latv) &
            + lonv*lonv*lonv*n*dcos(latv)*dcos(latv)*dcos(latv) &
            * (1.0 - dtan(latv)*dtan(latv) + e*e)/6.0 &
            + lonv*lonv*lonv*lonv*lonv*n &
            * dcos(latv)*dcos(latv)*dcos(latv)*dcos(latv)*dcos(latv) &
            * (5.0 - 18.0*dtan(latv)*dtan(latv) &
            + dtan(latv)*dtan(latv)*dtan(latv)*dtan(latv))/120.0) &
            * scale + deast
    end subroutine ll2ltm