/[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.4 - (hide annotations) (download)
Wed Feb 2 15:38:21 2005 UTC (20 years, 5 months ago) by enderton
Branch: MAIN
Changes since 1.3: +49 -14 lines
 o Added plotting options for actual temperature 'ActT' and moist potential temperature 'MoiPTc', though the later is calculated in a crude way.  Also changed 'exp' variables internally to 'trl' to avoid overriding exp() function.  Hopefully I got them all.

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 enderton 1.1 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 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     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 molod 1.2 [data,xax,yax,zax,months,time,dim] = DiagLoadGradsData(Grads,dad,fln,DiagDebug);
50 enderton 1.1 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 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
97 enderton 1.1 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 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
130 enderton 1.1 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 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,...
145 enderton 1.1 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 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
153 enderton 1.1 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 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
162 enderton 1.1 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 molod 1.3 else isequal(dat,'Ins')
170 enderton 1.1 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 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
189 enderton 1.1 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 molod 1.3 else isequal(dat,'Ins')
210 enderton 1.1 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 enderton 1.4
239    
240     % Actual temperature and moist potential temperature. Actual temperature
241     % is calculated as T=Theta*(p/po)^(R/cp). Moist potential temperature is
242     % calculated as Theta_e=Theta*e^(L*q*/cp*T) where q*=(R/Rv)*(es/p) and
243     % es=Ae^(beta*T).
244     elseif ismember(fln,{'ActT','MoiPTc'})
245     if isequal('flu','O'),
246     error('Calculation may only be done for atmosphere!');
247     end
248     [XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,...
249     HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
250     DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
251     if isequal(ddf,'MDS')
252     theta = DiagLoad_Local('T',dat,dad,grd,itr,ddf,...
253     filesuffix,mnchandle);
254     elseif isequal(ddf,'MNC')
255     theta = DiagLoad_Local('Temp',dat,dad,grd,itr,...
256     ddf,filesuffix,mnchandle);
257     end
258     theta = DiagAverage(theta,fln,avg,months,ddf,Dim);
259     if ismember(fln,{'ActT','MoiPTc'})
260     pres = NaN.*zeros(size(theta));
261     for iz=1:length(ZC), pres(:,:,iz)=ZC(iz); end
262     temp=theta.*(pres./presrefA).^(RdA/cpA);
263     if isequal(fln,'ActT'), data=temp; end
264     end
265     if isequal(fln,'MoiPTc')
266     es=A_CC.*exp(Beta_CC.*(temp-K2C));
267     qstar=(RdA/RvA).*(es./pres);
268     data=theta.*exp(LHvapA.*qstar./cpA./temp);
269     end
270     [data,xax,yax,pltslc] = ...
271     DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,...
272     avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
273    
274    
275 enderton 1.1 % Read in variabilities. Here we must make some simple calculations. For
276     % all the variablilty fields there is are (field)^2 fields from which we
277     % can calculate the desired variabilities: sqrt(field*field-(field)^2).
278     % Note that 'ETAstd' is scaled from [P] to [hPa]. For 'KEpri', we should
279     % use the values on either side of a grid, since it will be plotted from
280     % the grid center, but here we just use from one edge out of convenience.
281     elseif isequal(fln,'ETAstd')
282     ETA = DiagLoad_Local('ETA' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle)./100 ;
283     ETA2 = DiagLoad_Local('Eta2',dat,dad,grd,itr,ddf,filesuffix,mnchandle)./10000;
284     data = sqrt(abs(ETA.*ETA-ETA2));
285     data = DiagAverage(data,fln,avg,months,ddf,Dim);
286     [data,xax,yax,pltslc] = ...
287 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
288 enderton 1.1 flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
289     elseif isequal(fln,'Tstd')
290     T = DiagLoad_Local('T' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
291     TT = DiagLoad_Local('TT',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
292     data = sqrt(abs(T.*T-TT));
293     data = DiagAverage(data,fln,avg,months,ddf,Dim);
294     [data,xax,yax,pltslc] = ...
295 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
296 enderton 1.1 flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
297     elseif isequal(fln,'KEpri')
298     U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
299     V = DiagLoad_Local('vVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle);
300     UU = DiagLoad_Local('UU' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
301     VV = DiagLoad_Local('VV' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
302     if isequal(ddf,'MNC')
303     U = .5*( U(2:end,:,:,:) + U(1:end-1,:,:,:) );
304     V = .5*( V(:,2:end,:,:) + V(:,1:end-1,:,:) );
305     %UU = .5*( UU(2:end,:,:,:) + UU(1:end-1,:,:,:) ); % UU and VV might be incorrect.
306     %VV = .5*( VV(:,2:end,:,:) + VV(:,1:end-1,:,:) );
307     end
308     U = DiagAverage(U,fln,avg,months,ddf,Dim);
309     V = DiagAverage(V,fln,avg,months,ddf,Dim);
310     UU = DiagAverage(UU,fln,avg,months,ddf,Dim);
311     VV = DiagAverage(VV,fln,avg,months,ddf,Dim);
312     data = sqrt(abs((U.*U + V.*V) - (UU + VV)));
313     [data,xax,yax,pltslc] = ...
314 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
315 enderton 1.1 flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
316    
317     % If the field is not recognized, try best to try and open it.
318     else
319 enderton 1.4 disp([' Unknown field, attempting to load.']);
320 enderton 1.1 data = DiagLoad_Local(fln,dat,dad,grd,itr,ddf,filesuffix,mnchandle);
321     if DiagDebug, disp([' Debug -- ''data'' size after load: ',mat2str(size(data))]); end
322     if ~isequal(Index,0)
323     if isequal(Dim,2), data = squeeze(data(:,:,Index,:));
324     elseif isequal(Dim,3), data = squeeze(data(:,:,:,Index,:)); end
325     end
326     if DiagDebug, disp([' Debug -- ''data'' size after indexing: ',mat2str(size(data))]); end
327     data = DiagAverage(data,fln,avg,months,ddf,Dim);
328     if DiagDebug, disp([' Debug -- ''data'' size after average: ',mat2str(size(data))]); end
329     [data,xax,yax,pltslc] = ...
330 enderton 1.4 DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,...
331 enderton 1.1 flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
332     if DiagDebug, disp([' Debug -- ''data'' size after slice: ',mat2str(size(data))]); end
333     if DiagDebug, disp([' Debug -- ''xax'' size: ',mat2str(size(xax))]); end
334     if DiagDebug, disp([' Debug -- ''yax'' size: ',mat2str(size(yax))]); end
335     if DiagDebug, disp([' Debug -- ''pltslc'': ',pltslc]); end
336     end
337    
338    
339     %-------------------------------------------------------------------------%
340     % Local functions %
341     %-------------------------------------------------------------------------%
342    
343     % Thus is merely a local function that calls the load data functions
344     % according to the DataFormat for the data specified by dfm.
345     function data = DiagLoad_Local(fln,dat,dad,grd,itr,dfm,filesuffix,mnchandle)
346    
347     if isequal(dfm,'MDS')
348     data = rdmds([dad,'/',fln,filesuffix],itr);
349     elseif isequal(dfm,'MNC')
350     iter = rdmnc_mod([dad,mnchandle],'T');
351     iter = iter.T;
352     [dummy,indecies] = ismember(itr,iter);
353     data = rdmnc_mod([dad,mnchandle],[fln,filesuffix],'T',indecies);
354     iter = data.T;
355     if ~isequal(itr,iter'), error('Missing iterations in data!'); end
356     eval(['data = data.',fln,filesuffix,';']);
357     else
358     error(['Unrecognized data type: ',dfm]);
359     end

  ViewVC Help
Powered by ViewVC 1.1.22