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

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

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


Revision 1.1 - (hide annotations) (download)
Mon Jan 31 15:43:26 2005 UTC (20 years, 5 months ago) by enderton
Branch: MAIN
 o Initial check in.

1 enderton 1.1 function [data,xax,yax,time,pltslc] = ...
2     DiagLoad(fln,exp,dat,dad,grd,itr,tst,flu,ddf,gdf,avg,slc,pst,...
3     LoadGridData,DiagDebug,GridSuffix,ZcordFile,Index,Dim,...
4     Grads,Year0Iter,Months);
5    
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     % be monthly mean, and instantaneous data (dat=='Int'). There is also an
20     % 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     elseif ismember(dat,{'','Non','None'})
29     filesuffix = ''; mnchandle = '';
30     else
31     error(['DataType not recognized: ',datatype]);
32     end
33    
34     % Time stuff (itr = 0 is start of new 360 day year).
35     if isempty(Months)
36     months = (itr-Year0Iter)./(30.*24.*3600./tst);
37     else
38     months = Months;
39     end
40     time = months/12;
41    
42     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43     % Read in data %
44     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45    
46     % Open grads files
47     if ~isequal(Grads,0)
48     if DiagDebug, disp([' Debug -- Loading GRADS field.']); end
49     [data,xax,yax,zax,months,time,dim] = DiagLoadGradsData(Grads,dad,fln);
50     if DiagDebug, disp([' Debug -- ''data'' size after load: ',mat2str(size(data))]); end
51     data = DiagAverage(data,fln,avg,months,ddf,dim);
52     if DiagDebug, disp([' Debug -- ''data'' size after averaging: ',mat2str(size(data))]); end
53     if isequal(slc,'Sur')
54     data = data';
55     pltslc='lonlat';
56     else
57     error('Loading grads data only set to handle surface plots!');
58     end
59    
60    
61     % Load AIM physics data. The data is all stuffed into one file called
62     % 'aimPhy' which currently only comes in MDS format. When loading a field,
63     % call the variable name that is used in the actual MITgcm code. These
64     % variables are given in 'DiagGenParam'. Here, we check that to see if the
65     % parameter is an AIM parameter, then load the 'aimPhy' file and extract
66     % the data with the proper index as specified in 'DiagGenParam'. For total
67     % precipitation ('TotP') and evaporation minus precipitation ('EmP') some
68     % simple calculations of the AIM fields are performed and the U and V
69     % stresses (on atmosphere) must converted to a lat-lon grid.
70     elseif ismember(fln,aimparameters)
71     if isequal(ddf,'MNC')
72     error('Aim physics reader not set to hande ''MNC'' data format.')
73     end
74     data = DiagLoad_Local('aimPhy',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
75     if isequal(fln,'TotP')
76     data = data(:,:,PRECNV,:) ...
77     + data(:,:,PRECLS,:);
78     elseif isequal(fln,'EmP')
79     data = data(:,:,EVAP ,:) ...
80     - data(:,:,PRECNV,:) ...
81     - data(:,:,PRECLS,:);
82     elseif ismember(fln,{'USTR','VSTR'})
83     [XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
84     DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
85     eval(['USTRdata = data(:,:,',num2str(USTR),',:);']);
86     eval(['VSTRdata = data(:,:,',num2str(VSTR),',:);']);
87     [USTRdata,VSTRdata]=uvcube2latlon(XC,YC,USTRdata,VSTRdata,XL,YL);
88     if isequal(fln,'USTR'), data = USTRdata;
89     elseif isequal(fln,'VSTR'), data = VSTRdata; end
90     else
91     eval(['data = data(:,:,',fln,',:);']);
92     end
93     data = squeeze(data(:,:,1,:));
94     data = DiagAverage(data,fln,avg,months,ddf,Dim);
95     [data,xax,yax,pltslc] = ...
96     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...
97     ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
98    
99     % Load coupled fields. Like AIM, the data is all stuffed into one file
100     % called either 'cplFld.*' (MDS) or 'cpl_tave.*' (MNC). Variables within
101     % the MDS files given in 'DiagGenParam', accessed directly with MNC. TX
102     % and TY (stress on ocean) must be converted to lat-lon grid.
103     elseif ismember(fln,cplparameters)
104     if ismember(fln,{'TX','TY'})
105     [XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
106     DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
107     if isequal(ddf,'MNC')
108     TXdata = DiagLoad_Local('TX',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
109     TYdata = DiagLoad_Local('TY',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
110     elseif isequal(ddf,'MDS')
111     data = DiagLoad_Local('cplFld',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
112     eval(['TXdata = squeeze(data(:,:,',num2str(TX),',:));']);
113     eval(['TYdata = squeeze(data(:,:,',num2str(TY),',:));']);
114     end
115     [TXdata,TYdata]=uvcube2latlon(XC,YC,TXdata,TYdata,XL,YL);
116     if isequal(fln,'TX'), data = TXdata;
117     elseif isequal(fln,'TY'), data = TYdata; end
118     else
119     if isequal(ddf,'MNC')
120     data = DiagLoad_Local(fln,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(['data = squeeze(data(:,:,',fln,',:));']);
124     end
125     if isequal(fln,'SLP'), data = data ./ 100; end
126     end
127     data = DiagAverage(data,fln,avg,months,ddf,Dim);
128     [data,xax,yax,pltslc] = ...
129     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...
130     ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
131    
132     % Load ice parameters. Note that as the ice variable names appear to be in
133     % constant flux, this part quickly becomes out of date.
134     elseif ismember(fln,iceparameters)
135     if isequal(dat,'Ins')
136     data = DiagLoad_Local([fln,'-T'],dat,dad,grd,itr,ddf,filesuffix,mnchandle);
137     else
138     data = DiagLoad_Local('thSIce_',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
139     eval(['data = data(:,:,',fln,',:);']);
140     data = squeeze(data(:,:,1,:));
141     end
142     data = DiagAverage(data,fln,avg,months,ddf,Dim);
143     [data,xax,yax,pltslc] = ...
144     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,...
145     ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
146    
147     % Load generic fields where no modifications are needed.
148     elseif ismember(fln,{'S','T','Temp','aim_RH','phiHyd','Conv'})
149     data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
150     data = DiagAverage(data,fln,avg,months,ddf,Dim);
151     [data,xax,yax,pltslc] = ...
152     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
153     flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
154    
155     % Read vertical velocities, eta. Convert from [P] to [hPa] in atmosphere.
156     elseif ismember(fln,{'W','wVel','Eta','ETA'})
157     data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
158     if isequal(flu,'A'), data = data./100; end
159     data = DiagAverage(data,fln,avg,months,ddf,Dim);
160     [data,xax,yax,pltslc] = ...
161     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
162     flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
163    
164     % Load horizontal velocities, convert to lat-lon grid.
165     elseif ismember(fln,{'U','V','uVel','vVel'})
166     if isequal(dat,'Tav')
167     U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
168     V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
169     else isequal(dat,'Int')
170     U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
171     V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
172     end
173     if isequal(ddf,'MNC')
174     U = U(1:end-1,:,:,:);
175     V = V(:,1:end-1,:,:);
176     end
177     U = DiagAverage(U,fln,avg,months,ddf,Dim);
178     V = DiagAverage(V,fln,avg,months,ddf,Dim);
179     [XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
180     DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
181     [U,V]=uvcube2latlon(XC,YC,U,V,XL,YL);
182     if ismember(fln,{'U','uVel'})
183     data = U;
184     elseif ismember(fln,{'V','vVel'})
185     data = V;
186     end
187     [data,xax,yax,pltslc] = ...
188     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
189     flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
190    
191     % Meridional overturning.
192     elseif ismember(fln,{'Bol','Psi','Res'})
193     [XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
194     DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
195     if isequal(flu,'A'), delM = - 1000*drF ./ g;
196     elseif isequal(flu,'O'), delM = drF; end
197     if ismember(fln,{'Bol','Res'})
198     kwx = DiagLoad_Local('GM_Kwx-T',dat,dad,grd,itr,'MDS','',mnchandle);
199     kwy = DiagLoad_Local('GM_Kwy-T',dat,dad,grd,itr,'MDS','',mnchandle);
200     kwx = DiagAverage(kwx,fln,avg,months,ddf,Dim);
201     kwy = DiagAverage(kwy,fln,avg,months,ddf,Dim);
202     [ub,vb]=calc_Bolus_CS(RAC,dxC,dyC,dxG,dyG,drC ,kwx,kwy,HFacW,HFacS);
203     [psibol,mskG,ylat]=calc_PsiCube(delM,ub.*HFacW,vb.*HFacS,dxG,dyG);
204     end
205     if ismember(fln,{'Psi','Res'})
206     if isequal(dat,'Tav')
207     U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
208     V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
209     else isequal(dat,'Int')
210     U = DiagLoad_Local('U',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
211     V = DiagLoad_Local('V',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
212     end
213     if isequal(ddf,'MNC')
214     U = U(1:end-1,:,:,:);
215     V = V(:,1:end-1,:,:);
216     end
217     U = DiagAverage(U,fln,avg,months,ddf,Dim);
218     V = DiagAverage(V,fln,avg,months,ddf,Dim);
219     [psiavg,mskG,ylat]=calc_PsiCube(delM,U.*HFacW,V.*HFacS,dxG,dyG);
220     end
221     if isequal(fln,'Psi')
222     data = psiavg(2:end-1,:)';
223     elseif isequal(fln,'Bol')
224     data = psibol(2:end-1,:)';
225     elseif isequal(fln,'Res')
226     data = psibol(2:end-1,:)' + psiavg(2:end-1,:)';
227     end
228     if isequal(flu,'A'), data = data./1e6; end
229     xax = ylat; yax = ZG; pltslc = 'lathgt';
230     eval(['fixedrange = ',fln,'range',flu,';']);
231     datarange = [min(data(:)),max(data(:))];
232     if datarange(1) < fixedrange(1) | ...
233     datarange(2) > fixedrange(2)
234     disp(['***Warning*** Value out of range for ',fln]);
235     disp([' Calculated range: ',mat2str(datarange)]);
236     disp([' Given range: ',mat2str(fixedrange)]);
237     end
238    
239     % Read in variabilities. Here we must make some simple calculations. For
240     % all the variablilty fields there is are (field)^2 fields from which we
241     % can calculate the desired variabilities: sqrt(field*field-(field)^2).
242     % Note that 'ETAstd' is scaled from [P] to [hPa]. For 'KEpri', we should
243     % use the values on either side of a grid, since it will be plotted from
244     % the grid center, but here we just use from one edge out of convenience.
245     elseif isequal(fln,'ETAstd')
246     ETA = DiagLoad_Local('ETA' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle)./100 ;
247     ETA2 = DiagLoad_Local('Eta2',dat,dad,grd,itr,ddf,filesuffix,mnchandle)./10000;
248     data = sqrt(abs(ETA.*ETA-ETA2));
249     data = DiagAverage(data,fln,avg,months,ddf,Dim);
250     [data,xax,yax,pltslc] = ...
251     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
252     flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
253     elseif isequal(fln,'Tstd')
254     T = DiagLoad_Local('T' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
255     TT = DiagLoad_Local('TT',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
256     data = sqrt(abs(T.*T-TT));
257     data = DiagAverage(data,fln,avg,months,ddf,Dim);
258     [data,xax,yax,pltslc] = ...
259     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
260     flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
261     elseif isequal(fln,'KEpri')
262     U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
263     V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
264     UU = DiagLoad_Local('UU' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
265     VV = DiagLoad_Local('VV' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
266     if isequal(ddf,'MNC')
267     U = .5*( U(2:end,:,:,:) + U(1:end-1,:,:,:) );
268     V = .5*( V(:,2:end,:,:) + V(:,1:end-1,:,:) );
269     %UU = .5*( UU(2:end,:,:,:) + UU(1:end-1,:,:,:) ); % UU and VV might be incorrect.
270     %VV = .5*( VV(:,2:end,:,:) + VV(:,1:end-1,:,:) );
271     end
272     U = DiagAverage(U,fln,avg,months,ddf,Dim);
273     V = DiagAverage(V,fln,avg,months,ddf,Dim);
274     UU = DiagAverage(UU,fln,avg,months,ddf,Dim);
275     VV = DiagAverage(VV,fln,avg,months,ddf,Dim);
276     data = sqrt(abs((U.*U + V.*V) - (UU + VV)));
277     [data,xax,yax,pltslc] = ...
278     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
279     flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
280    
281     % If the field is not recognized, try best to try and open it.
282     else
283     if DiagDebug, disp([' Debug -- Unknown field, attempting to load.']); end
284     data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
285     if DiagDebug, disp([' Debug -- ''data'' size after load: ',mat2str(size(data))]); end
286     if ~isequal(Index,0)
287     if isequal(Dim,2), data = squeeze(data(:,:,Index,:));
288     elseif isequal(Dim,3), data = squeeze(data(:,:,:,Index,:)); end
289     end
290     if DiagDebug, disp([' Debug -- ''data'' size after indexing: ',mat2str(size(data))]); end
291     data = DiagAverage(data,fln,avg,months,ddf,Dim);
292     if DiagDebug, disp([' Debug -- ''data'' size after average: ',mat2str(size(data))]); end
293     [data,xax,yax,pltslc] = ...
294     DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,...
295     flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
296     if DiagDebug, disp([' Debug -- ''data'' size after slice: ',mat2str(size(data))]); end
297     if DiagDebug, disp([' Debug -- ''xax'' size: ',mat2str(size(xax))]); end
298     if DiagDebug, disp([' Debug -- ''yax'' size: ',mat2str(size(yax))]); end
299     if DiagDebug, disp([' Debug -- ''pltslc'': ',pltslc]); end
300     end
301    
302    
303     %-------------------------------------------------------------------------%
304     % Local functions %
305     %-------------------------------------------------------------------------%
306    
307     % Thus is merely a local function that calls the load data functions
308     % according to the DataFormat for the data specified by dfm.
309     function data = DiagLoad_Local(fln,dat,dad,grd,itr,dfm,filesuffix,mnchandle)
310    
311     if isequal(dfm,'MDS')
312     data = rdmds([dad,'/',fln,filesuffix],itr);
313     elseif isequal(dfm,'MNC')
314     %[dad,mnchandle]
315     iter = rdmnc_mod([dad,mnchandle],'T');
316     iter = iter.T;
317     [dummy,indecies] = ismember(itr,iter);
318     data = rdmnc_mod([dad,mnchandle],[fln,filesuffix],'T',indecies);
319     iter = data.T;
320     if ~isequal(itr,iter'), error('Missing iterations in data!'); end
321     eval(['data = data.',fln,filesuffix,';']);
322     else
323     error(['Unrecognized data type: ',dfm]);
324     end

  ViewVC Help
Powered by ViewVC 1.1.22