LL2LAEA_spherical Subroutine

    subroutine LL2LAEA_spherical(x, y, lon_in, lat_in, 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(in) :: lon_in, lat_in
        real, intent(out):: x, y

        ! Local variables
        real :: r
        real :: earth_radius
        real :: deg2rad, rad2deg, k_lambert
        real :: lat0, lat0_in, lon0, lon0_in
        real :: false_easting, false_northing
        real :: lon, lat

        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
        lon = lon_in*deg2rad
        lat = lat_in*deg2rad

        k_lambert = sqrt(2.0/(1.0 + sin(lat0)*sin(lat) + cos(lat0)*cos(lat)*cos(lon - lon0)))
        x = false_easting + r*k_lambert*cos(lat)*sin(lon - lon0)
        y = false_northing + r*k_lambert*(cos(lat0)*sin(lat) - sin(lat0)*cos(lat)*cos(lon - lon0))
    end subroutine LL2LAEA_spherical