geogrid/util/latlon_to_ij.F

   1 subroutine lltoxy(lat,lon,rx,ry,stagger)
   2
   3   implicit none
   4
   5   ! Arguments
   6   integer, intent(in) :: stagger
   7   real, intent(in) :: lat, lon
   8   real, intent(out) :: rx, ry
   9
  10   ! Local variables
  11   real :: dphd,dlmd !Grid increments, degrees
  12   integer :: ii,imt,jj,jmt,k,krows,ncol,nrow
  13   real :: glatd  !Geographic latitude, positive north
  14   real :: glond  !Geographic longitude, positive west
  15   real :: col,d1,d2,d2r,dlm,dlm1,dlm2,dph,glat,glon,    &
  16           pi,r2d,row,tlat,tlat1,tlat2,              &
  17           tlon,tlon1,tlon2,tph0,tlm0,x,y,z
  18
  19     glatd = lat
  20     glond = lon
  21
  22     dphd = phi/real((jydim(current_nest_number)-1)/2)
  23     dlmd = lambda/real(ixdim(current_nest_number)-1)
  24
  25     pi = dacos(-1.0)
  26     d2r = pi/180.
  27     r2d = 1./d2r
  28
  29     imt = 2*ixdim(current_nest_number)-1
  30     jmt = jydim(current_nest_number)/2+1
  31
  32     glat = glatd*d2r
  33     glon = glond*d2r
  34     dph = dphd*d2r
  35     dlm = dlmd*d2r
  36     tph0 = known_lat*d2r
  37     tlm0 = known_lon*d2r
  38
  39     x = cos(tph0)*cos(glat)*cos(glon-tlm0)+sin(tph0)*sin(glat)
  40     y = -cos(glat)*sin(glon-tlm0)
  41     z = cos(tph0)*sin(glat)-sin(tph0)*cos(glat)*cos(glon-tlm0)
  42     tlat = r2d*atan(z/sqrt(x*x+y*y))
  43     tlon = r2d*atan(y/x)
  44
  45     row = tlat/dphd+jmt
  46     col = tlon/dlmd+ixdim(current_nest_number)
  47     nrow = int(row)
  48     ncol = int(col)
  49     tlat = tlat*d2r
  50     tlon = tlon*d2r
  51
  52     if (stagger == HH) then
  53
  54       if ((mod(nrow,2) == 1 .and. mod(ncol,2) == 1) .or. &
  55           (mod(nrow,2) == 0 .and. mod(ncol,2) == 0)) then
  56         tlat1 = (nrow-jmt)*dph
  57         tlat2 = tlat1+dph
  58         tlon1 = (ncol-ixdim(current_nest_number))*dlm
  59         tlon2 = tlon1+dlm
  60         dlm1 = tlon-tlon1
  61         dlm2 = tlon-tlon2
  62         d1 = acos(cos(tlat)*cos(tlat1)*cos(dlm1)+sin(tlat)*sin(tlat1))
  63         d2 = acos(cos(tlat)*cos(tlat2)*cos(dlm2)+sin(tlat)*sin(tlat2))
  64
  65         if (d1 > d2) then
  66           nrow = nrow+1
  67           ncol = ncol+1
  68         end if
  69
  70       else
  71         tlat1 = (nrow+1-jmt)*dph
  72         tlat2 = tlat1-dph
  73         tlon1 = (ncol-ixdim(current_nest_number))*dlm
  74         tlon2 = tlon1+dlm
  75         dlm1 = tlon-tlon1
  76         dlm2 = tlon-tlon2
  77         d1 = acos(cos(tlat)*cos(tlat1)*cos(dlm1)+sin(tlat)*sin(tlat1))
  78         d2 = acos(cos(tlat)*cos(tlat2)*cos(dlm2)+sin(tlat)*sin(tlat2))
  79
  80         if (d1 < d2) then
  81           nrow = nrow+1
  82         else
  83           ncol = ncol+1
  84         end if
  85       end if
  86
  87     else if (stagger == VV) then
  88
  89       if ((mod(nrow,2) == 0 .and. mod(ncol,2) == 1) .or. &
  90           (mod(nrow,2) == 1 .and. mod(ncol,2) == 0)) then
  91         tlat1 = (nrow-jmt)*dph
  92         tlat2 = tlat1+dph
  93         tlon1 = (ncol-ixdim(current_nest_number))*dlm
  94         tlon2 = tlon1+dlm
  95         dlm1 = tlon-tlon1
  96         dlm2 = tlon-tlon2
  97         d1 = acos(cos(tlat)*cos(tlat1)*cos(dlm1)+sin(tlat)*sin(tlat1))
  98         d2 = acos(cos(tlat)*cos(tlat2)*cos(dlm2)+sin(tlat)*sin(tlat2))
  99
 100         if (d1 > d2) then
 101           nrow = nrow+1
 102           ncol = ncol+1
 103         end if
 104
 105       else
 106         tlat1 = (nrow+1-jmt)*dph
 107         tlat2 = tlat1-dph
 108         tlon1 = (ncol-ixdim(current_nest_number))*dlm
 109         tlon2 = tlon1+dlm
 110         dlm1 = tlon-tlon1
 111         dlm2 = tlon-tlon2
 112         d1 = acos(cos(tlat)*cos(tlat1)*cos(dlm1)+sin(tlat)*sin(tlat1))
 113         d2 = acos(cos(tlat)*cos(tlat2)*cos(dlm2)+sin(tlat)*sin(tlat2))
 114
 115         if (d1 < d2) then
 116           nrow = nrow+1
 117         else
 118           ncol = ncol+1
 119         end if
 120       end if
 121     end if
 122
 123     jj = nrow
 124     ii = ncol/2
 125     if (stagger == HH) then
 126       if (mod(jj,2) == 1) ii = ii+1
 127       krows = ((nrow-1)/2)*imt
 128       if (mod(nrow,2) == 1) then
 129         k = krows+(ncol+1)/2
 130       else
 131         k = krows+ixdim(current_nest_number)+ncol/2
 132       end if
 133
 134     else if (stagger == VV) then
 135       if (mod(jj,2) == 0) ii=ii+1
 136
 137       krows = ((nrow-1)/2)*imt
 138       if (mod(nrow,2) == 1) then
 139         k = krows+ncol/2
 140       else
 141         k = krows+ixdim(current_nest_number)-1+(ncol+1)/2
 142       end if
 143     end if
 144
 145     rx = real(ii)
 146     ry = real(jj)
 147
 148 end subroutine lltoxy