var/da/da_tools/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