function data = DiagLoadMNC(fln,dat,dad,grd,iter,filesuffix);

% Function: DiagLoad
% Author:   Daniel Enderton
%
% Input Fields:
%
%   Field       Type        (Brief) Description
%   -----------------------------------------------------------------------
%   fln         string      Field name
%   dat         string      Data type ('Tav', 'Int',...)
%   dad         string      Data directory.
%   grd         string      Grid data directory.
%   iter        string      Iterations for analysis.
%   filesuffix  string      Field name suffix ('tave','',...)
%
% Output Fields:
%
%   Field       Type        (Brief) Description
%   -----------------------------------------------------------------------
%   data        cell array  Loaded data.
%   datatime    array       Model time of data in months.
%
% Descripton:
%   This function loads MNC data.  The field to load is defined by 'fln'.
%   Typically, whatever 'fln' is set to is the fields that is read in,
%   though the AIM physics parameters in a noteable exception to this.  The
%   data must be all located in the folder specified by 'dad', and this
%   string must end with a '/', and similarly for the grid data.


% Load parameters (here only general diagnostics).
DiagGenParam;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             Read in data                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generally, the data is read in by using the 'rdmds' function (see MITgcm)
% with FNAME as [fln,filesuffix], and ITER = iter.  There are some
% exceptions, such as the AIM physics parameters.  Some of the fields then
% require some immediate manipulation.  Examples are U (instantaneous) and
% uVel (time-averaged) which are immediately converted to a lat-lon grid
% because the 'U' and 'V' files of the cube-sphere do NOT correspond to the
% typical zonal and meridional winds.  Some other fields, such as 'KEpri'
% require some simple calculations which are immediately done.

if isequal(dat,'Tav')
    if isequal(fln,'phiHyd')
        ncroot = 'phiHyd.';
        filesuffix = '';
    else
        ncroot = 'tave.';
    end
elseif isequal(dat,'Int')
    ncroot = 'state.';
else
    error(['Unrecognized DataType:  ',dat]);
end

if ismember(fln,fields2D)
    data = zeros([length(iter),hres,faces*hres]);
elseif ismember(fln,fields3D)
    if ismember(fln,{'U','uVel'})
        data = zeros([length(iter),vres,hres,faces*hres+1]);
    elseif ismember(fln,{'V','vVel'})
        data = zeros([length(iter),vres,hres+1,faces*hres]);
    else
        data = zeros([length(iter),vres,hres,faces*hres]);
    end
else
    error(['Unrecognized field name:  ',fln]);
end

ncdir = dir(dad);
datafound = zeros(size(iter));
for ifile = 1:length(ncdir)
    ncfile = ncdir(ifile).name;
    if ~isempty(strfind(ncfile,ncroot))
        nc=netcdf([dad,ncfile],'read');
        it = nc{'iter'}(:);
        face = nc.('tile_number')(:);
        xindex = [hres*(face-1)+1:hres*face];
        [mncindex,dataindex] = ismember(it,iter);
        if sum(mncindex) > 0
            dataindex = dataindex(find(dataindex));
            mncindex = find(mncindex);
            datafound(dataindex) = datafound(dataindex) + 1;
            mncdata = nc{[fln,filesuffix]};
            if ismember(fln,fields2D)
                data(dataindex,:,xindex) = mncdata(mncindex,1:hres,1:hres);
            elseif ismember(fln,fields3D)
                if ismember(fln,{'U','uVel'})
                    data = zeros([length(iter),vres,hres,faces*(hres+1)]);
                    data(dataindex,:,:,xindex) = mncdata(mncindex,:,1:hres,1:hres);
                elseif ismember(fln,{'V','vVel'})
                    data(dataindex,:,:,xindex) = mncdata(mncindex,:,:,:);
                else
                    data(dataindex,:,:,xindex) = mncdata(mncindex,:,1:hres,1:hres);
                end
            end
        end
        close(nc);
    end
end


if ~( isequal(datafound,  faces*ones(size(iter))) | ...
      isequal(datafound,2*faces*ones(size(iter)))       )  % Embarassing fix!!! 
    disp(['Iterations desired:   ',mat2str(iter)]);
    disp(['Faces per iteration:  ',mat2str(datafound)]);
    error(['Missing Data!  (See output above for details).']);
end

data = permute(data,[length(size(data)):-1:1]);