function [data,grd]=cf_tslice(uri,var,iTime, level)
% CF_TSLICE: Get data and coordinates from a CF-compliant file at a specific time step
% [data,grd]=cf_tslice(uri,var,[iTime], [level]);
% Inputs: uri= CF-compliant NetCDF, OpenDAP or NCML file
%         var = variable to slice
%         iTime = time step to extract data 
%         level = vertical level (if not supplied, volume data is retrieved.)
%            
% Outputs: data = data          
%          grd   = structure containing lon,lat,z,jdmat (Matlab datenum)
%           
%
%  Example: uri ='http://stellwagen.er.usgs.gov/cgi-bin/nph-dods/models/adria/roms_sed/bora_feb.nc';% NetCDF/OpenDAP/NcML file
%           [data,grd]=cf_tslice(uri,'temp',2, 14) - Retrieve data from level 14 at time step 2
%           [data,grd]=cf_tslice(uri,'temp',2) - Retrieve 3D data at time step 2 
%           [data,grd]=cf_tslice(uri,'h') - Retrieve 2D non time dependent array           

% 
% rsignel1@usgs.gov
% skbhate@ngi.msstate.edu
%
% import the NetCDF-Java methods 
import msstate.cstm.data.JDataset
import msstate.cstm.data.grid.JGeoGridDataset
import msstate.cstm.util.JTime

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

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

try 
    % open CF-compliant NetCDF File as a Common Data Model (CDM) "Grid Dataset"
        GridData = JDataset(uri);
        %associated geo grid 
        GeoGridData = JGeoGridDataset(GridData.getGridDataset(),var); 
        % get access to nc2.dt.grid.GridCoordSys;
        GridCoordSys = GeoGridData.getCoordinateSystem();
        
    switch nargin
      case 2
        myLevel = 0;
        myTime = 0;
        % read the volume data (3D). All times.
        myData = squeeze(GeoGridData.readVarData());   
      case 3         
        myLevel = 0;
        myTime = iTime;
        % read the volume data at given time index(3D)
        myData = squeeze(GeoGridData.readVarData(iTime-1)); % zero based!!!       
      case 4       
        myLevel = level;
        myTime = iTime;
         % read a Y-X 'horizontal slice' at the given time and specified vertical index (2D)
        myData = squeeze(GeoGridData.readYXData(iTime-1,myLevel-1));  % zero based!!! 
      otherwise, error('MATLAB:cf_tslice:Nargin',...
                        'Incorrect number of arguments'); 
    end

    switch nargout
      case 1
        data = myData;
      case 2
        data = myData;
        % get the coordinate data object for this geo grid:
        GridCoordData = GeoGridData.getGridCoordinateData(); 
        % get coordinate axes        
        lat=squeeze(GridCoordData.getLatAxis());
        lon=squeeze(GridCoordData.getLonAxis());
        % Stuff into grid structure
        grd.lon=lon;
        grd.lat=lat;

        % if 3D, get vertical coordinate
        if GridCoordSys.hasVerticalAxis
            vt=GridCoordSys.getVerticalTransform();
            if (~isempty(vt))
                if (myTime) 
                    z = squeeze(GridCoordData.getVerticalCoordinates(iTime-1)); % vertical transform
                else 
                    z = [];
                end
                % Has vertical axis but no transform. 1D coordinate axes.
            else       
                z = squeeze(GridCoordData.getVerticalAxis1D()); % 2D  
            end
            if myLevel
                grd.z=squeeze(z(myLevel,:,:));
            else
                grd.z=z;
            end     
        end
        % get times.
        if GridCoordSys.hasTimeAxis()        
            % read times from file as gregorian dates
            timeDate = GridCoordSys.getTimeDates();                        
            if (nargin > 2) 
                grd.jdmat=JTime.getSerialDateNum(timeDate, iTime-1); %zero based!!! java class   
            else
                grd.jdmat=JTime.getSerialDateNum(timeDate);  % get all times.
            end
        end
        
      otherwise, error('MATLAB:cf_tslice:Nargout',...
                        'Incorrect number of output arguments'); 
    end

    %cleanup
    GridData.close();
catch 
    %gets the last error generated 
    err = lasterror();    
    disp(err.message); 
end