da_xyll_ps.inc

1 subroutine da_xyll_ps(x, y, proj, lat, lon)
2 
3    ! This is the inverse subroutine da_of llij_ps.  It returns the 
4    ! latitude and longitude of an x/y point given the projection info 
5    ! structure.  
6 
7    implicit none
8 
9    real, intent(in)                    :: x    ! Column
10    real, intent(in)                    :: y    ! 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    if (trace_use_frequent) call da_trace_entry("da_xyll_ps")
23 
24    ! Compute the reference longitude by rotating 90 degrees to the east
25    ! to find the longitude line parallel to the positive x-axis.
26    reflon = proj%stdlon + 90.0
27    
28    ! Compute numerator term of map scale factor
29    scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
30 
31    ! Compute radius to point of interest
32    xx = x - proj%polei
33    yy = (y - proj%polej) * proj%hemi
34    r2 = xx**2 + yy**2
35 
36    ! Now the magic code
37    if (r2 == 0.0) then 
38       lat = proj%hemi * 90.0
39       lon = reflon
40    else
41       gi2 = (proj%rebydx * scale_top)**2.0
42       lat = deg_per_rad * proj%hemi * ASin((gi2-r2)/(gi2+r2))
43       arccos = ACOS(xx/sqrt(r2))
44       if (yy > 0) then
45          lon = reflon + deg_per_rad * arccos
46       else
47          lon = reflon - deg_per_rad * arccos
48       end if
49    end if
50   
51    ! Convert to a -180 -> 180 East convention
52    if (lon > 180.0) lon = lon - 360.0
53    if (lon < -180.0) lon = lon + 360.0
54 
55    if (trace_use_frequent) call da_trace_exit("da_xyll_ps")
56      
57 end subroutine da_xyll_ps
58 
59