var/da/da_tools/da_llxy_rotated_latlon.inc

   1 subroutine da_llxy_rotated_latlon(lat,lon, proj, x, y)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Compute the x/y location of a lat/lon on a rotated LATLON grid.
   5    ! Author :  Syed RH Rizvi,     MMM/NCAR
   6    !           06/01/2008
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    real, intent(in)             :: lat
  12    real, intent(in)             :: lon
  13    type(proj_info), intent(in)  :: proj
  14    real, intent(out)            :: x
  15    real, intent(out)            :: y
  16
  17    real                         :: rot_lat, rot_lon, deltalat,deltalon, lon360,latinc,loninc
  18    real                         :: xlat, xlon, cen_lat, cen_lon
  19
  20
  21    if (trace_use_frequent) call da_trace_entry("da_llxy_rotated_latlon")
  22    ! To account for issues around the dateline, convert the incoming
  23    ! longitudes to be 0->360.0
  24    if (lon < 0) then
  25       lon360 = lon + 360.0
  26    else
  27       lon360 = lon
  28    end if
  29
  30    xlat = deg_to_rad*lat
  31    xlon = deg_to_rad*lon360
  32    cen_lat = deg_to_rad*proj%lat1
  33    cen_lon = deg_to_rad*proj%lon1
  34    if (cen_lon < 0.) cen_lon = cen_lon + 360.
  35
  36    latinc = proj%latinc
  37    loninc = proj%loninc
  38
  39    rot_lon = rad_to_deg*atan( cos(xlat) * sin(xlon-cen_lon)/ &
  40              (cos(cen_lat)*cos(xlat)*cos(xlon-cen_lon) + sin(cen_lat)*sin(xlat)))
  41    rot_lat = rad_to_deg*asin(  cos(cen_lat)*sin(xlat) - sin(cen_lat)*cos(xlat)*cos(xlon-cen_lon))
  42
  43
  44    deltalat = rot_lat
  45    deltalon = rot_lon
  46
  47
  48    ! Compute x/y
  49    x = proj%knowni + deltalon/loninc + 1.0
  50    y = proj%knownj + deltalat/latinc + 1.0
  51
  52    if (trace_use_frequent) call da_trace_exit("da_llxy_rotated_latlon")
  53
  54 end subroutine da_llxy_rotated_latlon