<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_LLXY_LC_NEW'><A href='../../html_code/tools/da_llxy_lc_new.inc.html#DA_LLXY_LC_NEW' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_llxy_lc_new(proj, info) 1,2
!-----------------------------------------------------------------------
! 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