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

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

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


Revision 1.3 - (show annotations) (download)
Tue Feb 1 18:49:28 2005 UTC (20 years, 5 months ago) by molod
Branch: MAIN
Changes since 1.2: +3 -3 lines
Change all 'Int' to 'Ins' for instantaneous fields

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=='Ins'). 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,DiagDebug);
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,'Ins')
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,'Ins')
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