subroutine da_llxy_lc_new(proj, info)

   !-----------------------------------------------------------------------
   ! Purpose: compute the geographical latitude and longitude values
   ! to the cartesian x/y on a Lambert Conformal projection.
   !-----------------------------------------------------------------------
    
   implicit none

   type(proj_info), intent(in)    :: proj     ! Projection info structure
   type(infa_type), intent(inout) :: info


   real    :: tl1r
!   real    :: temp1
!   real    :: temp2
   real    :: ctl1r
   real    :: deltalon, rm, arg
   integer :: n

   if (trace_use) call da_trace_entry("da_llxy_lc_new")

! FAST
    
   ! Convert truelat1 to radian and compute COS for later use
!   tl1r = proj%truelat1 * rad_per_deg
!   ctl1r = COS(tl1r)    
!   temp1 = TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0)
!   temp2 = proj%rebydx * ctl1r/proj%cone
    
   ! Compute deltalon between known longitude and standard lon and ensure
   ! it is not in the cut zone

!   where (lon - proj%stdlon > +180.0)
!      x = proj%polei + proj%hemi * (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
!         * SIN(proj%cone*((lon - proj%stdlon-360.0)*rad_per_deg))
!      y = proj%polej - (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
!         * COS(proj%cone*((lon - proj%stdlon-360.0)*rad_per_deg))
!   elsewhere  (lon - proj%stdlon - -180.0)
!      x = proj%polei + proj%hemi * (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1) &
!         * SIN(proj%cone*((lon - proj%stdlon+360.0)*rad_per_deg))
!      y = proj%polej - (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
!         * COS(proj%cone*((lon - proj%stdlon+360.0)*rad_per_deg))
!   else
!      x = proj%polei + proj%hemi * (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1) &
!         * SIN(proj%cone*((lon - proj%stdlon)*rad_per_deg))
!      y = proj%polej - (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
!         * COS(proj%cone*((lon - proj%stdlon)*rad_per_deg))
!   end where

   ! Finally, if we are in the southern hemisphere, flip the i/j
   ! values to a coordinate system where (1,1) is the SW corner
   ! (what we assume) which is different than the original NCEP
   ! algorithms which used the NE corner as the origin in the 
   ! southern hemisphere (left-hand vs. right-hand coordinate?)
!   if (proj%hemi == -1.0) then
!      x(:,:) = 2.0 - i(:,:)
!      y(:,:) = 2.0 - j(:,:)
!   end if

! SLOW

   do n=lbound(info%lat,2), ubound(info%lat,2)
      ! Compute deltalon between known longitude and standard lon and ensure
      ! it is not in the cut zone
      deltalon = info%lon(1,n) - proj%stdlon
      if (deltalon > +180.0) deltalon = deltalon - 360.0
      if (deltalon < -180.0) deltalon = deltalon + 360.0

      ! Convert truelat1 to radian and compute COS for later use
      tl1r = proj%truelat1 * rad_per_deg
      ctl1r = COS(tl1r)     

      ! Radius to desired point
      rm = proj%rebydx * ctl1r/proj%cone * &
          (TAN((90.0*proj%hemi-info%lat(1,n))*rad_per_deg/2.0) / &
           TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone

      arg = proj%cone*(deltalon*rad_per_deg)
      info%x(:,n) = proj%polei + proj%hemi * rm * Sin(arg)
      info%y(:,n) = proj%polej - rm * COS(arg)

      ! Finally, if we are in the southern hemisphere, flip the i/j
      ! values to a coordinate system where (1,1) is the SW corner
      ! (what we assume) which is different than the original NCEP
      ! algorithms which used the NE corner as the origin in the 
      ! southern hemisphere (left-hand vs. right-hand coordinate?)
      if (proj%hemi == -1.0) then
         info%x(:,n) = 2.0 - info%x(:,n)  
	 info%y(:,n) = 2.0 - info%y(:,n)
      end if
   end do

   if (trace_use) call da_trace_exit("da_llxy_lc_new")

end subroutine da_llxy_lc_new