da_llij_ps.inc

1 subroutine da_llij_ps(lat,lon,proj,i,j)
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 i/j 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)              :: i !(x-index)
17    real, intent(out)              :: j !(y-index)
18    
19    real                           :: reflon
20    real                           :: scale_top
21    real                           :: ala
22    real                           :: alo
23    real                           :: rm
24 
25    reflon = proj%stdlon + 90.
26    
27    ! Compute numerator term of map scale factor
28 
29    scale_top = 1. + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
30 
31    ! Find radius to desired point
32    ala = lat * rad_per_deg
33    rm = proj%rebydx * COS(ala) * scale_top/(1. + proj%hemi *Sin(ala))
34    alo = (lon - reflon) * rad_per_deg
35    i = proj%polei + rm * COS(alo)
36    j = proj%polej + proj%hemi * rm * Sin(alo)
37  
38 end subroutine da_llij_ps
39 
40