var/da/da_tools/da_llxy_ps_new.inc

   1 subroutine da_llxy_ps_new(proj,info)
   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    type(proj_info), intent(in)    :: proj
  13    type(infa_type), intent(inout) :: info
  14
  15    real    :: reflon
  16    integer :: n
  17    real    :: scale_top, ala, rm, alo
  18
  19    if (trace_use) call da_trace_entry("da_llxy_ps_new")
  20
  21    reflon = proj%stdlon + 90.0
  22
  23    ! FAST
  24
  25 !   x(:,:) = proj%polei + proj%rebydx * COS(lat(:,:) * rad_per_deg) &
  26 !      * (1.0 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) &
  27 !      / (1.0 + proj%hemi *SIN(lat(:,:) * rad_per_deg)) &
  28 !      * COS((lon(:,:) - reflon) * rad_per_deg)
  29
  30 !   y(:,:) = proj%polej + proj%hemi * proj%rebydx * COS(lat(:,:) * rad_per_deg) &
  31 !      * (1.0 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) &
  32 !      / (1.0 + proj%hemi *SIN(lat(:,:) * rad_per_deg)) &
  33 !      * SIN((lon(:,:) - reflon) * rad_per_deg)
  34
  35 ! SLOW
  36
  37    do n=lbound(info%lat,2),ubound(info%lat,2)
  38       ! compute numerator term of map scale factor
  39
  40       scale_top = 1.0 + proj%hemi * sin(proj%truelat1 * rad_per_deg)
  41
  42       ! find radius to desired point
  43       ala = info%lat(1,n) * rad_per_deg
  44       rm = proj%rebydx * cos(ala) * scale_top/(1.0 + proj%hemi *sin(ala))
  45       alo = (info%lon(1,n) - reflon) * rad_per_deg
  46       info%x(:,n) = proj%polei + rm * cos(alo)
  47       info%y(:,n) = proj%polej + proj%hemi * rm * sin(alo)
  48    end do
  49
  50    if (trace_use) call da_trace_exit("da_llxy_ps_new")
  51
  52 end subroutine da_llxy_ps_new
  53
  54