subroutine da_xyll_ps(x, y, proj, lat, lon) ! This is the inverse subroutine da_of llij_ps. It returns the ! latitude and longitude of an x/y point given the projection info ! structure. implicit none real, intent(in) :: x ! Column real, intent(in) :: y ! Row type (proj_info), intent(in) :: proj real, intent(out) :: lat ! -90 -> 90 North real, intent(out) :: lon ! -180 -> 180 East real :: reflon real :: scale_top real :: xx,yy real :: gi2, r2 real :: arccos if (trace_use_frequent) call da_trace_entry("da_xyll_ps") ! Compute the reference longitude by rotating 90 degrees to the east ! to find the longitude line parallel to the positive x-axis. reflon = proj%stdlon + 90.0 ! Compute numerator term of map scale factor scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg) ! Compute radius to point of interest xx = x - proj%polei yy = (y - proj%polej) * proj%hemi r2 = xx**2 + yy**2 ! Now the magic code if (r2 == 0.0) then lat = proj%hemi * 90.0 lon = reflon else gi2 = (proj%rebydx * scale_top)**2.0 lat = deg_per_rad * proj%hemi * ASin((gi2-r2)/(gi2+r2)) arccos = ACOS(xx/sqrt(r2)) if (yy > 0) then lon = reflon + deg_per_rad * arccos else lon = reflon - deg_per_rad * arccos end if end if ! Convert to a -180 -> 180 East convention if (lon > 180.0) lon = lon - 360.0 if (lon < -180.0) lon = lon + 360.0 if (trace_use_frequent) call da_trace_exit("da_xyll_ps") end subroutine da_xyll_ps