var/da/da_tools/da_llxy_lc.inc

   1 subroutine da_llxy_lc(lat, lon, proj, x, y)
   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)             :: x        ! Cartesian X coordinate
  15    real, intent(out)             :: y        ! Cartesian Y coordinate
  16
  17    real                          :: arg
  18    real                          :: deltalon
  19    real                          :: tl1r
  20    real                          :: rm
  21    real                          :: ctl1r
  22
  23    if (trace_use_dull) call da_trace_entry("da_llxy_lc")
  24
  25    ! Compute deltalon between known longitude and standard lon and ensure
  26    ! it is not in the cut zone
  27    deltalon = lon - proj%stdlon
  28    if (deltalon > +180.0) deltalon = deltalon - 360.0
  29    if (deltalon < -180.0) deltalon = deltalon + 360.0
  30
  31    ! Convert truelat1 to radian and compute COS for later use
  32    tl1r = proj%truelat1 * rad_per_deg
  33    ctl1r = COS(tl1r)
  34
  35    ! Radius to desired point
  36    rm = proj%rebydx * ctl1r/proj%cone * &
  37        (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / &
  38         TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone
  39
  40    arg = proj%cone*(deltalon*rad_per_deg)
  41    x = proj%polei + proj%hemi * rm * Sin(arg)
  42    y = proj%polej - rm * COS(arg)
  43
  44    ! Finally, if we are in the southern hemisphere, flip the i/j
  45    ! values to a coordinate system where (1,1) is the SW corner
  46    ! (what we assume) which is different than the original NCEP
  47    ! algorithms which used the NE corner as the origin in the
  48    ! southern hemisphere (left-hand vs. right-hand coordinate?)
  49    if (proj%hemi == -1.0) then
  50       x = 2.0 - x
  51       y = 2.0 - y
  52    end if
  53
  54    if (trace_use_dull) call da_trace_exit("da_llxy_lc")
  55
  56 end subroutine da_llxy_lc
  57
  58