FUNCTION GET_GRID_INDEX, LAT, LON

;+
; NAME:
;    GET_GRID_INDEX
;
; PURPOSE:
;    For a given latitude and longitude, compute the grid index value
;    in a global equal area grid (25 km by 25 km).
;
;    The grid is defined as follows:
;    (1) Assume equatorial and polar circumference of Earth is 40000 km.
;    (2) Divide the latitude range from -90 to 90 into
;        800 bands spaced at 25 km (20000 km / 25 km = 800).
;    (3) Divide the longitude range at the equator from 0 to 360 into
;        1600 bands spaced at 25 km (40000 km / 25 km = 1600).
;    (4) For each latitude band north or south of the equator,
;        divide the longitude range from 0 to 360 into
;        (1600 * cosine(latitude)) bands.
;
;    Adapted from the UW MODIS shared/src function get_grid_index.f
;
; CATEGORY:
;    MODIS utilities.
;
; CALLING SEQUENCE:
;    RESULT = GET_GRID_INDEX(LAT, LON)
;
; INPUTS:
;    LAT               Array of latitudes (degrees, -90S to +90N).
;    LON               Array of longitudes (degrees, -180W to 180E, Greenwich=0).
;
; OPTIONAL INPUTS:
;    None.
;
; KEYWORD PARAMETERS:
;    None.
;
; OUTPUTS:
;    GET_GRID_INDEX    Array of index values in global grid (range 0-814879).
;
; OPTIONAL OUTPUTS:
;    None
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None.
;
; RESTRICTIONS:
;    Requires IDL 5.0 or higher (square bracket array syntax).
;
; EXAMPLE:
;
;print, get_grid_index(45.0, -90.0)
;
; MODIFICATION HISTORY:
; Liam.Gumley@ssec.wisc.edu
; http://cimss.ssec.wisc.edu/~gumley
; $Id: get_grid_index.pro,v 1.2 2000/01/07 19:05:11 gumley Exp $
;
; Copyright (C) 1999, 2000 Liam E. Gumley
;
; This program is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public License
; as published by the Free Software Foundation; either version 2
; of the License, or (at your option) any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
;-

rcs_id = '$Id: get_grid_index.pro,v 1.2 2000/01/07 19:05:11 gumley Exp $'

;- Check arguments
if (n_params() ne 2) then message, 'Usage: RESULT = GET_GRID_INDEX(LAT, LON)'
if (n_elements(lat) eq 0) then message, 'Argument LAT is undefined'
if (n_elements(lon) eq 0) then message, 'Argument LON is undefined'

;- Define common block
common get_grid_index_data, num_boxes, lon_inc

;- Read data file if common block elements are undefined
if (n_elements(num_boxes) ne 800) then begin
  message, 'Reading get_grid_index data', /continue
  openr, lun, 'get_grid_index.dat', /get_lun
  header = ''
  readf, lun, header, format='(a)'
  data = fltarr(2, 800)
  readf, lun, data
  free_lun, lun
  num_boxes = long(reform(data[0, *]))
  lon_inc = reform(data[1, *])
endif

;- Set constants for 25 km equal area grid
grdsiz     = 25.0
ce         = 40000.0
beg_lat    = 90.0
beg_lon    = 0.0
num_bands  = round((ce / grdsiz) * 0.5)
lat_inc    = 180.0 / num_bands
beg_ltindx = long((90.0 - beg_lat) / lat_inc)

;- Check input arguments
check = where((lat lt -90.0) or (lat gt 90.0), count)
if (count gt 0) then message, 'LAT values must be in the range [-90, 90]'
check = where((lon lt -180.0) or (lon gt 180.0), count)
if (count gt 0) then message, 'LON values must be in the range [-180, 180]'
      
;- Find latitude band
lat_band = (long((90.0 - lat) / lat_inc) + 1L - beg_ltindx) < 800L

;- Convert longitude (west negative, -180 to 180) to
;- McIDAS longitude (west positive, -180 to 180)
alon = -lon
      
;- Convert McIDAS longitude (west positive, -180 to 180) to
;- Greenwich system (Greenwich zero, 0 to 360)
flags = (alon le 0.0)
index = where(flags eq 1, count)
if (count gt 0) then alon[index] = abs(alon[index])
index = where(flags eq 0, count)
if (count gt 0) then alon[index] = 360.0 - alon[index]
       
;- Get longitude bin of the current latitude band
beg_lnindx = long(beg_lon / lon_inc[lat_band - 1L])
lon_bin = long( alon / lon_inc[lat_band - 1L]) + 1 - beg_lnindx

npboxes = num_boxes[(lat_band - 1 - 1) > 0L]
index = where(lat_band eq 1L, count)
if (count gt 0) then npboxes[index] = 0L

;- Compute grid index
return, npboxes + lon_bin - 1L
     
END