ll2utm Subroutine

public subroutine ll2utm(iutm, isone_in, lat, lon, utmn, utme)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: iutm

UTM coordinate system indicator
integer, intent(in) :: isone_in

UTM zone input?
real, intent(in) :: lat

Latitude in decimal degrees
real, intent(in) :: lon

Longitude in decimal degress
real, intent(out) :: utmn

UTM east-coordinate (y) (meter from west border)
real, intent(out) :: utme

UTM north-coordinate (x) (meter from equator)

Calls

proc~~ll2utm~~CallsGraph proc~ll2utm ll2utm dcos dcos proc~ll2utm->dcos dsin dsin proc~ll2utm->dsin dtan dtan proc~ll2utm->dtan

Called by

proc~~ll2utm~~CalledByGraph proc~ll2utm ll2utm proc~ll2proj LL2PROJ proc~ll2proj->proc~ll2utm proc~read_country_bounding_box_data read_country_bounding_box_data proc~read_country_bounding_box_data->proc~ll2utm proc~uemep_preaggregate_shipping_asi_data uEMEP_preaggregate_shipping_asi_data proc~uemep_preaggregate_shipping_asi_data->proc~ll2utm proc~uemep_read_industry_data uEMEP_read_industry_data proc~uemep_read_industry_data->proc~ll2utm proc~uemep_read_netcdf_population uEMEP_read_netcdf_population proc~uemep_read_netcdf_population->proc~ll2utm proc~uemep_read_receptor_data uEMEP_read_receptor_data proc~uemep_read_receptor_data->proc~ll2utm proc~uemep_read_roadlink_data_ascii uEMEP_read_roadlink_data_ascii proc~uemep_read_roadlink_data_ascii->proc~ll2utm proc~uemep_read_shipping_asi_data uEMEP_read_shipping_asi_data proc~uemep_read_shipping_asi_data->proc~ll2utm proc~uemep_set_subgrid_select_latlon_centre uEMEP_set_subgrid_select_latlon_centre proc~uemep_set_subgrid_select_latlon_centre->proc~ll2utm proc~uemep_calculate_emissions_for_emep uEMEP_calculate_emissions_for_EMEP proc~uemep_calculate_emissions_for_emep->proc~uemep_read_industry_data proc~uemep_calculate_emissions_for_emep->proc~uemep_read_roadlink_data_ascii proc~uemep_calculate_emissions_for_emep->proc~uemep_read_shipping_asi_data proc~uemep_region_mask_new uEMEP_region_mask_new proc~uemep_region_mask_new->proc~ll2proj proc~uemep_subgrid_emep_from_in_region uEMEP_subgrid_EMEP_from_in_region proc~uemep_subgrid_emep_from_in_region->proc~ll2proj program~uemep uEMEP program~uemep->proc~read_country_bounding_box_data program~uemep->proc~uemep_preaggregate_shipping_asi_data program~uemep->proc~uemep_read_industry_data program~uemep->proc~uemep_read_netcdf_population program~uemep->proc~uemep_read_receptor_data program~uemep->proc~uemep_read_roadlink_data_ascii program~uemep->proc~uemep_read_shipping_asi_data program~uemep->proc~uemep_set_subgrid_select_latlon_centre program~uemep->proc~uemep_calculate_emissions_for_emep program~uemep->proc~uemep_region_mask_new program~uemep->proc~uemep_subgrid_emep_from_in_region

Source Code

    subroutine ll2utm(iutm, isone_in, lat, lon, utmn, utme)
        integer, intent(in) :: iutm !! UTM coordinate system indicator
        integer, intent(in) :: isone_in !! UTM zone input?
        real, intent(in) :: lat !! Latitude in decimal degrees
        real, intent(in) :: lon !! Longitude in decimal degress
        real, intent(out) :: utmn !! UTM east-coordinate (y) (meter from west border)
        real, intent(out) :: utme !! UTM north-coordinate (x) (meter from equator)

        ! Local variables and parameters
        real :: isone ! UTM zone
        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) :: lon0 ! Central meridian of UTM zone
        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
        isone = abs(isone_in)
        latv = lat*pi/180.0
        lon0 = real(isone - 30.0)*6.0 - 3.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 ll2utm