ltm2ll Subroutine

public subroutine ltm2ll(iutm, isone_in, la0, utmn_in, utme, lat, lon)

Local Transverse Mecantor version of utm2ll

Arguments

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

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

UTM zone
real, intent(in) :: la0

Tangeringsmeridian
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

Calls

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

Called by

proc~~ltm2ll~~CalledByGraph proc~ltm2ll ltm2ll proc~proj2ll PROJ2LL proc~proj2ll->proc~ltm2ll proc~read_country_bounding_box_data read_country_bounding_box_data proc~read_country_bounding_box_data->proc~proj2ll proc~uemep_define_subgrid uEMEP_define_subgrid proc~uemep_define_subgrid->proc~proj2ll proc~uemep_define_subgrid_extent uEMEP_define_subgrid_extent proc~uemep_define_subgrid_extent->proc~proj2ll proc~uemep_grid_roads uEMEP_grid_roads proc~uemep_grid_roads->proc~proj2ll proc~uemep_preaggregate_shipping_asi_data uEMEP_preaggregate_shipping_asi_data proc~uemep_preaggregate_shipping_asi_data->proc~proj2ll proc~uemep_read_emep uEMEP_read_EMEP proc~uemep_read_emep->proc~proj2ll proc~uemep_read_landuse_rivm_data uEMEP_read_landuse_rivm_data proc~uemep_read_landuse_rivm_data->proc~proj2ll proc~uemep_read_meteo_nc uEMEP_read_meteo_nc proc~uemep_read_meteo_nc->proc~proj2ll proc~uemep_read_monthly_and_daily_shipping_asi_data uEMEP_read_monthly_and_daily_shipping_asi_data proc~uemep_read_monthly_and_daily_shipping_asi_data->proc~proj2ll proc~uemep_read_netcdf_landuse_latlon uEMEP_read_netcdf_landuse_latlon proc~uemep_read_netcdf_landuse_latlon->proc~proj2ll proc~uemep_read_netcdf_population_latlon uEMEP_read_netcdf_population_latlon proc~uemep_read_netcdf_population_latlon->proc~proj2ll proc~uemep_read_netcdf_shipping_latlon uEMEP_read_netcdf_shipping_latlon proc~uemep_read_netcdf_shipping_latlon->proc~proj2ll proc~uemep_read_roadlink_data_ascii uEMEP_read_roadlink_data_ascii proc~uemep_read_roadlink_data_ascii->proc~proj2ll proc~uemep_read_weekly_shipping_asi_data uEMEP_read_weekly_shipping_asi_data proc~uemep_read_weekly_shipping_asi_data->proc~proj2ll proc~uemep_region_mask_new uEMEP_region_mask_new proc~uemep_region_mask_new->proc~proj2ll proc~uemep_set_region_tile_grids uEMEP_set_region_tile_grids proc~uemep_set_region_tile_grids->proc~proj2ll proc~uemep_set_tile_grids uEMEP_set_tile_grids proc~uemep_set_tile_grids->proc~proj2ll proc~uemep_calculate_emissions_for_emep uEMEP_calculate_emissions_for_EMEP proc~uemep_calculate_emissions_for_emep->proc~uemep_grid_roads proc~uemep_calculate_emissions_for_emep->proc~uemep_read_meteo_nc proc~uemep_calculate_emissions_for_emep->proc~uemep_read_monthly_and_daily_shipping_asi_data proc~uemep_calculate_emissions_for_emep->proc~uemep_read_roadlink_data_ascii proc~uemep_calculate_emissions_for_emep->proc~uemep_read_weekly_shipping_asi_data program~uemep uEMEP program~uemep->proc~read_country_bounding_box_data program~uemep->proc~uemep_define_subgrid program~uemep->proc~uemep_define_subgrid_extent program~uemep->proc~uemep_grid_roads program~uemep->proc~uemep_preaggregate_shipping_asi_data program~uemep->proc~uemep_read_emep program~uemep->proc~uemep_read_landuse_rivm_data program~uemep->proc~uemep_read_meteo_nc program~uemep->proc~uemep_read_monthly_and_daily_shipping_asi_data program~uemep->proc~uemep_read_netcdf_landuse_latlon program~uemep->proc~uemep_read_netcdf_population_latlon program~uemep->proc~uemep_read_netcdf_shipping_latlon program~uemep->proc~uemep_read_roadlink_data_ascii program~uemep->proc~uemep_read_weekly_shipping_asi_data program~uemep->proc~uemep_region_mask_new program~uemep->proc~uemep_set_region_tile_grids program~uemep->proc~uemep_set_tile_grids program~uemep->proc~uemep_calculate_emissions_for_emep

Source Code

    subroutine ltm2ll(iutm, isone_in, la0, utmn_in, utme, lat, lon)
        !! Local Transverse Mecantor version of utm2ll
        integer, intent(in) :: iutm !! UTM coordinate system indicator
        integer, intent(in) :: isone_in !! UTM zone
        real, intent(in) :: la0 !! Tangeringsmeridian
        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) :: x ! Scaled north-coordinate
        real(dp) :: y ! Scaled east-coordinate
        real(dp) :: bb0, e, e2, fi, 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
         
        ! djust 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

        ! 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 ltm2ll