enderton/Diagnostics/DiagLoad.m

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,Vector,FieldName);

% Load coupled fields.  Like AIM, the data is all stuffed into one file
% called either 'cplFldtave.*' (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,Vector,FieldName);

% 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,Vector,FieldName);

% 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,Vector,FieldName);

% 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,Vector,FieldName);
    
% 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,Vector,FieldName);

% 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,Vector,FieldName);

              
% 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,Vector,FieldName);
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,Vector,FieldName);
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 = U(1:end-1,:,:,:); UU = UU(1:end-1,:,:,:);
        V = V(:,1:end-1,:,:); VV = 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,Vector,FieldName);
           
% 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,Vector,FieldName);
    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')
        data = rdmnc_mod2([dad,mnchandle],[fln,filesuffix],'iter','T',itr);
        if ~isequal(itr,data.iter'), error('Missing iterations in data!'); end
        eval(['data = data.',fln,filesuffix,';']);
        else
        error(['Unrecognized data type:  ',dfm]);
    end
1	enderton	1.1	function [data,xax,yax,time,pltslc] = ...
2	enderton	1.4	DiagLoad(fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,avg,slc,pst,...
3	molod	1.8	LoadGridData,DiagDebug,GridSuffix,ZcordFile,Index,Dim,Vector,Mate,...
4			Grads,Year0Iter,SecPerYear,Months,FieldName);
5	enderton	1.1
6
7			% Load parameters.
8			DiagGenParam;
9			DiagFieldParamA;
10			DiagFieldParamO;
11			DiagFieldParamI;
12			DiagFieldParamC;
13
14			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15			% Prepare for the various data types, set time %
16			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17
18			% Currently the function can handle time averages (dat=='Tav'), which must
19	molod	1.3	% be monthly mean, and instantaneous data (dat=='Ins'). There is also an
20	enderton	1.1	% option for generic data (dat=='Non') which, if the field name is exactly
21			% specified, is able to load any type file.
22			if isequal(dat,'Tav')
23			filesuffix = 'tave';
24			if isequal(flu,'C'), mnchandle = 'cpl_tave.*';
25			else , mnchandle = 'tave.*' ; end
26			elseif isequal(dat,'Ins')
27			filesuffix = ''; mnchandle = 'state.*';
28	enderton	1.6	elseif isequal(dat(1:2),'Mn')
29			filesuffix = '';
30			if ismember(dat,{'MnA','MnO'}), mnchandle = 'monitor.*'; end
31			if isequal(dat,'MnI'), mnchandle = 'monitor_sice.*'; end
32			if isequal(dat,'MnL'), mnchandle = 'monitor_land.*'; end
33	enderton	1.1	elseif ismember(dat,{'','Non','None'})
34			filesuffix = ''; mnchandle = '';
35			else
36			error(['DataType not recognized: ',datatype]);
37			end
38
39			% Time stuff (itr = 0 is start of new 360 day year).
40			if isempty(Months)
41			months = (itr-Year0Iter)./(30.24.3600./tst);
42			else
43			months = Months;
44			end
45			time = months/12;
46
47			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48			% Read in data %
49			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50
51			% Open grads files
52			if ~isequal(Grads,0)
53			if DiagDebug, disp([' Debug -- Loading GRADS field.']); end
54	enderton	1.6	[data,xax,yax,zax,months,time,dim] = ...
55			DiagLoadGradsData(Grads,dad,fln,DiagDebug);
56	enderton	1.1	if DiagDebug, disp([' Debug -- ''data'' size after load: ',mat2str(size(data))]); end
57			data = DiagAverage(data,fln,avg,months,ddf,dim);
58			if DiagDebug, disp([' Debug -- ''data'' size after averaging: ',mat2str(size(data))]); end
59			if isequal(slc,'Sur')
60			data = data';
61			pltslc='lonlat';
62			else
63			error('Loading grads data only set to handle surface plots!');
64			end
65	enderton	1.6
66			% Load monitor data.
67			elseif isequal(dat(1:2),'Mn')
68			[data,time] = ...
69			DiagLoadMonitor(fln,mnchandle,dad,itr,tst,SecPerYear,DiagDebug);
70			xax = time; yax = NaN; pltslc = 'timfld';
71	enderton	1.1
72			% Load AIM physics data. The data is all stuffed into one file called
73			% 'aimPhy' which currently only comes in MDS format. When loading a field,
74			% call the variable name that is used in the actual MITgcm code. These
75			% variables are given in 'DiagGenParam'. Here, we check that to see if the
76			% parameter is an AIM parameter, then load the 'aimPhy' file and extract
77			% the data with the proper index as specified in 'DiagGenParam'. For total
78			% precipitation ('TotP') and evaporation minus precipitation ('EmP') some
79			% simple calculations of the AIM fields are performed and the U and V
80			% stresses (on atmosphere) must converted to a lat-lon grid.
81			elseif ismember(fln,aimparameters)
82			if isequal(ddf,'MNC')
83	enderton	1.6	error('Aim physics reader not set to handle ''MNC'' data format.')
84	enderton	1.1	end
85			data = DiagLoad_Local('aimPhy',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
86			if isequal(fln,'TotP')
87			data = data(:,:,PRECNV,:) ...
88			+ data(:,:,PRECLS,:);
89			elseif isequal(fln,'EmP')
90			data = data(:,:,EVAP ,:) ...
91			- data(:,:,PRECNV,:) ...
92			- data(:,:,PRECLS,:);
93			elseif ismember(fln,{'USTR','VSTR'})
94			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
95			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
96			eval(['USTRdata = data(:,:,',num2str(USTR),',:);']);
97			eval(['VSTRdata = data(:,:,',num2str(VSTR),',:);']);
98			[USTRdata,VSTRdata]=uvcube2latlon(XC,YC,USTRdata,VSTRdata,XL,YL);
99			if isequal(fln,'USTR'), data = USTRdata;
100			elseif isequal(fln,'VSTR'), data = VSTRdata; end
101			else
102			eval(['data = data(:,:,',fln,',:);']);
103			end
104			data = squeeze(data(:,:,1,:));
105			data = DiagAverage(data,fln,avg,months,ddf,Dim);
106			[data,xax,yax,pltslc] = ...
107	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
108	molod	1.9	ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
109	enderton	1.1
110			% Load coupled fields. Like AIM, the data is all stuffed into one file
111	enderton	1.10	% called either 'cplFldtave.' (MDS) or 'cpl_tave.' (MNC). Variables within
112	enderton	1.1	% the MDS files given in 'DiagGenParam', accessed directly with MNC. TX
113			% and TY (stress on ocean) must be converted to lat-lon grid.
114			elseif ismember(fln,cplparameters)
115			if ismember(fln,{'TX','TY'})
116			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
117			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
118			if isequal(ddf,'MNC')
119			TXdata = DiagLoad_Local('TX',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
120			TYdata = DiagLoad_Local('TY',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
121			elseif isequal(ddf,'MDS')
122			data = DiagLoad_Local('cplFld',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
123			eval(['TXdata = squeeze(data(:,:,',num2str(TX),',:));']);
124			eval(['TYdata = squeeze(data(:,:,',num2str(TY),',:));']);
125			end
126			[TXdata,TYdata]=uvcube2latlon(XC,YC,TXdata,TYdata,XL,YL);
127			if isequal(fln,'TX'), data = TXdata;
128			elseif isequal(fln,'TY'), data = TYdata; end
129			else
130			if isequal(ddf,'MNC')
131			data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
132			elseif isequal(ddf,'MDS')
133			data = DiagLoad_Local('cplFld',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
134			eval(['data = squeeze(data(:,:,',fln,',:));']);
135			end
136			if isequal(fln,'SLP'), data = data ./ 100; end
137			end
138			data = DiagAverage(data,fln,avg,months,ddf,Dim);
139			[data,xax,yax,pltslc] = ...
140	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
141	molod	1.9	ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
142	enderton	1.1
143			% Load ice parameters. Note that as the ice variable names appear to be in
144			% constant flux, this part quickly becomes out of date.
145			elseif ismember(fln,iceparameters)
146			if isequal(dat,'Ins')
147			data = DiagLoad_Local([fln,'-T'],dat,dad,grd,itr,ddf,filesuffix,mnchandle);
148			else
149			data = DiagLoad_Local('thSIce_',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
150			eval(['data = data(:,:,',fln,',:);']);
151			data = squeeze(data(:,:,1,:));
152			end
153			data = DiagAverage(data,fln,avg,months,ddf,Dim);
154			[data,xax,yax,pltslc] = ...
155	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
156	molod	1.9	ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
157	enderton	1.1
158			% Load generic fields where no modifications are needed.
159			elseif ismember(fln,{'S','T','Temp','aim_RH','phiHyd','Conv'})
160			data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
161			data = DiagAverage(data,fln,avg,months,ddf,Dim);
162			[data,xax,yax,pltslc] = ...
163	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
164	molod	1.9	flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
165	enderton	1.1
166			% Read vertical velocities, eta. Convert from [P] to [hPa] in atmosphere.
167			elseif ismember(fln,{'W','wVel','Eta','ETA'})
168			data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
169			if isequal(flu,'A'), data = data./100; end
170			data = DiagAverage(data,fln,avg,months,ddf,Dim);
171			[data,xax,yax,pltslc] = ...
172	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
173	molod	1.9	flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
174	enderton	1.1
175			% Load horizontal velocities, convert to lat-lon grid.
176	molod	1.7	elseif ismember(fln,{'U','V','uVel','vVel','fizhi_U','fizhi_V'})
177	enderton	1.1	if isequal(dat,'Tav')
178			U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
179			V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
180	molod	1.3	else isequal(dat,'Ins')
181	molod	1.7	U = DiagLoad_Local([fln(1:end-1),'U'],dat,dad,grd,itr,ddf,filesuffix,mnchandle);
182			V = DiagLoad_Local([fln(1:end-1),'V'],dat,dad,grd,itr,ddf,filesuffix,mnchandle);
183	enderton	1.1	end
184			if isequal(ddf,'MNC')
185			U = U(1:end-1,:,:,:);
186			V = V(:,1:end-1,:,:);
187			end
188			U = DiagAverage(U,fln,avg,months,ddf,Dim);
189			V = DiagAverage(V,fln,avg,months,ddf,Dim);
190			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
191			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
192	enderton	1.10	%[U,V]=uvcube2latlon(XC,YC,U,V,XL,YL);
193	molod	1.7	if ismember(fln,{'U','uVel','fizhi_U'})
194	enderton	1.1	data = U;
195	molod	1.7	elseif ismember(fln,{'V','vVel','fizhi_V'})
196	enderton	1.1	data = V;
197			end
198			[data,xax,yax,pltslc] = ...
199	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
200	molod	1.9	flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
201	enderton	1.1
202			% Meridional overturning.
203			elseif ismember(fln,{'Bol','Psi','Res'})
204			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
205			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
206			if isequal(flu,'A'), delM = - 1000*drF ./ g;
207			elseif isequal(flu,'O'), delM = drF; end
208			if ismember(fln,{'Bol','Res'})
209			kwx = DiagLoad_Local('GM_Kwx-T',dat,dad,grd,itr,'MDS','',mnchandle);
210			kwy = DiagLoad_Local('GM_Kwy-T',dat,dad,grd,itr,'MDS','',mnchandle);
211			kwx = DiagAverage(kwx,fln,avg,months,ddf,Dim);
212			kwy = DiagAverage(kwy,fln,avg,months,ddf,Dim);
213			[ub,vb]=calc_Bolus_CS(RAC,dxC,dyC,dxG,dyG,drC ,kwx,kwy,HFacW,HFacS);
214			[psibol,mskG,ylat]=calc_PsiCube(delM,ub.HFacW,vb.HFacS,dxG,dyG);
215			end
216			if ismember(fln,{'Psi','Res'})
217			if isequal(dat,'Tav')
218			U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
219			V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
220	molod	1.3	else isequal(dat,'Ins')
221	enderton	1.1	U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
222			V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
223			end
224			if isequal(ddf,'MNC')
225			U = U(1:end-1,:,:,:);
226			V = V(:,1:end-1,:,:);
227			end
228			U = DiagAverage(U,fln,avg,months,ddf,Dim);
229			V = DiagAverage(V,fln,avg,months,ddf,Dim);
230			[psiavg,mskG,ylat]=calc_PsiCube(delM,U.HFacW,V.HFacS,dxG,dyG);
231			end
232			if isequal(fln,'Psi')
233			data = psiavg(2:end-1,:)';
234			elseif isequal(fln,'Bol')
235			data = psibol(2:end-1,:)';
236			elseif isequal(fln,'Res')
237			data = psibol(2:end-1,:)' + psiavg(2:end-1,:)';
238			end
239			if isequal(flu,'A'), data = data./1e6; end
240			xax = ylat; yax = ZG; pltslc = 'lathgt';
241			eval(['fixedrange = ',fln,'range',flu,';']);
242			datarange = [min(data(:)),max(data(:))];
243			if datarange(1) < fixedrange(1) \| ...
244			datarange(2) > fixedrange(2)
245			disp(['*Warning* Value out of range for ',fln]);
246			disp([' Calculated range: ',mat2str(datarange)]);
247			disp([' Given range: ',mat2str(fixedrange)]);
248			end
249	enderton	1.4
250
251			% Actual temperature and moist potential temperature. Actual temperature
252			% is calculated as T=Theta*(p/po)^(R/cp). Moist potential temperature is
253			% calculated as Theta_e=Thetae^(Lq/cpT) where q=(R/Rv)(es/p) and
254			% es=Ae^(beta*T).
255	enderton	1.5	elseif ismember(fln,{'ActT','ThetaEc'})
256	enderton	1.4	if isequal('flu','O'),
257			error('Calculation may only be done for atmosphere!');
258			end
259			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,...
260			HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
261			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
262			if isequal(ddf,'MDS')
263			theta = DiagLoad_Local('T',dat,dad,grd,itr,ddf,...
264			filesuffix,mnchandle);
265			elseif isequal(ddf,'MNC')
266			theta = DiagLoad_Local('Temp',dat,dad,grd,itr,...
267			ddf,filesuffix,mnchandle);
268			end
269			theta = DiagAverage(theta,fln,avg,months,ddf,Dim);
270	enderton	1.5	if ismember(fln,{'ActT','ThetaEc'})
271	enderton	1.4	pres = NaN.*zeros(size(theta));
272			for iz=1:length(ZC), pres(:,:,iz)=ZC(iz); end
273			temp=theta.*(pres./presrefA).^(RdA/cpA);
274			if isequal(fln,'ActT'), data=temp; end
275			end
276	enderton	1.5	if isequal(fln,'ThetaEc')
277	enderton	1.4	es=A_CC.exp(Beta_CC.(temp-K2C));
278			qstar=(RdA/RvA).*(es./pres);
279			data=theta.exp(LHvapA.qstar./cpA./temp);
280			end
281			[data,xax,yax,pltslc] = ...
282			DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,...
283	molod	1.9	avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
284	enderton	1.4
285
286	enderton	1.1	% Read in variabilities. Here we must make some simple calculations. For
287			% all the variablilty fields there is are (field)^2 fields from which we
288			% can calculate the desired variabilities: sqrt(field*field-(field)^2).
289			% Note that 'ETAstd' is scaled from [P] to [hPa]. For 'KEpri', we should
290			% use the values on either side of a grid, since it will be plotted from
291			% the grid center, but here we just use from one edge out of convenience.
292			elseif isequal(fln,'ETAstd')
293			ETA = DiagLoad_Local('ETA' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle)./100 ;
294			ETA2 = DiagLoad_Local('Eta2',dat,dad,grd,itr,ddf,filesuffix,mnchandle)./10000;
295			data = sqrt(abs(ETA.*ETA-ETA2));
296			data = DiagAverage(data,fln,avg,months,ddf,Dim);
297			[data,xax,yax,pltslc] = ...
298	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
299	molod	1.9	flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
300	enderton	1.1	elseif isequal(fln,'Tstd')
301			T = DiagLoad_Local('T' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
302			TT = DiagLoad_Local('TT',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
303			data = sqrt(abs(T.*T-TT));
304			data = DiagAverage(data,fln,avg,months,ddf,Dim);
305			[data,xax,yax,pltslc] = ...
306	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
307	molod	1.9	flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
308	enderton	1.1	elseif isequal(fln,'KEpri')
309			U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
310			V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
311			UU = DiagLoad_Local('UU' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
312			VV = DiagLoad_Local('VV' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
313			if isequal(ddf,'MNC')
314	enderton	1.10	U = U(1:end-1,:,:,:); UU = UU(1:end-1,:,:,:);
315			V = V(:,1:end-1,:,:); VV = VV(:,1:end-1,:,:);
316	enderton	1.1	end
317	enderton	1.10	U = DiagAverage(U ,fln,avg,months,ddf,Dim);
318			V = DiagAverage(V ,fln,avg,months,ddf,Dim);
319	enderton	1.1	UU = DiagAverage(UU,fln,avg,months,ddf,Dim);
320			VV = DiagAverage(VV,fln,avg,months,ddf,Dim);
321			data = sqrt(abs((U.U + V.V) - (UU + VV)));
322			[data,xax,yax,pltslc] = ...
323	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
324	molod	1.9	flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
325	enderton	1.1
326			% If the field is not recognized, try best to try and open it.
327			else
328	enderton	1.4	disp([' Unknown field, attempting to load.']);
329	enderton	1.1	data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
330			if DiagDebug, disp([' Debug -- ''data'' size after load: ',mat2str(size(data))]); end
331	molod	1.8	if isequal(Vector,1) \| isequal(Vector,2) \| isequal(Vector,3)
332			if isequal(Vector,1) \| isequal(Vector,3)
333			if ~isequal(Index,0)
334			if isequal(Dim,2), U = squeeze(data(:,:,Index,:));
335			elseif isequal(Dim,3), U = squeeze(data(:,:,:,Index,:)); end
336			end
337			if ~isequal(Mate,0)
338			if isequal(Dim,2), V = squeeze(data(:,:,Mate,:));
339			elseif isequal(Dim,3), V = squeeze(data(:,:,:,Mate,:)); end
340			end
341			elseif isequal(Vector,2)
342			if ~isequal(Index,0)
343			if isequal(Dim,2), V = squeeze(data(:,:,Index,:));
344			elseif isequal(Dim,3), V = squeeze(data(:,:,:,Index,:)); end
345			end
346			if ~isequal(Mate,0)
347			if isequal(Dim,2), U = squeeze(data(:,:,Mate,:));
348			elseif isequal(Dim,3), U = squeeze(data(:,:,:,Mate,:)); end
349			end
350			end
351			U = DiagAverage(U,fln,avg,months,ddf,Dim);
352			V = DiagAverage(V,fln,avg,months,ddf,Dim);
353			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
354			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
355			if isequal(Vector,1) \| isequal(Vector,2)
356			[U,V]=uvcube2latlon(XC,YC,U,V,XL,YL);
357			end
358			if isequal(Vector,1)
359			data = U;
360			elseif isequal(Vector,2)
361			data = V;
362			elseif isequal(Vector,3)
363			if isequal(flu,'A'), delM = - 1000*drF ./ g;
364			elseif isequal(flu,'O'), delM = drF; end
365			[psiavg,mskG,ylat]=calc_PsiCube(delM,U.HFacW,V.HFacS,dxG,dyG);
366			data = psiavg(2:end-1,:)';
367			if isequal(flu,'A'), data = data./1e6; end
368			xax = ylat; yax = ZG; pltslc = 'lathgt';
369			end
370			if DiagDebug, disp([' Debug -- ''data'' size after vector manipulation: ',mat2str(size(data))]); end
371			else
372			% Sequence here for Vector = 0
373			if ~isequal(Index,0)
374	enderton	1.1	if isequal(Dim,2), data = squeeze(data(:,:,Index,:));
375			elseif isequal(Dim,3), data = squeeze(data(:,:,:,Index,:)); end
376	molod	1.8	end
377			if DiagDebug, disp([' Debug -- ''data'' size after indexing: ',mat2str(size(data))]); end
378			data = DiagAverage(data,fln,avg,months,ddf,Dim);
379			% Actual Temperature -- When the FieldName says ActT we assume that
380			% the field that is pointed to in Index is potential temperature
381			if isequal(FieldName,'ActT')
382			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
383			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
384			pres = NaN.*zeros(size(data));
385			for iz=1:length(ZC), pres(:,:,iz)=ZC(iz); end
386			temp=data.*(pres./presrefA).^(RdA/cpA);
387			data=temp;
388			end
389
390			if DiagDebug, disp([' Debug -- ''data'' size after average: ',mat2str(size(data))]); end
391	enderton	1.1	end
392	molod	1.8	if ~isequal(Vector,3)
393	enderton	1.1	[data,xax,yax,pltslc] = ...
394	enderton	1.4	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
395	molod	1.9	flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile,Vector,FieldName);
396	enderton	1.1	if DiagDebug, disp([' Debug -- ''data'' size after slice: ',mat2str(size(data))]); end
397			if DiagDebug, disp([' Debug -- ''xax'' size: ',mat2str(size(xax))]); end
398			if DiagDebug, disp([' Debug -- ''yax'' size: ',mat2str(size(yax))]); end
399			if DiagDebug, disp([' Debug -- ''pltslc'': ',pltslc]); end
400	molod	1.8	end
401	enderton	1.1	end
402
403			%-------------------------------------------------------------------------%
404			% Local functions %
405			%-------------------------------------------------------------------------%
406
407			% Thus is merely a local function that calls the load data functions
408			% according to the DataFormat for the data specified by dfm.
409			function data = DiagLoad_Local(fln,dat,dad,grd,itr,dfm,filesuffix,mnchandle)
410			if isequal(dfm,'MDS')
411			data = rdmds([dad,'/',fln,filesuffix],itr);
412			elseif isequal(dfm,'MNC')
413	enderton	1.10	data = rdmnc_mod2([dad,mnchandle],[fln,filesuffix],'iter','T',itr);
414			if ~isequal(itr,data.iter'), error('Missing iterations in data!'); end
415	enderton	1.1	eval(['data = data.',fln,filesuffix,';']);
416			else
417			error(['Unrecognized data type: ',dfm]);
418			end