LAEA2LL_spherical Subroutine

    subroutine LAEA2LL_spherical(x, y, lon, lat, projection_attr)
        ! https://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html
        ! 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

        ! Local variables
        real :: r, rho, c
        real :: earth_radius
        real :: deg2rad, rad2deg
        real :: lat0, lat0_in, lon0, lon0_in
        real :: false_easting, false_northing
        
        lon0_in = projection_attr(1)
        lat0_in = projection_attr(2)
        false_easting = projection_attr(3)
        false_northing = projection_attr(4)
        earth_radius = projection_attr(5)

        deg2rad = PI/180.0
        rad2deg = 180.0/PI
        r = earth_radius

        lat0 = lat0_in*deg2rad
        lon0 = lon0_in*deg2rad

        rho = sqrt((x - false_easting)*(x - false_easting) + (y - false_northing)*(y - false_northing))
        c = 2.0*asin(rho*0.5/r)
        lat = asin(cos(c)*sin(lat0)+(y-false_northing)/rho*sin(c)*cos(lat0))
        lon = lon0 + atan((x - false_easting)*sin(c)/(rho*cos(lat0)*cos(c) - (y - false_northing)*sin(lat0)*sin(c)))
        lat = lat*rad2deg
        lon = lon*rad2deg
    end subroutine LAEA2LL_spherical