function [data,xax,yax,pltslc] = ...
    DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,...
              avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);

% Function: DiagSlice
% Author:   Daniel Enderton
%
% Input Fields:
%
%   Field        Type        (Brief) Description
%   -----------------------------------------------------------------------
%   data         array       Field name
%   fln          string      Field name
%   dad          string      Data directory.
%   grd          string      Grid data directory.
%   avg          string      Averaging scheme ('Ann','DJF', 'JJA',...)
%   slc          string      Slice type ('Zon','k=#',...)
%   pst          string      Plot style ('Con','Sur',...)
%   flu          string      Fluid ('A','O')
%   dfm          string      Data format ('MDS' or 'MNC')
%   LoadGridData 0/1         Optionally load grid data
%
% Output Fields:
%
%   Field       Type        (Brief) Description
%   -----------------------------------------------------------------------
%   data        array       Sliced data.
%
% Descripton:
%   This function slices the data according to SliceType 'slc'.  Depending
%   on the SliceType and PlotStyle, the data may be converted to a lat-lon
%   grid.  At the end there is also a range check.  If the value is out of
%   the user specified range given in 'DaigFieldParamA'.

% Load diagnostics parameters, grid data, mark data size
DiagGenParam;
DiagFieldParamA;
DiagFieldParamO;
DiagFieldParamC;
DiagFieldParamI;
[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
    DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
datasize = size(data);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Surface plot                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% All surface plots are direct cube sphere outputs, and therefore must have
% the dimensions defined in 'hres' in 'DiagLoadGridData'.  This is commonly
% [192,32].  If the 'pst' is 'Grd' or 'Int', everything is already all set
% since those handle the cubed sphere data.  If the 'pst' is 'Con' or 'Cnf',
% convert data to a lat-lon grid so that it can easily be contours.
if isequal(slc,'Sur')
    if ismember(fln,{'TX','TY','USTR','VSTR'}) | ~isequal(Vector,0)
        data = data'; xax = XL; yax = YL;
    elseif ~isequal(datasize,[faces*hres,hres]) 
        error('Incorrect dimensions for slc:  ',slc);
    else
        if ismember(pst,{'Grd','Int'})
            if     isequal(pst,'Grd'),
                xax = XG(1:faces*hres,1:hres);
                yax = YG(1:faces*hres,1:hres);
            elseif isequal(pst,'Int'),
                xax = XC(1:faces*hres,1:hres);
                yax = YC(1:faces*hres,1:hres);
            end
        elseif ismember(pst,{'Con','Cnf'})
            data = cube2latlon(XC,YC,data,XL,YL)';
            xax = XL; yax = YL;
        else
            error(['slc = ''Sur'' can not use pst:  ',pst]);
        end
    end
    pltslc='lonlat';
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Zonal average plot                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% When computing a zonal average, there are many different options for the
% initial data.  If the data is cube-sphere, it could be of the size
% [hres,z-axis length] or just [hres].  In either case, 'calc_ZonAv_CS' (a
% nifty zonally averaging script from JMC) is called to compute the zonal
% average for raw cube sphere data.  For horizontal velocities, the data
% could be of the size [length(XL),length(YL),z-axis length], in which case
% you just need to take the mean over the longitude axis (always 1).
elseif isequal(slc,'Zon')
    if isequal(datasize(1:2),[faces*hres,hres])
        if isequal(flu,'O'), nBas = 0; end
        if isequal(flu,'A'), nBas = 0; end
	    [data,dump1,dump2] = ...
	        calc_ZonAv_CS(data,kpr,kwr,nBas,XC,YC,XG,YG,RAC,dad,HFacC);
        if isequal(avg,'Tse')
            data=data(:,:,1);  xax=years;  yax=Ylat; pltslc='timlat';
        elseif isequal(avg,'Tyr')
            data=data(:,:,1);  xax=[0:12]; yax=Ylat; pltslc='timlat';
        elseif isequal(pst,'Lin')
            data=data(:,:,1)'; xax=Ylat;   yax=NaN;  pltslc='latfld';
        else
            data=data(:,:,1)'; xax=Ylat;   yax=ZC;   pltslc='lathgt';
        end
    elseif isequal(datasize(1:2),[length(XL),length(YL)])
        if ~isequal(pst,'Lin')
            data = squeeze(mean(data,1))'; xax=YL; yax=ZC;  pltslc='lathgt';
        else
            data = squeeze(mean(data,1))'; xax=YL; yax=NaN; pltslc='latfld';
        end
    else
        error('Incorrect dimensions for slc = ''Zon''');
    end

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                   i,j,k=#                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% These slicing settings just take a slice on one of the axes.  If the
% slice option is 'i=' or 'j=' and the data is cube-sphere, it is converted
% to lat-lon first and indexed then (thus the index number should be based
% off of the XL and YL interpolated grids defined in 'DiagGenParam'.  It
% would be nice to add functions that can slice by longitude, latitude, and
% height values directly.  Note that there is transposing of all the data
% to prepare it for contour plots (harmless is line plots).
elseif isequal(slc(1:2),'i=')
    ii = str2num(slc(3:end));
    data = cube2latlon(XC,YC,data,XL,YL);
    if ismember(fln,fields2D)
        data=squeeze(data(ii,:));
        xax=YL; yax=NaN; pltslc='latfld';
    elseif ismember(fln,fields3D)
        data=squeeze(data(ii,:,:))';
        xax=YL; yax=ZC;  pltslc='lathgt';
    end
    
elseif isequal(slc(1:2),'j=')
    jj = str2num(slc(3:end));
    data = cube2latlon(XC,YC,data,XL,YL);
    if ismember(fln,fields2D)
        data = squeeze(data(:,jj));
        xax=XL; yax=NaN; pltslc='lonfld';
    elseif ismember(fln,fields3D)
        data = squeeze(data(:,jj,:))';
        xax=XL; yax=ZC;  pltslc='lonhgt';
    end
    
elseif isequal(slc(1:2),'k=')
    kk = str2num(slc(3:end));
    data = squeeze(data(:,:,kk));
    if ismember(fln,{'U','V','uVel','vVel','fizhi_U','fizhi_V'})
        data = data'; xax = XL; yax = YL;
    elseif ismember(pst,{'Grd','Int'})
        if isequal(pst,'Grd'),
            xax = XG(1:faces*hres,1:hres);
            yax = YG(1:faces*hres,1:hres);
        elseif isequal(pst,'Int'),
            xax = XC(1:faces*hres,1:hres);
            yax = YC(1:faces*hres,1:hres);
        end
    elseif ismember(pst,{'Con','Cnf'})
        data = cube2latlon(XC,YC,data,XL,YL)'; xax = XL; yax = YL;
    end
    pltslc='lonlat';
    
else
    error(['Unrecognized SliceType:  ',slc]);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Check range                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load fixed and data ranges, throw a warning is data out of range.
datarange = [min(data(:)),max(data(:))];
try
    eval(['fixedrange = ',fln,'range',flu,';']);
    if datarange(1) < fixedrange(1) | ...
       datarange(2) > fixedrange(2)
        disp(['***Warning***  Value out of range for ',fln]);
        disp(['               Data range:  ',mat2str(datarange)]);
        disp(['               Fixed range:  ',mat2str(fixedrange)]);
    end
catch
    disp(['***Warning***  No range information found for ',fln]);
    disp(['               Data range:  ',mat2str(datarange)]);
end