da_llxy_ps.inc

1 subroutine da_llxy_ps(lat,lon,proj,x,y)
2 
3    !-----------------------------------------------------------------------
4    ! Purpose: Given latitude (-90 to 90), longitude (-180 to 180), and the
5    ! standard polar-stereographic projection information via the 
6    ! public proj structure, this routine returns the x/y indices which
7    ! if within the domain range from 1->nx and 1->ny, respectively.
8    !-----------------------------------------------------------------------
9 
10    implicit none
11 
12    real, intent(in)               :: lat
13    real, intent(in)               :: lon
14    type(proj_info),intent(in)     :: proj
15 
16    real, intent(out)              :: x !(x-index)
17    real, intent(out)              :: y !(y-index)
18    
19    real                           :: reflon
20    real                           :: scale_top
21    real                           :: ala
22    real                           :: alo
23    real                           :: rm
24 
25    if (trace_use_frequent) call da_trace_entry("da_llxy_ps")
26 
27    reflon = proj%stdlon + 90.0
28    
29    ! Compute numerator term of map scale factor
30 
31    scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
32 
33    ! Find radius to desired point
34    ala = lat * rad_per_deg
35    rm = proj%rebydx * COS(ala) * scale_top/(1.0 + proj%hemi *Sin(ala))
36    alo = (lon - reflon) * rad_per_deg
37    x = proj%polei + rm * COS(alo)
38    y = proj%polej + proj%hemi * rm * Sin(alo)
39 
40    if (trace_use_frequent) call da_trace_exit("da_llxy_ps")
41  
42 end subroutine da_llxy_ps
43 
44