da_ijll_lc.inc

1 subroutine da_ijll_lc(i, j, proj, lat, lon)
2 
3    ! subroutine da_to convert from the (i,j) cartesian coordinate to the 
4    ! geographical latitude and longitude for a Lambert Conformal projection.
5 
6    implicit none
7 
8    real, intent(in)              :: i        ! Cartesian X coordinate
9    real, intent(in)              :: j        ! 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    chi1 = (90. - proj%hemi*proj%truelat1)*rad_per_deg
25    chi2 = (90. - proj%hemi*proj%truelat2)*rad_per_deg
26 
27    ! See if we are in the southern hemispere and flip the indices
28    ! if we are. 
29    if (proj%hemi == -1.0) then 
30       inew = -i + 2.0
31       jnew = -j + 2.0
32    else
33       inew = i
34       jnew = j
35    end if
36 
37    ! Compute radius**2 to i/j location
38    xx = inew - proj%polei
39    yy = proj%polej - jnew
40    r2 = (xx*xx + yy*yy)
41    r = sqrt(r2)/proj%rebydx
42    
43    ! Convert to lat/lon
44    if (r2 == 0.0) then
45       lat = proj%hemi * 90.0
46       lon = proj%stdlon
47    else
48        
49       ! Longitude
50       lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
51       lon = MOD(lon+360., 360.)
52 
53       ! Latitude.  Latitude determined by solving an equation adapted 
54       ! from:
55       !  Maling, D.H., 1973: Coordinate Systems and Map Projections
56       ! Equations #20 in Appendix I.  
57         
58       if (chi1 == chi2) then
59          chi = 2.0*ATAN((r/TAN(chi1))**(1./proj%cone) * TAN(chi1*0.5))
60       else
61          chi = 2.0*ATAN((r*proj%cone/Sin(chi1))**(1./proj%cone) * TAN(chi1*0.5)) 
62       end if
63       lat = (90.0-chi*deg_per_rad)*proj%hemi
64    end if
65 
66    if (lon > +180.) lon = lon - 360.
67    if (lon < -180.) lon = lon + 360.
68 
69 end subroutine da_ijll_lc
70 
71