WPS/geogrid/util/ij_to_latlon.F

   1 subroutine xytoll(rx,ry,lat,lon,stagger)
   2
   3   implicit none
   4
   5   ! Arguments
   6   integer, intent(in) :: stagger
   7   real, intent(in) :: rx, ry
   8   real, intent(out) :: lat, lon
   9
  10   ! Local variables
  11   integer :: ih,jh
  12   integer :: midcol,midrow,ncol,iadd1,iadd2,imt,jh2,knrow,krem,kv,nrow
  13   real :: dphd,dlmd !Grid increments, degrees
  14   real :: arg1,arg2,d2r,fctr,glatr,glatd,glond,pi, &
  15           r2d,tlatd,tlond,tlatr,tlonr,tlm0,tph0
  16
  17     ih = nint(rx)
  18     jh = nint(ry)
  19
  20     dphd = phi/real((jydim(current_nest_number)-1)/2)
  21     dlmd = lambda/real(ixdim(current_nest_number)-1)
  22
  23     pi = dacos(-1.0)
  24     d2r = pi/180.
  25     r2d = 1./d2r
  26     tph0 = known_lat*d2r
  27     tlm0 = known_lon*d2r
  28
  29     midrow = (jydim(current_nest_number)+1)/2
  30     midcol = ixdim(current_nest_number)
  31
  32     if (stagger == HH) then
  33       ncol = 2*ih-1+mod(jh+1,2)
  34       tlatd = (jh-midrow)*dphd
  35       tlond = (ncol-midcol)*dlmd
  36     else if (stagger == VV) then
  37       imt = 2*ixdim(current_nest_number)-1
  38       jh2 = jh/2
  39       iadd1 = 0
  40       iadd2 = 0
  41
  42       if (2*jh2 == jh) then
  43         iadd1 = -1
  44         iadd2 = ixdim(current_nest_number)-1
  45       end if
  46
  47       kv = (jh2+iadd1)*imt+iadd2+ih
  48
  49       nrow = 2*((kv-1)/imt)
  50       knrow = imt*nrow/2
  51       krem = kv-knrow
  52
  53       if (krem <= ixdim(current_nest_number)-1) then
  54         nrow = nrow+1
  55         ncol = 2*krem
  56       else
  57         nrow = nrow+2
  58         ncol = 2*(krem-ixdim(current_nest_number))+1
  59       end if
  60       tlatd = (nrow-(jydim(current_nest_number)+1)/2)*dphd
  61       tlond = (ncol-ixdim(current_nest_number))*dlmd
  62     end if
  63
  64     tlatr = tlatd*d2r
  65     tlonr = tlond*d2r
  66     arg1 = sin(tlatr)*cos(tph0)+cos(tlatr)*sin(tph0)*cos(tlonr)
  67     glatr = asin(arg1)
  68
  69     glatd = glatr*r2d
  70
  71     arg2 = dcos(tlatr)*dcos(tlonr)/(dcos(glatr)*dcos(tph0))-dtan(glatr)*dtan(tph0)
  72     if (abs(arg2) > 1.) arg2 = abs(arg2)/arg2
  73     fctr = 1.
  74     if (tlond > 0.) fctr = -1.
  75
  76     glond = known_lon+fctr*dacos(arg2)*r2d
  77
  78     xlat = glatd
  79     xlon = glond
  80
  81 end subroutine xytoll