function [data,grd]=cf_subsetGrid(uri,var,lonLatRange,dn1,dn2)
% CF_subsetGrid: Subset grid based on lat-lon bounding box and time
% [data, grd]=cf_subsetGrid(uri,var,[lonLatRange], dn1, dn2);
%
% Inputs: uri= CF-compliant NetCDF, OpenDAP or NCML file
%         var = variable to subset
%         lonLatRange = [minLon maxLon minLat maxLat]  % matlab 'axes' function      
%         dn1 = Matlab datenum, datestr or datevec  ex: [1990 4 5 0 0 0] or '5-Apr-1990 00:00'        
%         dn2 = Matlab datenum, datestr or datevec (optional)  

% Outputs: data = data % subset data  based on lonlat and time range
%          grd =  structure containing lon,lat

%  Example: uri='http://www.gri.msstate.edu/rsearch_data/nopp/bora_feb.nc';
%           var = 'temp';
%           lonLatRange = [13.0 16.0 41.0 42.0];   [minlon maxlon minlat maxlat]
%           dn1 = '14-Feb-2003 12:00:00';
%           dn2 = '[2003 2 16 14 0 0]';
%           [data, grd]=cf_subsetGrid(uri,var,lonLatRange,dn1, dn2)     or
%           [data, grd]=cf_subsetGrid(uri,var,lonLatRange,dn1)          or
%           [data, grd]=cf_subsetGrid(uri,var,lonLatRange)              or
%           [data, grd]=cf_subsetGrid(uri,var)                          or
%
% skbhate@ngi.msstate.edu
%
% import the NetCDF-Java methods 
import msstate.cstm.data.JDataset
import msstate.cstm.data.grid.JGeoGridDataset
import msstate.cstm.data.grid.JGeoGridUtil


if nargin < 2 
    disp('Check input arguments!')
    help cf_subsetGrid
    return
end

%initialize
%data (volume or subset)
data=[]; 
%structure containing lon,lat
grd.lat=[];
grd.lon=[];
grd.jdmat=[];

jdmat = [];

try 
    % open CF-compliant NetCDF File as a Common Data Model (CDM) "Grid Dataset"
        GridData = JDataset(uri);
        GeoGridData = JGeoGridDataset(GridData.getGridDataset(),var);
        if ( ~isa(GeoGridData.getGeoGrid(), 'ucar.nc2.dt.grid.GeoGrid') )  % make sure it's gridded data so that we can subset
            disp('Attempt to subset non-gridded variable.');
            help cf_subsetGrid
            return
        end
        GridCoordSys = GeoGridData.getCoordinateSystem();
        %get original geo grid associate with variable
        origGrid = GeoGridData.getGeoGrid();
        GridUtil = JGeoGridUtil(origGrid);
        
    switch nargin
      case 2        
        % read the volume data (3D). All times.
        myData = squeeze(GeoGridData.readVarData());     
        
      case 3 
        tindex = int32([]);   
         %subset the orginal geoGrid based on lat-lon range. tindex = null;
        subGrid = GridUtil.subsetGrid(lonLatRange,[]);       
        %set subset grid as new GeoGrid
        GeoGridData.setGeoGrid(subGrid);
        %get the data associated with new subGrid.
        myData = squeeze(GeoGridData.readVarData());  
        
      case 4 
        doSubset(dn1);        
       case 5 
        doSubset(dn1,dn2);
        
      otherwise, error('MATLAB:cf_subsetGrid:Nargin',...
                        'Incorrect number of arguments'); 
    end
    
     
    switch nargout
      case 1
        data = myData;
      case 2
        data = myData;
        % get coordinate axes associated with subGrid
        GridCoordData = GeoGridData.getGridCoordinateData();       
        lat=squeeze(GridCoordData.getLatAxis());
        lon=squeeze(GridCoordData.getLonAxis()); 
        %get time dates
        if ( GridCoordSys.hasTimeAxis() && GridCoordSys.isDate() ) 
            jdmat = nj_time(GeoGridData);  % get the subset of time dates
        end
        % Stuff into grid structure       
         grd.lon=lon;
         grd.lat=lat;
         grd.jdmat = jdmat;
      otherwise, error('MATLAB:cf_subsetGrid:Nargout',...
                        'Incorrect number of output arguments');
    end
    
    %cleanup
    GridData.close();
    
catch
    %gets the last error generated 
    err = lasterror();    
    disp(err.message); 
end

    %subset function based on time range.
    function doSubset(dn1,dn2)
        tindex = int32([]);        
        %subset the orginal geoGrid based on lat-lon and time range.       
        if ( GridCoordSys.hasTimeAxis() && GridCoordSys.isDate() ) 
            jdmat = nj_time(GeoGridData);  % get all the time dates
            switch nargin
                case 1
                    tindex = date_index(jdmat, dn1);
                case 2
                    tindex = date_index(jdmat, dn1,dn2);      
            end            
            if ~isempty(tindex)
                subGrid = GridUtil.subsetGrid(lonLatRange,tindex-1);
            else
                disp('No time range found.');
                subGrid = GridUtil.subsetGrid(lonLatRange);
            end
        else
            disp('No time axis associated with the dataset');
            subGrid = GridUtil.subsetGrid(lonLatRange);
        end        
        %set subset grid as new GeoGrid
        GeoGridData.setGeoGrid(subGrid);
        %get the data associated with new subGrid.
        myData = squeeze(GeoGridData.readVarData());
        % get new subset time dates
        %jdmat = nj_time(GridData, var);  % get all the time dates
    end
    
end