da_xyll_lc.inc

1 subroutine da_xyll_lc(x, y, proj, lat, lon)
2 
3    ! subroutine da_to convert from the (x,y) cartesian coordinate to the 
4    ! geographical latitude and longitude for a Lambert Conformal projection.
5 
6    implicit none
7 
8    real, intent(in)              :: x        ! Cartesian X coordinate
9    real, intent(in)              :: y        ! Cartesian Y coordinate
10    type(proj_info),intent(in)    :: proj     ! Projection info structure
11 
12                 
13    real, intent(out)             :: lat      ! Latitude (-90->90 deg N)
14    real, intent(out)             :: lon      ! Longitude (-180->180 E)
15 
16    real                          :: inew
17    real                          :: jnew
18    real                          :: r
19    real                          :: chi,chi1,chi2
20    real                          :: r2
21    real                          :: xx
22    real                          :: yy
23 
24    if (trace_use_dull) call da_trace_entry("da_xyll_lc")
25 
26 
27    chi1 = (90.0 - proj%hemi*proj%truelat1)*rad_per_deg
28    chi2 = (90.0 - proj%hemi*proj%truelat2)*rad_per_deg
29 
30    ! See if we are in the southern hemispere and flip the indices
31    ! if we are. 
32    if (proj%hemi == -1.0) then 
33       inew = -x + 2.0
34       jnew = -y + 2.0
35    else
36       inew = x
37       jnew = y
38    end if
39 
40    ! Compute radius**2 to i/j location
41    xx = inew - proj%polei
42    yy = proj%polej - jnew
43    r2 = (xx*xx + yy*yy)
44    r = sqrt(r2)/proj%rebydx
45    
46    ! Convert to lat/lon
47    if (r2 == 0.0) then
48       lat = proj%hemi * 90.0
49       lon = proj%stdlon
50    else
51        
52       ! Longitude
53       lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
54       lon = MOD(lon+360.0, 360.0)
55 
56       ! Latitude.  Latitude determined by solving an equation adapted 
57       ! from:
58       !  Maling, D.H., 1973: Coordinate Systems and Map Projections
59       ! Equations #20 in Appendix I.  
60         
61       if (chi1 == chi2) then
62          chi = 2.0*ATAN((r/TAN(chi1))**(1.0/proj%cone) * TAN(chi1*0.5))
63       else
64          chi = 2.0*ATAN((r*proj%cone/Sin(chi1))**(1.0/proj%cone) * TAN(chi1*0.5)) 
65       end if
66       lat = (90.0-chi*deg_per_rad)*proj%hemi
67    end if
68 
69    if (lon > +180.0) lon = lon - 360.0
70    if (lon < -180.0) lon = lon + 360.0
71 
72    if (trace_use_dull) call da_trace_exit("da_xyll_lc")
73 
74 end subroutine da_xyll_lc
75 
76