da_llij_lc.inc

1 subroutine da_llij_lc(lat, lon, proj, i, j)
2 
3    !-----------------------------------------------------------------------
4    ! Purpose: compute the geographical latitude and longitude values
5    ! to the cartesian x/y on a Lambert Conformal projection.
6    !-----------------------------------------------------------------------
7     
8    implicit none
9 
10    real, intent(in)              :: lat      ! Latitude (-90->90 deg N)
11    real, intent(in)              :: lon      ! Longitude (-180->180 E)
12    type(proj_info),intent(in)      :: proj     ! Projection info structure
13 
14    real, intent(out)             :: i        ! Cartesian X coordinate
15    real, intent(out)             :: j        ! Cartesian Y coordinate
16 
17    real                          :: arg
18    real                          :: deltalon
19    real                          :: tl1r
20    real                          :: rm
21    real                          :: ctl1r
22     
23    ! Compute deltalon between known longitude and standard lon and ensure
24    ! it is not in the cut zone
25    deltalon = lon - proj%stdlon
26    if (deltalon > +180.0) deltalon = deltalon - 360.0
27    if (deltalon < -180.0) deltalon = deltalon + 360.0
28     
29    ! Convert truelat1 to radian and compute COS for later use
30    tl1r = proj%truelat1 * rad_per_deg
31    ctl1r = COS(tl1r)     
32    
33    ! Radius to desired point
34    rm = proj%rebydx * ctl1r/proj%cone * &
35        (TAN((90.*proj%hemi-lat)*rad_per_deg/2.) / &
36         TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone
37 
38    arg = proj%cone*(deltalon*rad_per_deg)
39    i = proj%polei + proj%hemi * rm * Sin(arg)
40    j = proj%polej - rm * COS(arg)
41 
42    ! Finally, if we are in the southern hemisphere, flip the i/j
43    ! values to a coordinate system where (1,1) is the SW corner
44    ! (what we assume) which is different than the original NCEP
45    ! algorithms which used the NE corner as the origin in the 
46    ! southern hemisphere (left-hand vs. right-hand coordinate?)
47    if (proj%hemi == -1.0) then
48       i = 2.0 - i  
49       j = 2.0 - j
50    end if
51 
52 end subroutine da_llij_lc
53 
54