subroutine da_llxy_ps_new(proj,info)

   !-----------------------------------------------------------------------
   ! Purpose: Given latitude (-90 to 90), longitude (-180 to 180), and the
   ! standard polar-stereographic projection information via the 
   ! public proj structure, this routine returns the x/y indices which
   ! if within the domain range from 1->nx and 1->ny, respectively.
   !-----------------------------------------------------------------------

   implicit none

   type(proj_info), intent(in)    :: proj
   type(infa_type), intent(inout) :: info
   
   real    :: reflon
   integer :: n
   real    :: scale_top, ala, rm, alo

   if (trace_use) call da_trace_entry("da_llxy_ps_new")

   reflon = proj%stdlon + 90.0

   ! FAST

!   x(:,:) = proj%polei + proj%rebydx * COS(lat(:,:) * rad_per_deg) &
!      * (1.0 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) &
!      / (1.0 + proj%hemi *SIN(lat(:,:) * rad_per_deg)) &
!      * COS((lon(:,:) - reflon) * rad_per_deg)

!   y(:,:) = proj%polej + proj%hemi * proj%rebydx * COS(lat(:,:) * rad_per_deg) &
!      * (1.0 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) &
!      / (1.0 + proj%hemi *SIN(lat(:,:) * rad_per_deg)) &
!      * SIN((lon(:,:) - reflon) * rad_per_deg)

! SLOW

   do n=lbound(info%lat,2),ubound(info%lat,2)
      ! compute numerator term of map scale factor

      scale_top = 1.0 + proj%hemi * sin(proj%truelat1 * rad_per_deg)

      ! find radius to desired point
      ala = info%lat(1,n) * rad_per_deg
      rm = proj%rebydx * cos(ala) * scale_top/(1.0 + proj%hemi *sin(ala))
      alo = (info%lon(1,n) - reflon) * rad_per_deg
      info%x(:,n) = proj%polei + rm * cos(alo)
      info%y(:,n) = proj%polej + proj%hemi * rm * sin(alo)
   end do

   if (trace_use) call da_trace_exit("da_llxy_ps_new")
 
end subroutine da_llxy_ps_new