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


% 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 isequal(dat(1:2),'Mn')
    filesuffix = '';
    if ismember(dat,{'MnA','MnO'}), mnchandle = 'monitor.*'; end
    if isequal(dat,'MnI'), mnchandle = 'monitor_sice.*'; end
    if isequal(dat,'MnL'), mnchandle = 'monitor_land.*'; end
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 monitor data.
elseif isequal(dat(1:2),'Mn')
    [data,time] = ...
        DiagLoadMonitor(fln,mnchandle,dad,itr,tst,SecPerYear,DiagDebug);
    xax = time; yax = NaN; pltslc = 'timfld';

% 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 handle ''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,trl,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,trl,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,trl,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,trl,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,trl,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','fizhi_U','fizhi_V'})
    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([fln(1:end-1),'U'],dat,dad,grd,itr,ddf,filesuffix,mnchandle);
        V = DiagLoad_Local([fln(1:end-1),'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','fizhi_U'})
        data = U;
    elseif ismember(fln,{'V','vVel','fizhi_V'})
        data = V;
    end
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,trl,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
    
    
% Actual temperature and moist potential temperature.  Actual temperature
% is calculated as T=Theta*(p/po)^(R/cp).  Moist potential temperature is
% calculated as Theta_e=Theta*e^(L*q*/cp*T) where q*=(R/Rv)*(es/p) and
% es=Ae^(beta*T).
elseif ismember(fln,{'ActT','ThetaEc'})
    if isequal('flu','O'),
        error('Calculation may only be done for atmosphere!');
    end
    [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,'MDS')
        theta = DiagLoad_Local('T',dat,dad,grd,itr,ddf,...
                               filesuffix,mnchandle);
    elseif isequal(ddf,'MNC')
        theta = DiagLoad_Local('Temp',dat,dad,grd,itr,...
                               ddf,filesuffix,mnchandle);
    end
    theta = DiagAverage(theta,fln,avg,months,ddf,Dim);
    if ismember(fln,{'ActT','ThetaEc'})
        pres = NaN.*zeros(size(theta));
        for iz=1:length(ZC), pres(:,:,iz)=ZC(iz); end
        temp=theta.*(pres./presrefA).^(RdA/cpA);
        if isequal(fln,'ActT'), data=temp; end
    end
    if isequal(fln,'ThetaEc')
        es=A_CC.*exp(Beta_CC.*(temp-K2C));
        qstar=(RdA/RvA).*(es./pres);
        data=theta.*exp(LHvapA.*qstar./cpA./temp);
    end
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,...
                  avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);

              
% 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,trl,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,trl,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,trl,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
    disp(['  Unknown field, attempting to load.']);
    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(Vector,1) | isequal(Vector,2) | isequal(Vector,3)
     if isequal(Vector,1) | isequal(Vector,3)
      if ~isequal(Index,0)
          if     isequal(Dim,2), U = squeeze(data(:,:,Index,:));
          elseif isequal(Dim,3), U = squeeze(data(:,:,:,Index,:)); end
      end
      if ~isequal(Mate,0)
          if     isequal(Dim,2), V = squeeze(data(:,:,Mate,:));
          elseif isequal(Dim,3), V = squeeze(data(:,:,:,Mate,:)); end
      end
     elseif isequal(Vector,2)
      if ~isequal(Index,0)
          if     isequal(Dim,2), V = squeeze(data(:,:,Index,:));
          elseif isequal(Dim,3), V = squeeze(data(:,:,:,Index,:)); end
      end
      if ~isequal(Mate,0)
          if     isequal(Dim,2), U = squeeze(data(:,:,Mate,:));
          elseif isequal(Dim,3), U = squeeze(data(:,:,:,Mate,:)); end
      end
     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);
     if isequal(Vector,1) | isequal(Vector,2)
	[U,V]=uvcube2latlon(XC,YC,U,V,XL,YL);
     end
     if isequal(Vector,1)
      data = U;
     elseif isequal(Vector,2)
      data = V;
     elseif isequal(Vector,3)
      if isequal(flu,'A'),     delM = - 1000*drF ./ g;
      elseif isequal(flu,'O'), delM = drF;             end
      [psiavg,mskG,ylat]=calc_PsiCube(delM,U.*HFacW,V.*HFacS,dxG,dyG);
      data = psiavg(2:end-1,:)';
      if isequal(flu,'A'), data = data./1e6; end
      xax = ylat; yax = ZG; pltslc = 'lathgt';
     end
     if DiagDebug, disp(['  Debug -- ''data'' size after vector manipulation:  ',mat2str(size(data))]); end
    else
% Sequence here for Vector = 0
     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);
%  Actual Temperature -- When the FieldName says ActT we assume that
%  the field that is pointed to in Index is potential temperature
     if isequal(FieldName,'ActT')
	[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
        DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
      pres = NaN.*zeros(size(data));
      for iz=1:length(ZC), pres(:,:,iz)=ZC(iz); end
      temp=data.*(pres./presrefA).^(RdA/cpA);
      data=temp; 
     end

     if DiagDebug, disp(['  Debug -- ''data'' size after average:  ',mat2str(size(data))]); end
    end
    if ~isequal(Vector,3)
    [data,xax,yax,pltslc] = ...
        DiagSlice(data,fln,trl,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
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')
        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