da_ijll_ps.inc

1 subroutine da_ijll_ps(i, j, proj, lat, lon)
2 
3    ! This is the inverse subroutine da_of llij_ps.  It returns the 
4    ! latitude and longitude of an i/j point given the projection info 
5    ! structure.  
6 
7    implicit none
8 
9    real, intent(in)                    :: i    ! Column
10    real, intent(in)                    :: j    ! Row
11    type (proj_info), intent(in)        :: proj
12     
13    real, intent(out)                   :: lat     ! -90 -> 90 North
14    real, intent(out)                   :: lon     ! -180 -> 180 East
15 
16    real                                :: reflon
17    real                                :: scale_top
18    real                                :: xx,yy
19    real                                :: gi2, r2
20    real                                :: arccos
21 
22    ! Compute the reference longitude by rotating 90 degrees to the east
23    ! to find the longitude line parallel to the positive x-axis.
24    reflon = proj%stdlon + 90.
25    
26    ! Compute numerator term of map scale factor
27    scale_top = 1. + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
28 
29    ! Compute radius to point of interest
30    xx = i - proj%polei
31    yy = (j - proj%polej) * proj%hemi
32    r2 = xx**2 + yy**2
33 
34    ! Now the magic code
35    if (r2 == 0.) then 
36       lat = proj%hemi * 90.
37       lon = reflon
38    else
39       gi2 = (proj%rebydx * scale_top)**2.
40       lat = deg_per_rad * proj%hemi * ASin((gi2-r2)/(gi2+r2))
41       arccos = ACOS(xx/sqrt(r2))
42       if (yy > 0) then
43          lon = reflon + deg_per_rad * arccos
44       else
45          lon = reflon - deg_per_rad * arccos
46       end if
47    end if
48   
49    ! Convert to a -180 -> 180 East convention
50    if (lon > 180.) lon = lon - 360.
51    if (lon < -180.) lon = lon + 360.
52    
53   
54 end subroutine da_ijll_ps
55 
56