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