LAEA2LL Subroutine

private subroutine LAEA2LL(x, y, lon, lat, projection_attr)

Arguments

Type IntentOptional Attributes Name
real, intent(in) :: x
real, intent(in) :: y
real, intent(out) :: lon
real, intent(out) :: lat
double precision, intent(in) :: projection_attr(10)

Called by

proc~~laea2ll~~CalledByGraph proc~laea2ll LAEA2LL proc~proj2ll PROJ2LL proc~proj2ll->proc~laea2ll 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 LAEA2LL(x, y, lon, lat, projection_attr)
        ! www.epsg.org
        ! grid_mapping_name = lambert_azimuthal_equal_area
        ! Map parameters:
        ! longitude_of_projection_origin
        ! latitude_of_projection_origin
        ! false_easting - This parameter is optional (default is 0)
        ! false_northing - This parameter is optional (default is 0)lat_stand1_lambert=projection_attributes(1)
        double precision, intent(in) :: projection_attr(10)
        real, intent(out) :: lon, lat
        real, intent(in) :: x, y
        real :: rho
        real :: semi_major_axis
        real :: deg2rad,rad2deg
        real :: lat0, lat0_in, lon0, lon0_in
        real :: false_easting, false_northing
        real :: inv_f, f, a, e, q_p, q_0, beta_0, R_q, D
        real :: C, beta_d
        
        lon0_in = projection_attr(1)
        lat0_in = projection_attr(2)
        false_easting = projection_attr(3)
        false_northing = projection_attr(4)
        semi_major_axis = projection_attr(5)
        inv_f = projection_attr(6) ! flattening

        deg2rad = PI/180.0
        rad2deg = 180.0/PI
        a = semi_major_axis
        f = 1.0/inv_f
        e = sqrt(2.0*f - f*f)

        lat0 = lat0_in*deg2rad
        lon0 = lon0_in*deg2rad

        q_p = (1.0 - e*e)*(1.0/(1.0 - e*e) - 1.0/(2.0*e)*log((1.0 - e)/(1.0 + e)))
        q_0 = (1.0 - e*e)*(sin(lat0)/(1.0 - e*e*sin(lat0)**2) - 1.0/(2.0*e)*log((1.0 - e*sin(lat0))/(1.0 + e*sin(lat0))))
        beta_0 = asin(q_0/q_p)
        R_q = a*sqrt(q_p/2.0)
        D = a*(cos(lat0)/sqrt(1.0 - e*e*sin(lat0)**2)/(R_q*cos(beta_0)))
        rho = sqrt((x - false_easting)*(x - false_easting)/D/D+D*D*(y - false_northing)*(y - false_northing))
        C = 2.0*asin(rho/2.0/R_q)
        beta_d = asin(cos(C)*sin(beta_0) + D*(y - false_northing)*sin(C)*cos(beta_0)/rho)

        lon = lon0 + atan2((x - false_easting)*sin(C), &
            (D*rho*cos(beta_0)*cos(C) - D*D*(y - false_northing)*sin(beta_0)*sin(C)))
        lat = beta_d + e**2*(1.0/3.0+31.0/180.0*e**2.0 + 517.0/5040.0*e**4)*sin(2.0*beta_d) &
            + e**4.0*(23.0/360.0+251.0/3780.0*e**2)*sin(4.0*beta_d)+(761.0/45360.0*e**6)*sin(6.0*beta_d)

        lat = lat*rad2deg
        lon = lon*rad2deg
    end subroutine LAEA2LL