--- MITgcm_contrib/enderton/Diagnostics/DiagLoad.m 2005/01/31 15:43:26 1.1 +++ MITgcm_contrib/enderton/Diagnostics/DiagLoad.m 2005/02/02 16:20:12 1.5 @@ -1,5 +1,5 @@ function [data,xax,yax,time,pltslc] = ... - 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,... LoadGridData,DiagDebug,GridSuffix,ZcordFile,Index,Dim,... Grads,Year0Iter,Months); @@ -16,7 +16,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Currently the function can handle time averages (dat=='Tav'), which must -% be monthly mean, and instantaneous data (dat=='Int'). There is also an +% 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') @@ -46,7 +46,7 @@ % 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); + [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 @@ -93,7 +93,7 @@ data = squeeze(data(:,:,1,:)); data = DiagAverage(data,fln,avg,months,ddf,Dim); [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,... ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); % Load coupled fields. Like AIM, the data is all stuffed into one file @@ -126,7 +126,7 @@ end data = DiagAverage(data,fln,avg,months,ddf,Dim); [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,... ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); % Load ice parameters. Note that as the ice variable names appear to be in @@ -141,7 +141,7 @@ end data = DiagAverage(data,fln,avg,months,ddf,Dim); [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,flu,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,... ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); % Load generic fields where no modifications are needed. @@ -149,7 +149,7 @@ 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,exp,dat,dad,grd,itr,tst,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,... flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); % Read vertical velocities, eta. Convert from [P] to [hPa] in atmosphere. @@ -158,7 +158,7 @@ if isequal(flu,'A'), data = data./100; end data = DiagAverage(data,fln,avg,months,ddf,Dim); [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,... flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); % Load horizontal velocities, convert to lat-lon grid. @@ -166,7 +166,7 @@ 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,'Int') + 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 @@ -185,7 +185,7 @@ data = V; end [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,... flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); % Meridional overturning. @@ -206,7 +206,7 @@ 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,'Int') + 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 @@ -235,7 +235,43 @@ 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); + + % 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). @@ -248,7 +284,7 @@ data = sqrt(abs(ETA.*ETA-ETA2)); data = DiagAverage(data,fln,avg,months,ddf,Dim); [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,... flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); elseif isequal(fln,'Tstd') T = DiagLoad_Local('T' ,dat,dad,grd,itr,ddf,filesuffix,mnchandle); @@ -256,7 +292,7 @@ data = sqrt(abs(T.*T-TT)); data = DiagAverage(data,fln,avg,months,ddf,Dim); [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,... flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); elseif isequal(fln,'KEpri') U = DiagLoad_Local('uVel',dat,dad,grd,itr,ddf,filesuffix,mnchandle); @@ -275,12 +311,12 @@ VV = DiagAverage(VV,fln,avg,months,ddf,Dim); data = sqrt(abs((U.*U + V.*V) - (UU + VV))); [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,... flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); % If the field is not recognized, try best to try and open it. else - if DiagDebug, disp([' Debug -- Unknown field, attempting to load.']); end + 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(Index,0) @@ -291,7 +327,7 @@ data = DiagAverage(data,fln,avg,months,ddf,Dim); if DiagDebug, disp([' Debug -- ''data'' size after average: ',mat2str(size(data))]); end [data,xax,yax,pltslc] = ... - DiagSlice(data,fln,exp,dat,dad,grd,itr,tst,... + DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,... flu,ddf,gdf,avg,slc,pst,LoadGridData,GridSuffix,ZcordFile); if DiagDebug, disp([' Debug -- ''data'' size after slice: ',mat2str(size(data))]); end if DiagDebug, disp([' Debug -- ''xax'' size: ',mat2str(size(xax))]); end @@ -311,7 +347,6 @@ if isequal(dfm,'MDS') data = rdmds([dad,'/',fln,filesuffix],itr); elseif isequal(dfm,'MNC') - %[dad,mnchandle] iter = rdmnc_mod([dad,mnchandle],'T'); iter = iter.T; [dummy,indecies] = ismember(itr,iter);