/[MITgcm]/MITgcm_contrib/enderton/Diagnostics/DiagLoad.m
ViewVC logotype

Diff of /MITgcm_contrib/enderton/Diagnostics/DiagLoad.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by molod, Mon Jan 31 21:26:32 2005 UTC revision 1.6 by enderton, Mon Feb 7 05:31:45 2005 UTC
# Line 1  Line 1 
1  function [data,xax,yax,time,pltslc] = ...  function [data,xax,yax,time,pltslc] = ...
2      DiagLoad(fln,exp,dat,dad,grd,itr,tst,flu,ddf,gdf,avg,slc,pst,...      DiagLoad(fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,avg,slc,pst,...
3               LoadGridData,DiagDebug,GridSuffix,ZcordFile,Index,Dim,...               LoadGridData,DiagDebug,GridSuffix,ZcordFile,Index,Dim,...
4               Grads,Year0Iter,Months);               Grads,Year0Iter,SecPerYear,Months);
5    
6    
7  % Load parameters.  % Load parameters.
# Line 16  DiagFieldParamC; Line 16  DiagFieldParamC;
16  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17    
18  % Currently the function can handle time averages (dat=='Tav'), which must  % Currently the function can handle time averages (dat=='Tav'), which must
19  % be monthly mean, and instantaneous data (dat=='Int').  There is also an  % be monthly mean, and instantaneous data (dat=='Ins').  There is also an
20  % option for generic data (dat=='Non') which, if the field name is exactly  % option for generic data (dat=='Non') which, if the field name is exactly
21  % specified, is able to load any type file.  % specified, is able to load any type file.
22  if isequal(dat,'Tav')  if isequal(dat,'Tav')
# Line 25  if isequal(dat,'Tav') Line 25  if isequal(dat,'Tav')
25      else               , mnchandle  = 'tave.*'    ; end      else               , mnchandle  = 'tave.*'    ; end
26  elseif isequal(dat,'Ins')  elseif isequal(dat,'Ins')
27      filesuffix = '';  mnchandle  = 'state.*';      filesuffix = '';  mnchandle  = 'state.*';
28    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  elseif ismember(dat,{'','Non','None'})  elseif ismember(dat,{'','Non','None'})
34      filesuffix = '';  mnchandle  = '';      filesuffix = '';  mnchandle  = '';
35  else  else
# Line 46  time = months/12; Line 51  time = months/12;
51  % Open grads files  % Open grads files
52  if ~isequal(Grads,0)  if ~isequal(Grads,0)
53      if DiagDebug, disp(['  Debug -- Loading GRADS field.']); end      if DiagDebug, disp(['  Debug -- Loading GRADS field.']); end
54      [data,xax,yax,zax,months,time,dim] = DiagLoadGradsData(Grads,dad,fln,DiagDebug);      [data,xax,yax,zax,months,time,dim] = ...
55            DiagLoadGradsData(Grads,dad,fln,DiagDebug);
56      if DiagDebug, disp(['  Debug -- ''data'' size after load:  ',mat2str(size(data))]); end      if DiagDebug, disp(['  Debug -- ''data'' size after load:  ',mat2str(size(data))]); end
57      data = DiagAverage(data,fln,avg,months,ddf,dim);      data = DiagAverage(data,fln,avg,months,ddf,dim);
58      if DiagDebug, disp(['  Debug -- ''data'' size after averaging:  ',mat2str(size(data))]); end      if DiagDebug, disp(['  Debug -- ''data'' size after averaging:  ',mat2str(size(data))]); end
# Line 56  if ~isequal(Grads,0) Line 62  if ~isequal(Grads,0)
62      else      else
63          error('Loading grads data only set to handle surface plots!');          error('Loading grads data only set to handle surface plots!');
64      end      end
65        
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    
72  % Load AIM physics data.  The data is all stuffed into one file called  % 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,  % 'aimPhy' which currently only comes in MDS format.  When loading a field,
# Line 69  if ~isequal(Grads,0) Line 80  if ~isequal(Grads,0)
80  % stresses (on atmosphere) must converted to a lat-lon grid.  % stresses (on atmosphere) must converted to a lat-lon grid.
81  elseif ismember(fln,aimparameters)  elseif ismember(fln,aimparameters)
82      if isequal(ddf,'MNC')      if isequal(ddf,'MNC')
83          error('Aim physics reader not set to hande ''MNC'' data format.')          error('Aim physics reader not set to handle ''MNC'' data format.')
84      end      end
85      data = DiagLoad_Local('aimPhy',dat,dad,grd,itr,ddf,filesuffix,mnchandle);      data = DiagLoad_Local('aimPhy',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
86      if isequal(fln,'TotP')      if isequal(fln,'TotP')
# Line 93  elseif ismember(fln,aimparameters) Line 104  elseif ismember(fln,aimparameters)
104      data = squeeze(data(:,:,1,:));      data = squeeze(data(:,:,1,:));
105      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
106      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
107          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
108                    ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
109    
110  % Load coupled fields.  Like AIM, the data is all stuffed into one file  % Load coupled fields.  Like AIM, the data is all stuffed into one file
# Line 126  elseif ismember(fln,cplparameters) Line 137  elseif ismember(fln,cplparameters)
137      end      end
138      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
139      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
140          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
141                    ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
142    
143  % Load ice parameters.  Note that as the ice variable names appear to be in  % Load ice parameters.  Note that as the ice variable names appear to be in
# Line 141  elseif ismember(fln,iceparameters) Line 152  elseif ismember(fln,iceparameters)
152      end      end
153      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
154      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
155          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
156                    ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
157    
158  % Load generic fields where no modifications are needed.  % Load generic fields where no modifications are needed.
# Line 149  elseif ismember(fln,{'S','T','Temp','aim Line 160  elseif ismember(fln,{'S','T','Temp','aim
160      data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);      data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
161      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
162      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
163          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
164                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
165    
166  % Read vertical velocities, eta.  Convert from [P] to [hPa] in atmosphere.  % Read vertical velocities, eta.  Convert from [P] to [hPa] in atmosphere.
# Line 158  elseif ismember(fln,{'W','wVel','Eta','E Line 169  elseif ismember(fln,{'W','wVel','Eta','E
169      if isequal(flu,'A'), data = data./100; end      if isequal(flu,'A'), data = data./100; end
170      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
171      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
172          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
173                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
174            
175  % Load horizontal velocities, convert to lat-lon grid.  % Load horizontal velocities, convert to lat-lon grid.
# Line 166  elseif ismember(fln,{'U','V','uVel','vVe Line 177  elseif ismember(fln,{'U','V','uVel','vVe
177      if isequal(dat,'Tav')      if isequal(dat,'Tav')
178          U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);          U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
179          V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);          V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
180      else isequal(dat,'Int')      else isequal(dat,'Ins')
181          U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);          U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
182          V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);          V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
183      end      end
# Line 185  elseif ismember(fln,{'U','V','uVel','vVe Line 196  elseif ismember(fln,{'U','V','uVel','vVe
196          data = V;          data = V;
197      end      end
198      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
199          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
200                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
201    
202  % Meridional overturning.  % Meridional overturning.
# Line 206  elseif ismember(fln,{'Bol','Psi','Res'}) Line 217  elseif ismember(fln,{'Bol','Psi','Res'})
217          if isequal(dat,'Tav')          if isequal(dat,'Tav')
218              U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);              U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
219              V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);              V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
220          else isequal(dat,'Int')          else isequal(dat,'Ins')
221              U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);              U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
222              V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);              V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
223          end          end
# Line 235  elseif ismember(fln,{'Bol','Psi','Res'}) Line 246  elseif ismember(fln,{'Bol','Psi','Res'})
246          disp(['               Calculated range:  ',mat2str(datarange)]);          disp(['               Calculated range:  ',mat2str(datarange)]);
247          disp(['               Given range:  ',mat2str(fixedrange)]);          disp(['               Given range:  ',mat2str(fixedrange)]);
248      end      end
249          
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=Theta*e^(L*q*/cp*T) where q*=(R/Rv)*(es/p) and
254    % es=Ae^(beta*T).
255    elseif ismember(fln,{'ActT','ThetaEc'})
256        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        if ismember(fln,{'ActT','ThetaEc'})
271            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        if isequal(fln,'ThetaEc')
277            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                      avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
284    
285                  
286  % Read in variabilities.  Here we must make some simple calculations.  For  % 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  % all the variablilty fields there is are (field)^2 fields from which we
288  % can calculate the desired variabilities:  sqrt(field*field-(field)^2).  % can calculate the desired variabilities:  sqrt(field*field-(field)^2).
# Line 248  elseif isequal(fln,'ETAstd') Line 295  elseif isequal(fln,'ETAstd')
295      data = sqrt(abs(ETA.*ETA-ETA2));      data = sqrt(abs(ETA.*ETA-ETA2));
296      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
297      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
298          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
299                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
300  elseif isequal(fln,'Tstd')  elseif isequal(fln,'Tstd')
301      T  = DiagLoad_Local('T' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);      T  = DiagLoad_Local('T' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
# Line 256  elseif isequal(fln,'Tstd') Line 303  elseif isequal(fln,'Tstd')
303      data = sqrt(abs(T.*T-TT));      data = sqrt(abs(T.*T-TT));
304      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
305      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
306          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
307                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
308  elseif isequal(fln,'KEpri')  elseif isequal(fln,'KEpri')
309      U  = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);      U  = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
# Line 275  elseif isequal(fln,'KEpri') Line 322  elseif isequal(fln,'KEpri')
322      VV = DiagAverage(VV,fln,avg,months,ddf,Dim);      VV = DiagAverage(VV,fln,avg,months,ddf,Dim);
323      data = sqrt(abs((U.*U + V.*V) - (UU + VV)));      data = sqrt(abs((U.*U + V.*V) - (UU + VV)));
324      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
325          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
326                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
327                        
328  % If the field is not recognized, try best to try and open it.  % If the field is not recognized, try best to try and open it.
329  else  else
330      if DiagDebug, disp(['  Debug -- Unknown field, attempting to load.']); end      disp(['  Unknown field, attempting to load.']);
331      data  = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);      data  = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
332      if DiagDebug, disp(['  Debug -- ''data'' size after load:  ',mat2str(size(data))]); end      if DiagDebug, disp(['  Debug -- ''data'' size after load:  ',mat2str(size(data))]); end
333      if ~isequal(Index,0)      if ~isequal(Index,0)
# Line 291  else Line 338  else
338      data = DiagAverage(data,fln,avg,months,ddf,Dim);      data = DiagAverage(data,fln,avg,months,ddf,Dim);
339      if DiagDebug, disp(['  Debug -- ''data'' size after average:  ',mat2str(size(data))]); end      if DiagDebug, disp(['  Debug -- ''data'' size after average:  ',mat2str(size(data))]); end
340      [data,xax,yax,pltslc] = ...      [data,xax,yax,pltslc] = ...
341          DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...          DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
342                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);                    flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
343      if DiagDebug, disp(['  Debug -- ''data'' size after slice:  ',mat2str(size(data))]); end      if DiagDebug, disp(['  Debug -- ''data'' size after slice:  ',mat2str(size(data))]); end
344      if DiagDebug, disp(['  Debug -- ''xax'' size:  ',mat2str(size(xax))]); end      if DiagDebug, disp(['  Debug -- ''xax'' size:  ',mat2str(size(xax))]); end
# Line 311  function data = DiagLoad_Local(fln,dat,d Line 358  function data = DiagLoad_Local(fln,dat,d
358          if isequal(dfm,'MDS')          if isequal(dfm,'MDS')
359          data = rdmds([dad,'/',fln,filesuffix],itr);          data = rdmds([dad,'/',fln,filesuffix],itr);
360          elseif isequal(dfm,'MNC')          elseif isequal(dfm,'MNC')
         %[dad,mnchandle]  
361          iter = rdmnc_mod([dad,mnchandle],'T');          iter = rdmnc_mod([dad,mnchandle],'T');
362          iter = iter.T;          iter = iter.T;
363          [dummy,indecies] = ismember(itr,iter);          [dummy,indecies] = ismember(itr,iter);
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22