var/da/da_tools/da_llxy_ps.inc

   1 subroutine da_llxy_ps(lat,lon,proj,x,y)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Given latitude (-90 to 90), longitude (-180 to 180), and the
   5    ! standard polar-stereographic projection information via the
   6    ! public proj structure, this routine returns the x/y indices which
   7    ! if within the domain range from 1->nx and 1->ny, respectively.
   8    !-----------------------------------------------------------------------
   9
  10    implicit none
  11
  12    real, intent(in)               :: lat
  13    real, intent(in)               :: lon
  14    type(proj_info),intent(in)     :: proj
  15
  16    real, intent(out)              :: x !(x-index)
  17    real, intent(out)              :: y !(y-index)
  18
  19    real                           :: reflon
  20    real                           :: scale_top
  21    real                           :: ala
  22    real                           :: alo
  23    real                           :: rm
  24
  25    if (trace_use_frequent) call da_trace_entry("da_llxy_ps")
  26
  27    reflon = proj%stdlon + 90.0
  28
  29    ! Compute numerator term of map scale factor
  30
  31    scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
  32
  33    ! Find radius to desired point
  34    ala = lat * rad_per_deg
  35    rm = proj%rebydx * COS(ala) * scale_top/(1.0 + proj%hemi *Sin(ala))
  36    alo = (lon - reflon) * rad_per_deg
  37    x = proj%polei + rm * COS(alo)
  38    y = proj%polej + proj%hemi * rm * Sin(alo)
  39
  40    if (trace_use_frequent) call da_trace_exit("da_llxy_ps")
  41
  42 end subroutine da_llxy_ps
  43
  44