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