function [data,xax,yax,time,pltslc] = ...
    DiagLoad(fln,exp,dat,dad,grd,itr,tst,flu,ddf,gdf,avg,slc,pst,...
             LoadGridData,DiagDebug,GridSuffix,ZcordFile,Index,Dim,...
             Grads,Year0Iter,Months);


% Load parameters.
DiagGenParam;
DiagFieldParamA;
DiagFieldParamO;
DiagFieldParamI;
DiagFieldParamC;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Prepare for the various data types, set time                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Currently the function can handle time averages (dat=='Tav'), which must
% be monthly mean, and instantaneous data (dat=='Ins').  There is also an
% option for generic data (dat=='Non') which, if the field name is exactly
% specified, is able to load any type file.
if isequal(dat,'Tav')
    filesuffix = 'tave';
    if isequal(flu,'C'), mnchandle  = 'cpl_tave.*';
    else               , mnchandle  = 'tave.*'    ; end
elseif isequal(dat,'Ins')
    filesuffix = '';  mnchandle  = 'state.*';
elseif ismember(dat,{'','Non','None'})
    filesuffix = '';  mnchandle  = '';
else
    error(['DataType not recognized:  ',datatype]);
end

% Time stuff (itr = 0 is start of new 360 day year).
if isempty(Months)
    months = (itr-Year0Iter)./(30.*24.*3600./tst);
else
    months = Months;
end
time = months/12;

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

% Open grads files
if ~isequal(Grads,0)
    if DiagDebug, disp(['  Debug -- Loading GRADS field.']); end
    [data,xax,yax,zax,months,time,dim] = DiagLoadGradsData(Grads,dad,fln,DiagDebug);
    if DiagDebug, disp(['  Debug -- ''data'' size after load:  ',mat2str(size(data))]); end
    data = DiagAverage(data,fln,avg,months,ddf,dim);
    if DiagDebug, disp(['  Debug -- ''data'' size after averaging:  ',mat2str(size(data))]); end
    if isequal(slc,'Sur')
        data = data';
        pltslc='lonlat';
    else
        error('Loading grads data only set to handle surface plots!');
    end
    

% Load AIM physics data.  The data is all stuffed into one file called
% 'aimPhy' which currently only comes in MDS format.  When loading a field,
% call the variable name that is used in the actual MITgcm code.  These
% variables are given in 'DiagGenParam'.  Here, we check that to see if the
% parameter is an AIM parameter, then load the 'aimPhy' file and extract
% the data with the proper index as specified in 'DiagGenParam'.  For total
% precipitation ('TotP') and evaporation minus precipitation ('EmP') some
% simple calculations of the AIM fields are performed and the U and V
% stresses (on atmosphere) must converted to a lat-lon grid.
elseif ismember(fln,aimparameters)
    if isequal(ddf,'MNC')
        error('Aim physics reader not set to hande ''MNC'' data format.')
    end
    data = DiagLoad_Local('aimPhy',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    if isequal(fln,'TotP')
        data =   data(:,:,PRECNV,:) ...
               + data(:,:,PRECLS,:);
    elseif isequal(fln,'EmP')
        data =   data(:,:,EVAP  ,:) ...
               - data(:,:,PRECNV,:) ...
               - data(:,:,PRECLS,:);
    elseif ismember(fln,{'USTR','VSTR'})
        [XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
            DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
        eval(['USTRdata = data(:,:,',num2str(USTR),',:);']);
        eval(['VSTRdata = data(:,:,',num2str(VSTR),',:);']);
        [USTRdata,VSTRdata]=uvcube2latlon(XC,YC,USTRdata,VSTRdata,XL,YL);
        if     isequal(fln,'USTR'), data = USTRdata; 
        elseif isequal(fln,'VSTR'), data = VSTRdata; end
    else
        eval(['data = data(:,:,',fln,',:);']);
    end
    data = squeeze(data(:,:,1,:));
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...
                  ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);

% Load coupled fields.  Like AIM, the data is all stuffed into one file
% called either 'cplFld.*' (MDS) or 'cpl_tave.*' (MNC).  Variables within
% the MDS files given in 'DiagGenParam', accessed directly with MNC.  TX
% and TY (stress on ocean) must be converted to lat-lon grid.
elseif ismember(fln,cplparameters)
    if ismember(fln,{'TX','TY'})
		[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
        DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
        if isequal(ddf,'MNC')
            TXdata = DiagLoad_Local('TX',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
            TYdata = DiagLoad_Local('TY',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        elseif isequal(ddf,'MDS')
            data = DiagLoad_Local('cplFld',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
            eval(['TXdata = squeeze(data(:,:,',num2str(TX),',:));']);
            eval(['TYdata = squeeze(data(:,:,',num2str(TY),',:));']);
        end
        [TXdata,TYdata]=uvcube2latlon(XC,YC,TXdata,TYdata,XL,YL);
        if     isequal(fln,'TX'), data = TXdata; 
        elseif isequal(fln,'TY'), data = TYdata; end
    else
        if isequal(ddf,'MNC')
            data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        elseif isequal(ddf,'MDS')
            data = DiagLoad_Local('cplFld',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
            eval(['data = squeeze(data(:,:,',fln,',:));']);
        end
        if isequal(fln,'SLP'), data = data ./ 100; end
    end
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...
                  ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);

% Load ice parameters.  Note that as the ice variable names appear to be in
% constant flux, this part quickly becomes out of date.
elseif ismember(fln,iceparameters)
    if isequal(dat,'Ins')
        data = DiagLoad_Local([fln,'-T'],dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    else
        data = DiagLoad_Local('thSIce_',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        eval(['data = data(:,:,',fln,',:);']);
        data = squeeze(data(:,:,1,:));
    end
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...
                  ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);

% Load generic fields where no modifications are needed.
elseif ismember(fln,{'S','T','Temp','aim_RH','phiHyd','Conv'})
    data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
                  flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);

% Read vertical velocities, eta.  Convert from [P] to [hPa] in atmosphere.
elseif ismember(fln,{'W','wVel','Eta','ETA'})
    data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    if isequal(flu,'A'), data = data./100; end
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
                  flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
    
% Load horizontal velocities, convert to lat-lon grid.
elseif ismember(fln,{'U','V','uVel','vVel'})
    if isequal(dat,'Tav')
        U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    else isequal(dat,'Ins')
        U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    end
    if isequal(ddf,'MNC')
        U = U(1:end-1,:,:,:);
        V = V(:,1:end-1,:,:);
    end
    U = DiagAverage(U,fln,avg,months,ddf,Dim);
    V = DiagAverage(V,fln,avg,months,ddf,Dim);
	[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
        DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
	[U,V]=uvcube2latlon(XC,YC,U,V,XL,YL);
    if ismember(fln,{'U','uVel'})
        data = U;
    elseif ismember(fln,{'V','vVel'})
        data = V;
    end
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
                  flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);

% Meridional overturning.
elseif ismember(fln,{'Bol','Psi','Res'})
    [XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
        DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
    if isequal(flu,'A'),     delM = - 1000*drF ./ g;
    elseif isequal(flu,'O'), delM = drF;             end
    if ismember(fln,{'Bol','Res'})
        kwx = DiagLoad_Local('GM_Kwx-T',dat,dad,grd,itr,'MDS','',mnchandle);
        kwy = DiagLoad_Local('GM_Kwy-T',dat,dad,grd,itr,'MDS','',mnchandle);
        kwx = DiagAverage(kwx,fln,avg,months,ddf,Dim);
        kwy = DiagAverage(kwy,fln,avg,months,ddf,Dim);
        [ub,vb]=calc_Bolus_CS(RAC,dxC,dyC,dxG,dyG,drC ,kwx,kwy,HFacW,HFacS);
        [psibol,mskG,ylat]=calc_PsiCube(delM,ub.*HFacW,vb.*HFacS,dxG,dyG);
    end
    if ismember(fln,{'Psi','Res'})
        if isequal(dat,'Tav')
            U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
            V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        else isequal(dat,'Ins')
            U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
            V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        end
        if isequal(ddf,'MNC')
            U = U(1:end-1,:,:,:);
            V = V(:,1:end-1,:,:);
        end
        U = DiagAverage(U,fln,avg,months,ddf,Dim);
        V = DiagAverage(V,fln,avg,months,ddf,Dim);
        [psiavg,mskG,ylat]=calc_PsiCube(delM,U.*HFacW,V.*HFacS,dxG,dyG);
    end
    if isequal(fln,'Psi')
        data = psiavg(2:end-1,:)';
    elseif isequal(fln,'Bol')
        data = psibol(2:end-1,:)';
    elseif isequal(fln,'Res')
        data = psibol(2:end-1,:)' + psiavg(2:end-1,:)';
    end
    if isequal(flu,'A'), data = data./1e6; end
    xax = ylat; yax = ZG; pltslc = 'lathgt';
    eval(['fixedrange = ',fln,'range',flu,';']);
    datarange = [min(data(:)),max(data(:))];
    if datarange(1) < fixedrange(1) | ...
        datarange(2) > fixedrange(2)
        disp(['***Warning***  Value out of range for ',fln]);
        disp(['               Calculated range:  ',mat2str(datarange)]);
        disp(['               Given range:  ',mat2str(fixedrange)]);
    end
  
% Read in variabilities.  Here we must make some simple calculations.  For
% all the variablilty fields there is are (field)^2 fields from which we
% can calculate the desired variabilities:  sqrt(field*field-(field)^2).
% Note that 'ETAstd' is scaled from [P] to [hPa].  For 'KEpri', we should
% use the values on either side of a grid, since it will be plotted from
% the grid center, but here we just use from one edge out of convenience.
elseif isequal(fln,'ETAstd')
    ETA  = DiagLoad_Local('ETA' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle)./100  ;
    ETA2 = DiagLoad_Local('Eta2',dat,dad,grd,itr,ddf,filesuffix,mnchandle)./10000;
    data = sqrt(abs(ETA.*ETA-ETA2));
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
                  flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
elseif isequal(fln,'Tstd')
    T  = DiagLoad_Local('T' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    TT = DiagLoad_Local('TT',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    data = sqrt(abs(T.*T-TT));
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
                  flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
elseif isequal(fln,'KEpri')
    U  = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    V  = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    UU = DiagLoad_Local('UU'  ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    VV = DiagLoad_Local('VV'  ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    if isequal(ddf,'MNC')
        U  = .5*(  U(2:end,:,:,:) +  U(1:end-1,:,:,:) );
        V  = .5*(  V(:,2:end,:,:) +  V(:,1:end-1,:,:) );
        %UU = .5*( UU(2:end,:,:,:) + UU(1:end-1,:,:,:) );  % UU and VV might be incorrect. 
        %VV = .5*( VV(:,2:end,:,:) + VV(:,1:end-1,:,:) );
    end
    U = DiagAverage(U,fln,avg,months,ddf,Dim);
    V = DiagAverage(V,fln,avg,months,ddf,Dim);
    UU = DiagAverage(UU,fln,avg,months,ddf,Dim);
    VV = DiagAverage(VV,fln,avg,months,ddf,Dim);
    data = sqrt(abs((U.*U + V.*V) - (UU + VV)));
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
                  flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
           
% If the field is not recognized, try best to try and open it.
else
    if DiagDebug, disp(['  Debug -- Unknown field, attempting to load.']); end
    data  = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
    if DiagDebug, disp(['  Debug -- ''data'' size after load:  ',mat2str(size(data))]); end
    if ~isequal(Index,0)
        if     isequal(Dim,2), data = squeeze(data(:,:,Index,:));
        elseif isequal(Dim,3), data = squeeze(data(:,:,:,Index,:)); end
    end
    if DiagDebug, disp(['  Debug -- ''data'' size after indexing:  ',mat2str(size(data))]); end
    data = DiagAverage(data,fln,avg,months,ddf,Dim);
    if DiagDebug, disp(['  Debug -- ''data'' size after average:  ',mat2str(size(data))]); end
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
                  flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
    if DiagDebug, disp(['  Debug -- ''data'' size after slice:  ',mat2str(size(data))]); end
    if DiagDebug, disp(['  Debug -- ''xax'' size:  ',mat2str(size(xax))]); end
    if DiagDebug, disp(['  Debug -- ''yax'' size:  ',mat2str(size(yax))]); end
    if DiagDebug, disp(['  Debug -- ''pltslc'':  ',pltslc]); end
end


%-------------------------------------------------------------------------%
%                            Local functions                              %
%-------------------------------------------------------------------------%

% Thus is merely a local function that calls the load data functions
% according to the DataFormat for the data specified by dfm.
function data = DiagLoad_Local(fln,dat,dad,grd,itr,dfm,filesuffix,mnchandle)
    
	if isequal(dfm,'MDS')
        data = rdmds([dad,'/',fln,filesuffix],itr);
	elseif isequal(dfm,'MNC')
        %[dad,mnchandle]
        iter = rdmnc_mod([dad,mnchandle],'T');
        iter = iter.T;
        [dummy,indecies] = ismember(itr,iter);
        data = rdmnc_mod([dad,mnchandle],[fln,filesuffix],'T',indecies);
        iter = data.T;
        if ~isequal(itr,iter'), error('Missing iterations in data!'); end
        eval(['data = data.',fln,filesuffix,';']);
	else
        error(['Unrecognized data type:  ',dfm]);
    end