[Dart-dev] [5688] DART/branches/development: Support for unstructured grids (mpas_atm) and replacing many
nancy at ucar.edu
nancy at ucar.edu
Tue Apr 10 21:38:39 MDT 2012
Revision: 5688
Author: thoar
Date: 2012-04-10 21:38:37 -0600 (Tue, 10 Apr 2012)
Log Message:
-----------
Support for unstructured grids (mpas_atm) and replacing many
one-off (hard to support) functions with get_hyperslab().
Each model with a RunAllTests.m script exercises all of
the state-space diagnostic scripts.
The plot_total_err() script no longer plots the total error
and spread for each level separately, but performs a simple
average over all available levels. Then, a spatial weighting
to account for the area of each grid cell is applied, and the
square root is applied ... Root(weighted-mean(squared(error))).
For the 'higher-order' models, each variable is plotted separately
in its own figure window.
Modified Paths:
--------------
DART/branches/development/matlab/CheckModel.m
DART/branches/development/matlab/CheckModelCompatibility.m
DART/branches/development/matlab/GetBgridInfo.m
DART/branches/development/matlab/GetCamInfo.m
DART/branches/development/matlab/GetMITgcm_oceanInfo.m
DART/branches/development/matlab/GetNCindices.m
DART/branches/development/matlab/GetPe2lyrInfo.m
DART/branches/development/matlab/GetTIEGCMInfo.m
DART/branches/development/matlab/GetWRFInfo.m
DART/branches/development/matlab/ParseAlphaNumerics.m
DART/branches/development/matlab/PlotBins.m
DART/branches/development/matlab/PlotCEnsErrSpread.m
DART/branches/development/matlab/PlotCorrel.m
DART/branches/development/matlab/PlotEnsErrSpread.m
DART/branches/development/matlab/PlotEnsMeanTimeSeries.m
DART/branches/development/matlab/PlotEnsTimeSeries.m
DART/branches/development/matlab/PlotJeffCorrel.m
DART/branches/development/matlab/PlotPhaseSpace.m
DART/branches/development/matlab/PlotSawtooth.m
DART/branches/development/matlab/PlotTotalErr.m
DART/branches/development/matlab/PlotVarVarCorrel.m
DART/branches/development/matlab/SetCopyID.m
DART/branches/development/matlab/SetVariableID.m
DART/branches/development/matlab/SimpleMap.m
DART/branches/development/matlab/get_ens_series.m
DART/branches/development/matlab/get_state_copy.m
DART/branches/development/matlab/get_var_series.m
DART/branches/development/matlab/jeff_correl.m
DART/branches/development/matlab/plot_bins.m
DART/branches/development/matlab/plot_correl.m
DART/branches/development/matlab/plot_ens_err_spread.m
DART/branches/development/matlab/plot_ens_mean_time_series.m
DART/branches/development/matlab/plot_ens_time_series.m
DART/branches/development/matlab/plot_jeff_correl.m
DART/branches/development/matlab/plot_phase_space.m
DART/branches/development/matlab/plot_reg_factor.m
DART/branches/development/matlab/plot_sawtooth.m
DART/branches/development/matlab/plot_smoother_err.m
DART/branches/development/matlab/plot_total_err.m
DART/branches/development/matlab/plot_var_var_correl.m
DART/branches/development/matlab/total_err.m
DART/branches/development/models/9var/matlab/RunAllTests.m
DART/branches/development/models/bgrid_solo/matlab/RunAllTests.m
DART/branches/development/models/forced_lorenz_96/matlab/RunAllTests.m
DART/branches/development/models/ikeda/matlab/RunAllTests.m
DART/branches/development/models/lorenz_04/matlab/RunAllTests.m
DART/branches/development/models/lorenz_63/matlab/RunAllTests.m
DART/branches/development/models/lorenz_84/matlab/RunAllTests.m
DART/branches/development/models/lorenz_96/matlab/RunAllTests.m
DART/branches/development/models/lorenz_96_2scale/matlab/RunAllTests.m
DART/branches/development/models/simple_advection/matlab/RunAllTests.m
Added Paths:
-----------
DART/branches/development/matlab/BgridTotalError.m
DART/branches/development/matlab/CAMTotalError.m
DART/branches/development/matlab/GetMPAS_ATMInfo.m
DART/branches/development/matlab/WRFTotalError.m
DART/branches/development/matlab/continents.m
DART/branches/development/matlab/get_ensemble_indices.m
DART/branches/development/matlab/get_hyperslab.m
DART/branches/development/models/cam/matlab/RunAllTests.m
DART/branches/development/models/pe2lyr/matlab/
DART/branches/development/models/pe2lyr/matlab/RunAllTests.m
DART/branches/development/models/wrf/matlab/RunAllTests.m
-------------- next part --------------
Added: DART/branches/development/matlab/BgridTotalError.m
===================================================================
--- DART/branches/development/matlab/BgridTotalError.m (rev 0)
+++ DART/branches/development/matlab/BgridTotalError.m 2012-04-11 03:38:37 UTC (rev 5688)
@@ -0,0 +1,148 @@
+function BgridTotalError( pinfo )
+%% -------------------------------------------------------------------
+% Plot the total area-weighted error for each variable.
+%---------------------------------------------------------------------
+
+%% DART software - Copyright 2004 - 2011 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% <next few lines under version control, do not edit>
+% $URL: https://proxy.subversion.ucar.edu/DAReS/DART/branches/mpas/matlab/BgridTotalError.m $
+% $Id: BgridTotalError.m 5655 2012-04-05 23:17:16Z thoar $
+% $Revision: 5655 $
+% $Date: 2012-04-05 17:17:16 -0600 (Thu, 05 Apr 2012) $
+
+% Get the indices for the true state, ensemble mean and spread
+% The metadata is queried to determine which "copy" is appropriate.
+truth_index = get_copy_index(pinfo.truth_file, 'true state');
+ens_mean_index = get_copy_index(pinfo.diagn_file, 'ensemble mean');
+ens_spread_index = get_copy_index(pinfo.diagn_file, 'ensemble spread');
+
+%----------------------------------------------------------------------
+%
+%----------------------------------------------------------------------
+
+for ivar=1:pinfo.num_state_vars,
+
+ fprintf('Processing %s ...\n', pinfo.vars{ivar} )
+
+ rmse = zeros(pinfo.time_series_length,1);
+ sprd = zeros(pinfo.time_series_length,1);
+ varunits = nc_attget(pinfo.truth_file, pinfo.vars{ivar}, 'units');
+
+ % determine what grid the variable lives on
+ % determine the number of levels
+
+ nlevels = 1;
+
+ varinfo = nc_getvarinfo(pinfo.diagn_file,pinfo.vars{ivar});
+
+ for idim = 1:length(varinfo.Dimension),
+ dimname = varinfo.Dimension{idim};
+ dimlength = varinfo.Size(idim);
+ switch lower(dimname)
+ case {'tmpj', 'velj'}
+ latitudes = nc_varget(pinfo.diagn_file, dimname);
+ case {'tmpi', 'veli'}
+ longitudes = nc_varget(pinfo.diagn_file, dimname);
+ case {'lev'}
+ nlevels = dimlength;
+ end
+ end
+
+ % Calculate weights for area-averaging.
+ weights = SphereWeights(latitudes, longitudes);
+
+ for itime=1:pinfo.time_series_length,
+
+ truth = get_hyperslab('fname',pinfo.truth_file, 'varname',pinfo.vars{ivar}, ...
+ 'copyindex',truth_index, 'timeindex',pinfo.truth_time(1)+itime-1);
+ ens = get_hyperslab('fname',pinfo.diagn_file, 'varname',pinfo.vars{ivar}, ...
+ 'copyindex',ens_mean_index, 'timeindex',pinfo.diagn_time(1)+itime-1);
+ spread = get_hyperslab('fname',pinfo.diagn_file, 'varname',pinfo.vars{ivar}, ...
+ 'copyindex',ens_spread_index, 'timeindex',pinfo.diagn_time(1)+itime-1);
+
+ %% Calculate the weighted mean squared error for each level.
+ % tensors come back [nlev,nlat,nlon] - or - [nlat,nlon]
+
+ sqerr = (truth - ens).^2;
+ sqsprd = spread .^2;
+
+ if (nlevels > 1) % take the mean over the first dimension
+ sqerr = squeeze(mean(sqerr ,1));
+ sqsprd = squeeze(mean(sqsprd,1));
+ end
+
+ %% Create the (weighted) mean squared error
+
+ ms_err = sum(sqerr(:) .* weights);
+ ms_spread = sum(sqsprd(:) .* weights);
+
+ %% Take the square root of the mean squared error
+ rmse(itime) = sqrt(ms_err);
+ sprd(itime) = sqrt(ms_spread);
+
+ end % loop over time
+
+ %-------------------------------------------------------------------
+ % Each variable in its own figure window
+ %-------------------------------------------------------------------
+ figure(ivar); clf;
+ plot(pinfo.time,rmse,'-', pinfo.time,sprd,'--')
+
+ s{1} = sprintf('time-mean Ensemble Mean error = %f', mean(rmse));
+ s{2} = sprintf('time-mean Ensemble Spread = %f', mean(sprd));
+
+ h = legend(s); legend(h,'boxoff')
+ grid on;
+ xdates(pinfo.time)
+ ylabel(sprintf('global-area-weighted rmse (%s)',varunits))
+ s1 = sprintf('%s %s Ensemble Mean', pinfo.model,pinfo.vars{ivar});
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
+
+end % loop around variables
+
+clear truth ens spread err XY_spread
+
+
+
+
+function weights = SphereWeights(lats,lons)
+%% SphereWeights creates [nlat*nlon,1] matrix of weights based on latitude
+%
+% lats,lons must be 1D arrays (in degrees)
+
+nlats = length(lats);
+nlons = length(lons);
+
+if ( numel(lats) ~= nlats )
+ disp('latitude array is of higher dimension than anticipated.')
+ error('Must be a vector.')
+end
+if ( numel(lons) ~= nlons )
+ disp('longitude array is of higher dimension than anticipated.')
+ error('Must be a vector.')
+end
+
+rads = zeros(nlats,1); % Ensure lats is a column vector,
+rads(:) = pi*lats/180.0; % and convert to radians.
+wts = cos( rads ) * ones(1,nlons); % Results in a [nlat-x-nlon] matrix.
+wts = wts ./ sum(wts(:)); % Normalize to unity.
+weights = wts(:);
+
+
+
+
+function xdates(dates)
+if (length(get(gca,'XTick')) > 6)
+ datetick('x','mm.dd.HH','keeplimits'); % 'mm/dd'
+ monstr = datestr(dates(1),31);
+ xlabelstring = sprintf('month/day/HH - %s start',monstr);
+else
+ datetick('x',31,'keeplimits'); %'yyyy-mm-dd HH:MM:SS'
+ monstr = datestr(dates(1),31);
+ xlabelstring = sprintf('%s start',monstr);
+end
+xlabel(xlabelstring)
+
Added: DART/branches/development/matlab/CAMTotalError.m
===================================================================
--- DART/branches/development/matlab/CAMTotalError.m (rev 0)
+++ DART/branches/development/matlab/CAMTotalError.m 2012-04-11 03:38:37 UTC (rev 5688)
@@ -0,0 +1,120 @@
+function CAMTotalError( pinfo )
+%% -------------------------------------------------------------------
+% Plot the total area-weighted error for each variable.
+%---------------------------------------------------------------------
+
+%% DART software - Copyright 2004 - 2011 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% <next few lines under version control, do not edit>
+% $URL: https://proxy.subversion.ucar.edu/DAReS/DART/branches/mpas/matlab/CAMTotalError.m $
+% $Id: CAMTotalError.m 5655 2012-04-05 23:17:16Z thoar $
+% $Revision: 5655 $
+% $Date: 2012-04-05 17:17:16 -0600 (Thu, 05 Apr 2012) $
+
+% Since the models are "compatible", get the info from either one.
+lons = nc_varget(pinfo.truth_file, 'lon');
+gw = nc_varget(pinfo.truth_file, 'gw');
+num_lons = length(lons);
+
+% make a matrix of weights for each horizontal slice
+% ensure the gaussian weights sum to unity.
+[~,weights] = meshgrid(ones(1,num_lons),gw);
+weights = weights/sum(weights(:));
+wts = weights(:);
+
+% Get the indices for the true state, ensemble mean and spread
+% The metadata is queried to determine which "copy" is appropriate.
+truth_index = get_copy_index(pinfo.truth_file, 'true state');
+ens_mean_index = get_copy_index(pinfo.diagn_file, 'ensemble mean');
+ens_spread_index = get_copy_index(pinfo.diagn_file, 'ensemble spread');
+
+%----------------------------------------------------------------------
+%
+%----------------------------------------------------------------------
+
+for ivar=1:pinfo.num_state_vars,
+
+ varname = pinfo.vars{ivar};
+
+ rmse = zeros(pinfo.time_series_length,1);
+ sprd = zeros(pinfo.time_series_length,1);
+
+ for itime=1:pinfo.time_series_length,
+
+ fprintf('Processing %s timestep %d of %d ...\n', ...
+ varname, itime, pinfo.time_series_length)
+
+ truth = get_hyperslab('fname',pinfo.truth_file, 'varname',varname, ...
+ 'copyindex',truth_index, 'timeindex',pinfo.truth_time(1)+itime-1);
+ ens = get_hyperslab('fname',pinfo.diagn_file, 'varname',varname, ...
+ 'copyindex',ens_mean_index, 'timeindex',pinfo.diagn_time(1)+itime-1);
+ spread = get_hyperslab('fname',pinfo.diagn_file, 'varname',varname, ...
+ 'copyindex',ens_spread_index, 'timeindex',pinfo.diagn_time(1)+itime-1);
+
+ if (length(size(truth)) == 2)
+ nlev = 1;
+ elseif (length(size(truth)) == 3)
+ nlev = size(truth,3);
+ else
+ error('Dang, this cannot happen in CAM.')
+ end
+
+ %% Calculate the weighted mean squared error for each level.
+
+ msqe_Z = zeros(nlev,1);
+ sprd_Z = zeros(nlev,1);
+
+ for ilevel=1:nlev,
+
+ slabS2E = (truth(:,:,ilevel) - ens(:,:,ilevel)).^2; % OK even if 2D iff ilevel = 1
+ XY_err = sum(slabS2E(:) .* wts);
+ slabS2E = spread(:,:,ilevel).^2;
+ XY_spread = sum(slabS2E(:) .* wts);
+
+ msqe_Z(ilevel) = XY_err;
+ sprd_Z(ilevel) = XY_spread;
+
+ end % loop over levels
+
+ %% Take the square root of the mean of all levels
+ rmse(itime) = sqrt(mean(msqe_Z));
+ sprd(itime) = sqrt(mean(sprd_Z));
+
+ end % loop over time
+
+ %-------------------------------------------------------------------
+ % Each variable in its own figure window
+ %-------------------------------------------------------------------
+ figure(ivar); clf;
+ varunits = nc_attget(pinfo.truth_file, pinfo.vars{ivar}, 'units');
+
+ plot(pinfo.time,rmse,'-', pinfo.time,sprd,'--')
+
+ s{1} = sprintf('time-mean Ensemble Mean error = %f', mean(rmse));
+ s{2} = sprintf('time-mean Ensemble Spread = %f', mean(sprd));
+
+ h = legend(s); legend(h,'boxoff')
+ grid on;
+ xdates(pinfo.time)
+ ylabel(sprintf('global-area-weighted rmse (%s)',varunits))
+ s1 = sprintf('%s %s Ensemble Mean', pinfo.model,pinfo.vars{ivar});
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
+
+end % loop around variables
+
+clear truth ens spread err XY_spread
+
+
+function xdates(dates)
+if (length(get(gca,'XTick')) > 6)
+ datetick('x','mm.dd.HH','keeplimits'); % 'mm/dd'
+ monstr = datestr(dates(1),31);
+ xlabelstring = sprintf('month/day/HH - %s start',monstr);
+else
+ datetick('x',31,'keeplimits'); %'yyyy-mm-dd HH:MM:SS'
+ monstr = datestr(dates(1),31);
+ xlabelstring = sprintf('%s start',monstr);
+end
+xlabel(xlabelstring)
Modified: DART/branches/development/matlab/CheckModel.m
===================================================================
--- DART/branches/development/matlab/CheckModel.m 2012-04-11 03:29:46 UTC (rev 5687)
+++ DART/branches/development/matlab/CheckModel.m 2012-04-11 03:38:37 UTC (rev 5688)
@@ -1,4 +1,4 @@
-function vars = CheckModel(fname);
+function vars = CheckModel(fname)
%% CheckModel tries to ensure that a netcdf file has what we expect.
%
% vars is a structure containing a minimal amount of metadata about the netCDF file.
@@ -20,17 +20,22 @@
if ( exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
% Get some information from the file
-model = nc_attget(fname,nc_global,'model');
-
+model = nc_attget(fname,nc_global,'model');
num_copies = dim_length(fname,'copy'); % determine # of ensemble members
-num_times = dim_length(fname,'time'); % determine # of output times
+[ens_size, ens_indices] = get_ensemble_indices(fname);
+times = nc_varget(fname,'time');
+timeunits = nc_attget(fname,'time','units');
+timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+dates = times + timeorigin;
+num_times = length(dates);
+clear times timeunits timebase timeorigin
+
if (isempty(model))
error('%s has no ''model'' global attribute.',fname)
end
-copy = nc_varget(fname,'copy');
-
switch lower(model)
case {'9var','lorenz_63','lorenz_84','ikeda'}
@@ -44,12 +49,13 @@
vars = struct('model',model, ...
'def_var','state', ...
'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
+ 'num_copies',num_copies, ...
+ 'num_ens_members',ens_size, ...
+ 'ensemble_indices',ens_indices, ...
+ 'time',dates, ...
'time_series_length',num_times, ...
'min_state_var',min(StateVariable), ...
'max_state_var',max(StateVariable), ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy), ...
'def_state_vars',def_state_vars);
vars.fname = fname;
@@ -67,12 +73,13 @@
vars = struct('model',model, ...
'def_var','state', ...
'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
+ 'num_copies',num_copies, ...
+ 'num_ens_members',ens_size, ...
+ 'ensemble_indices',ens_indices, ...
+ 'time',dates, ...
'time_series_length',num_times, ...
'min_state_var',min(StateVariable), ...
'max_state_var',max(StateVariable), ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy), ...
'def_state_vars',def_state_vars);
vars.fname = fname;
@@ -87,8 +94,7 @@
time_step_seconds = nc_attget(fname, nc_global, 'model_time_step_seconds');
num_model_vars = nc_attget(fname, nc_global, 'model_num_state_vars');
- num_vars = dim_length(fname,'StateVariable'); % determine # of state varbls
- StateVariable = nc_varget(fname,'StateVariable');
+ num_vars = dim_length(fname,'StateVariable'); % determine # of state varbls
% The only trick is to pick an equally-spaced subset of state
% variables for the default.
@@ -101,14 +107,20 @@
'num_state_vars',num_vars, ...
'num_model_vars',num_model_vars, ...
'num_force_vars',num_vars - num_model_vars, ...
- 'num_ens_members',num_copies, ...
+ 'num_copies',num_copies, ...
+ 'num_ens_members',ens_size, ...
+ 'ensemble_indices',ens_indices, ...
+ 'time',dates, ...
'time_series_length',num_times, ...
'min_state_var', 1, 'max_state_var', num_vars, ...
'min_model_var', 1, 'max_model_var', num_model_vars, ...
'min_force_var', 1, 'max_force_var', num_vars - num_model_vars, ...
- 'min_ens_mem',min(copy), 'max_ens_mem', max(copy), ...
'def_state_vars',def_state_vars, ...
- 'def_force_vars',def_force_vars);
+ 'def_force_vars',def_force_vars, ...
+ 'forcing',forcing, ...
+ 'delta_t',delta_t, ...
+ 'time_step_days', time_step_days, ...
+ 'time_step_seconds', time_step_seconds);
vars.fname = fname;
@@ -129,13 +141,16 @@
vars = struct('model',model, ...
'def_var','X', ...
'num_state_vars',num_X, ...
- 'num_ens_members',num_copies, ...
+ 'num_copies',num_copies, ...
+ 'num_ens_members',ens_size, ...
+ 'ensemble_indices',ens_indices, ...
+ 'time',dates, ...
'time_series_length',num_times, ...
'min_state_var',min(Xdim), 'max_state_var',max(Xdim), ...
'min_X_var', min(Xdim), 'max_X_var', max(Xdim), ...
'min_Y_var', min(Ydim), 'max_Y_var', max(Ydim), ...
- 'min_ens_mem', min(copy), 'max_ens_mem', max(copy), ...
- 'def_state_vars',def_X_inds);
+ 'def_state_vars',def_X_inds, ...
+ 'def_Y_inds', def_Y_inds);
vars.fname = fname;
@@ -148,160 +163,71 @@
varnames = {'state'};
def_inds = [1 13 27];
else
- varnames = {'concentration','source','wind', ...
- 'mean_source','source_phase'};
+ varnames = get_DARTvars(fname);
def_inds = round([1 , num_locs/3 , 2*num_locs/3]);
end
vars = struct('model' ,model, ...
'loc1d' ,loc1d, ...
- 'num_ens_members' ,num_copies, ...
- 'min_ens_mem' ,min(copy), ...
- 'max_ens_mem' ,max(copy), ...
+ 'num_copies' ,num_copies, ...
+ 'num_ens_members' ,ens_size, ...
+ 'ensemble_indices' ,ens_indices, ...
+ 'time' ,dates, ...
'time_series_length',num_times, ...
'model_size' ,length(varnames)*length(loc1d), ...
'def_var' ,varnames{1}, ...
'min_state_var' ,1, ...
'max_state_var' ,num_locs, ...
+ 'num_state_vars' ,num_locs, ...
'def_state_vars' ,def_inds, ...
'num_vars' ,length(varnames));
vars.vars = varnames;
vars.fname = fname;
- case 'fms_bgrid'
+ case 'wrf'
- % A more robust way would be to use the netcdf low-level ops:
- % bob = var(f); % bob is a cell array of ncvars
- % name(bob{1}) % is the variable name string
- % bob{1}(:) % is the value of the netcdf variable (no offset/scale)
+ % requires a 'domain' and 'bottom_top_d01' dimension.
+ % without both of these, it will fail in an ugly fashion.
- varnames = {'ps','t','u','v'};
- num_vars = length(varnames);
- nlevels = dim_length(fname,'lev'); % determine # of state variables
+ varnames = get_DARTvars(fname);
+ num_vars = length(varnames);
+ num_domains = dim_length(fname,'domain');
+ num_levels = dim_length(fname,'bottom_top_d01');
-% times = nc_varget(fname,'time');
-% TmpI = nc_varget(fname,'TmpI'); % longitude
-% TmpJ = nc_varget(fname,'TmpJ'); % latitude
-% levels = nc_varget(fname,'level');
-% VelI = nc_varget(fname,'VelI'); % longitude
-% VelJ = nc_varget(fname,'VelJ'); % latitude
-
vars = struct('model',model, ...
+ 'def_var',varnames{1}, ...
+ 'def_state_vars',[], ...
'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
+ 'num_copies',num_copies, ...
+ 'num_ens_members',ens_size, ...
+ 'ensemble_indices',ens_indices, ...
+ 'time',dates, ...
'time_series_length',num_times, ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy));
+ 'num_unstaggered_levels',num_levels, ...
+ 'num_domains',num_domains);
vars.vars = varnames;
vars.fname = fname;
+
+ case {'cam','tiegcm','fms_bgrid','pe2lyr','mitgcm_ocean','pbl_1d','mpas_atm'}
- case {'cam','tiegcm'}
-
varnames = get_DARTvars(fname);
num_vars = length(varnames);
- nlevels = dim_length(fname,'lev'); % determine # of state variables
vars = struct('model',model, ...
+ 'def_var',varnames{1}, ...
+ 'def_state_vars',[], ...
'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
- 'time_series_length',num_times, ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy) );
+ 'num_copies',num_copies, ...
+ 'num_ens_members',ens_size, ...
+ 'ensemble_indices',ens_indices, ...
+ 'time',dates, ...
+ 'time_series_length',num_times);
vars.vars = varnames;
vars.fname = fname;
-
- case 'pbl_1d'
-
- % A more robust way would be to use the netcdf low-level ops:
- % bob = var(f); % bob is a cell array of ncvars
- % name(bob{1}) % is the variable name string
- % bob{1}(:) % is the value of the netcdf variable (no offset/scale)
-
- num_vars = 22; % ps, t, u, v
- z_level = dim_length(fname, 'z_level'); % determine # of state variables
- sl_level = dim_length(fname,'sl_level'); % determine # of state variables
- times = nc_varget(fname,'time');
- z_level = nc_varget(fname,'z_level');
- sl_level = nc_varget(fname,'sl_level');
-
- vars = struct('model',model, ...
- 'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
- 'time_series_length',num_times, ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy));
-
- vars.fname = fname;
-
- case 'pe2lyr'
-
- % Since this is a 3D model, only the most rudimentary information
- % is gathered here. Each plot requires different information,
- % so there is a separate function (GetPe2lyrInfo.m) that gets
- % the information for each specific plot type.
-
- varnames = {'u','v','z'};
- num_vars = length(varnames);
-
- vars = struct('model',model, ...
- 'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
- 'time_series_length',num_times, ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy) );
-
- vars.vars = varnames;
- vars.fname = fname;
-
- case 'mitgcm_ocean'
-
- % A more robust way would be to use the netcdf low-level ops:
- % bob = var(f); % bob is a cell array of ncvars
- % name(bob{1}) % is the variable name string
- % bob{1}(:) % is the value of the netcdf variable (no offset/scale)
- % have not yet figured out a way to only use non-coordinate variables.
-
- varnames = {'S','T','U','V','SSH'};
- num_vars = length(varnames);
- nlevels = dim_length(fname,'ZG'); % determine # of state variables
-
- vars = struct('model',model, ...
- 'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
- 'time_series_length',num_times, ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy) );
-
- vars.vars = varnames;
- vars.fname = fname;
-
- case 'wrf'
-
- % requires a 'domain' and 'bottom_top_d01' dimension.
- % without both of these, it will fail in an ugly fashion.
-
- varnames = get_DARTvars(fname);
- num_vars = length(varnames);
- dinfo = nc_getdiminfo(fname,'domain');
- num_domains = dinfo.Length;
- dinfo = nc_getdiminfo(fname,'bottom_top_d01');
- num_levels = dinfo.Length;
-
- vars = struct('model',model, ...
- 'num_state_vars',num_vars, ...
- 'num_ens_members',num_copies, ...
- 'time_series_length',num_times, ...
- 'num_unstaggered_levels',num_levels, ...
- 'num_domains',num_domains, ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy));
-
- vars.vars = varnames;
- vars.fname = fname;
-
+
otherwise
error('model %s unknown',model)
@@ -309,14 +235,20 @@
end
-
-
function x = dim_length(fname,dimname)
+% Check for the existence of the named dimension and return it
+% if it exists. If it does not, error out with a useful message.
-y = nc_isvar(fname,dimname);
-if (y < 1)
- error('%s has no %s dimension/coordinate variable',fname,dimname)
+info = nc_info(fname);
+n = length(dimname);
+x = [];
+for i = 1:length(info.Dimension),
+ if ( strncmp(info.Dimension(i).Name, dimname, n) > 0 )
+ x = info.Dimension(i).Length;
+ break
+ end
end
-bob = nc_getdiminfo(fname,dimname);
-x = bob.Length;
+if isempty(x)
+ error('%s has no dimension named %s',fname,dimname)
+end
Modified: DART/branches/development/matlab/CheckModelCompatibility.m
===================================================================
--- DART/branches/development/matlab/CheckModelCompatibility.m 2012-04-11 03:29:46 UTC (rev 5687)
+++ DART/branches/development/matlab/CheckModelCompatibility.m 2012-04-11 03:38:37 UTC (rev 5688)
@@ -18,7 +18,6 @@
% $Revision$
% $Date$
-
if (nargin == 1) % better be a pinfo struct with at least these fields
file1 = arg1.truth_file; % string
file2 = arg1.diagn_file; % string
@@ -46,7 +45,7 @@
error('%s has no ''model'' global attribute.',file1)
end
-tnum_copies = dim_length(file1,'copy');
+tvars = get_DARTvars(file1);
tnum_times = dim_length(file1,'time');
times = nc_varget( file1,'time');
timeunits = nc_attget( file1,'time','units');
@@ -54,7 +53,7 @@
timeorigin = datenum(timebase(1),timebase(2),timebase(3));
ttimes = times + timeorigin;
-[tnum_vars,tdims] = ModelDimension(file1,tmodel);
+[tnum_vars,~] = ModelDimension(file1,tmodel);
if ( tnum_vars <= 0 )
error('Unable to determine resolution of %s.',file1)
end
@@ -66,7 +65,7 @@
error('%s has no ''model'' global attribute.',file2)
end
-dnum_copies = dim_length(file2,'copy');
+dvars = get_DARTvars(file2);
dnum_times = dim_length(file2,'time');
times = nc_varget( file2,'time');
timeunits = nc_attget( file2,'time','units');
@@ -74,12 +73,12 @@
timeorigin = datenum(timebase(1),timebase(2),timebase(3));
dtimes = times + timeorigin;
-[dnum_vars,ddims] = ModelDimension(file2,dmodel);
+[dnum_vars,~] = ModelDimension(file2,dmodel);
if ( dnum_vars <= 0 )
error('Unable to determine resolution of %s.',file2)
end
-% rudimentary bulletproofing
+%% rudimentary bulletproofing
if (strcmp(tmodel,dmodel) ~= 1)
fprintf('%s has model %s\n',file1,tmodel)
fprintf('%s has model %s\n',file2,dmodel)
@@ -93,13 +92,17 @@
error('no No NO ... both files must have same shape of state variables.')
end
-% if the lengths of the time arrays did not match, this used to be an
-% error. now we call a function to try to find any overlapping regions
-% in the time arrays and pass them back up to the called in the pinfo struct.
-% they then get used to extract the corresponding hyperslabs of data for
-% the matching times.
+% find the variables common to both files.
+vars = cell(0);
+for i = 1:length(dvars),
+ if any(strcmpi(dvars(i),tvars))
+ vars(length(vars)+1) = dvars(i);
+ end
+end
+pinfo_out.vars = vars;
+pinfo_out.num_state_vars = length(vars);
-% construct the pinfo struct in this function
+% Call a function to find the indices of theln - times common to both files.
pinfo_out = timearray_intersect(pinfo_out, file1, file2, ttimes, dtimes);
% fail here if the times had nothing in common.
@@ -114,15 +117,22 @@
error('These files have no timesteps in common')
end
+% Set the array of the times common to both files.
+T1 = pinfo_out.truth_time(1);
+T2 = pinfo_out.truth_time(1) + pinfo_out.truth_time(2) - 1;
+pinfo_out.time = ttimes(T1:T2);
+pinfo_out.time_series_length = pinfo_out.truth_time(2);
-%----------------
+
+
+function pret = timearray_intersect(pinfo, file1, file2, times1, times2)
% min1,max1 and min2,max2 are the index numbers of the intersection of the
% two input arrays. -1s in those numbers means no intersection. 1, length()
% means identical (could add a separate flag to simplify the calling code).
-function pret = timearray_intersect(pinfo, file1, file2, times1, times2)
-% for floating point comparisons, must be within this (single precision)
-% roundoff
+%% for floating point comparisons, must be within this (single precision)
+% roundoff
+
epsilon = 0.0000001;
% default return; no intersection
@@ -136,10 +146,10 @@
% a constant delta or not? compute delta array and validate those match?
% (to within an epsilon with floating pt roundoff)
-% check for the no-brainer case - identical time arrays.
-% watch out for the floating point compares, and the min/max are probably
-% redundant with the (1) and (l) comparisons, but until we put in checks
-% for monotonicity, it's a cheap safety check.
+%% check for the no-brainer case - identical time arrays.
+% watch out for the floating point compares, and the min/max are probably
+% redundant with the (1) and (l) comparisons, but until we put in checks
+% for monotonicity, it's a cheap safety check.
len = length(times1);
if ( (length(times1) == length(times2)) ...
&& (abs(min(times1) - min(times2)) < epsilon) ...
@@ -148,11 +158,12 @@
&& (times1(len) == times2(len)))
pret.truth_time = [1,len]; % start/count
pret.diagn_time = [1,len]; % start/count
+ pret.time = times1; % the common times in datenum-compatible form.
return
end
-% A is whichever array has the lower min. this reduces the number of
-% cases below we have to check for.
+%% A is whichever array has the lower min. this reduces the number of
+% cases below we have to check for.
if (min(times1) < min(times2))
A = times1;
B = times2;
@@ -161,9 +172,9 @@
A = times2;
end
-% precompute the data max, min, lengths using the A,B assignments
-% also, if differences are < epsilon, force equality to simplify
-% the comparison code below
+%% precompute the data max, min, lengths using the A,B assignments
+% also, if differences are < epsilon, force equality to simplify
+% the comparison code below
lenA = length(A);
lenB = length(B);
minA = min(A);
@@ -174,14 +185,14 @@
if (abs(maxA - minB) < epsilon) , minB = maxA; end
if (abs(maxA - maxB) < epsilon) , maxB = maxA; end
-% case 1: disjoint regions; simply return here because
-% return struct was initialized to the 'no intersection' case.
+%% case 1: disjoint regions; simply return here because
+% return struct was initialized to the 'no intersection' case.
if ((minA < minB) && (maxA < minB))
return
end
-% case 2: B fully contained in A; return corresponding index nums of overlap
-% include equal start & end points in this case.
+%% case 2: B fully contained in A; return corresponding index nums of overlap
+% include equal start & end points in this case.
if ((minA <= minB) && (maxB <= maxA))
minI = find(abs(A - minB) < epsilon);
maxI = find(abs(A - maxB) < epsilon);
@@ -195,8 +206,8 @@
maxJ = find(abs(B - maxA) < epsilon);
end
-% now map back to the original input order - this test must match exactly
-% the one used initially to assign A and B above.
+%% now map back to the original input order - this test must match exactly
+% the one used initially to assign A and B above.
if (min(times1) < min(times2))
min1 = minI;
max1 = maxI;
@@ -211,33 +222,40 @@
% now put the indices in the return struct and we are done.
pret.truth_time = [min1, max1-min1+1]; % start,count
-pret.diagn_time = [min2, min2-max2+1]; % start,count
+pret.diagn_time = [min2, max2-min2+1]; % start,count
+pret.time = times1(min1:max1); % the common times in datenum-compatible form.
% return here
function x = dim_length(fname,dimname)
+%% Check for the existence of the named dimension and return it
+% if it exists. If it does not, error out with a useful message.
-y = nc_isvar(fname,dimname);
-if (y < 1)
- error('%s has no %s dimension/coordinate variable',fname,dimname)
+info = nc_info(fname);
+n = length(dimname);
+x = [];
+for i = 1:length(info.Dimension),
+ if ( strncmp(info.Dimension(i).Name, dimname, n) > 0 )
+ x = info.Dimension(i).Length;
+ break
+ end
end
-bob = nc_getdiminfo(fname,dimname);
-x = bob.Length;
+if isempty(x)
+ error('%s has no dimension named %s',fname,dimname)
+end
function [x,y] = ModelDimension(fname,modelname)
-% Check the base geometry of the grid
-x = 0;
-y = NaN;
+%% Check the base geometry of the grid
switch lower(modelname)
case 'wrf'
- diminfo = nc_getdiminfo(fname, 'west_east_d01'); dnum_lons = diminfo.Length;
- diminfo = nc_getdiminfo(fname,'south_north_d01'); dnum_lats = diminfo.Length;
- diminfo = nc_getdiminfo(fname, 'bottom_top_d01'); dnum_lvls = diminfo.Length;
+ dnum_lons = dim_length(fname, 'west_east_d01');
+ dnum_lats = dim_length(fname,'south_north_d01');
+ dnum_lvls = dim_length(fname, 'bottom_top_d01');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
@@ -269,6 +287,13 @@
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
+ case 'mpas_atm'
+
+ dnum_cells = dim_length(fname,'nCells');
+ dnum_lvls = dim_length(fname,'nVertLevels');
+ x = 2;
+ y = [dnum_cells dnum_lvls];
+
case 'lorenz_96_2scale'
dnum_X = dim_length(fname,'Xdim');
dnum_Y = dim_length(fname,'Ydim');
Modified: DART/branches/development/matlab/GetBgridInfo.m
===================================================================
--- DART/branches/development/matlab/GetBgridInfo.m 2012-04-11 03:29:46 UTC (rev 5687)
+++ DART/branches/development/matlab/GetBgridInfo.m 2012-04-11 03:38:37 UTC (rev 5688)
@@ -1,4 +1,4 @@
-function pinfo = GetBgridInfo(pinfo_in,fname,routine);
+function pinfo = GetBgridInfo(pinfo_in,fname,routine)
%% GetBgridInfo prepares a structure of information needed by the subsequent "routine"
% The information is gathered via rudimentary "input" routines.
%
@@ -23,159 +23,141 @@
pinfo = pinfo_in;
model = nc_attget(fname, nc_global, 'model');
-if strcmp(lower(model),'fms_bgrid') ~= 1
+if strcmpi(model,'fms_bgrid') ~= 1
error('Not so fast, this is not a bgrid model.')
end
-copy = nc_varget(fname,'copy');
-times = nc_varget(fname,'time');
levels = nc_varget(fname,'lev');
TmpI = nc_varget(fname,'TmpI'); % temperature/pressure grid longitude
TmpJ = nc_varget(fname,'TmpJ'); % temperature/pressure grid latitude
VelI = nc_varget(fname,'VelI'); % velocity grid longitude
VelJ = nc_varget(fname,'VelJ'); % velocity grid latitude
-prognostic_vars = {'ps','t','u','v'};
-
-% Coordinate between time types and dates
-
-timeunits = nc_attget(fname,'time','units');
-timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
-timeorigin = datenum(timebase(1),timebase(2),timebase(3));
-dates = times + timeorigin;
-
switch lower(deblank(routine))
case {'plotbins','plotenserrspread','plotensmeantimeseries','plotenstimeseries'}
- pgvar = GetVar(prognostic_vars); % Determine prognostic variable
- [level, lvlind] = GetLevel(pgvar,levels); % Determine level and index
- [lat , latind] = GetLatitude( pgvar,TmpJ,VelJ);
- [lon , lonind] = GetLongitude(pgvar,TmpI,VelI);
+ pgvar = GetVar(pinfo.vars); % Determine prognostic variable
+ [level, lvlind] = GetLevel( pgvar, levels); % Determine level and index
+ [lat , latind] = GetLatitude( pgvar, TmpJ, VelJ);
+ [lon , lonind] = GetLongitude(pgvar, TmpI, VelI);
- pinfo = setfield(pinfo, 'model', model);
- pinfo = setfield(pinfo, 'fname', fname);
- pinfo = setfield(pinfo, 'times', dates);
- pinfo = setfield(pinfo, 'var', pgvar);
- pinfo = setfield(pinfo, 'level', level);
- pinfo = setfield(pinfo, 'levelindex', lvlind);
- pinfo = setfield(pinfo, 'longitude', lon);
- pinfo = setfield(pinfo, 'lonindex', lonind);
- pinfo = setfield(pinfo, 'latitude', lat);
- pinfo = setfield(pinfo, 'latindex',latind);
+ pinfo.fname = fname;
+ pinfo.var = pgvar;
+ pinfo.level = level;
+ pinfo.levelindex = lvlind;
+ pinfo.longitude = lon;
+ pinfo.lonindex = lonind;
+ pinfo.latitude = lat;
+ pinfo.latindex = latind;
case 'plotcorrel'
disp('Getting information for the ''base'' variable.')
- base_var = GetVar(prognostic_vars);
- [base_time, base_tmeind] = GetTime( base_var,dates);
- [base_lvl, base_lvlind] = GetLevel( base_var,levels);
- [base_lat, base_latind] = GetLatitude( base_var,TmpJ,VelJ);
- [base_lon, base_lonind] = GetLongitude(base_var,TmpI,VelI);
+ base_var = GetVar(pinfo.vars);
+ [base_time, base_tmeind] = GetTime(pinfo.time);
+ [base_lvl, base_lvlind] = GetLevel( base_var, levels);
+ [base_lat, base_latind] = GetLatitude( base_var, TmpJ, VelJ);
+ [base_lon, base_lonind] = GetLongitude(base_var, TmpI, VelI);
disp('Getting information for the ''comparison'' variable.')
- comp_var = GetVar(prognostic_vars, base_var);
- [comp_lvl, comp_lvlind] = GetLevel( comp_var,levels, base_lvl);
+ comp_var = GetVar(pinfo.vars, base_var);
+ [comp_lvl, comp_lvlind] = GetLevel(comp_var,levels, base_lvl);
- pinfo = setfield(pinfo, 'model', model);
- pinfo = setfield(pinfo, 'fname', fname);
- pinfo = setfield(pinfo, 'times', dates);
- pinfo = setfield(pinfo, 'base_var', base_var);
- pinfo = setfield(pinfo, 'comp_var', comp_var);
- pinfo = setfield(pinfo, 'base_time', base_time);
- pinfo = setfield(pinfo, 'base_tmeind', base_tmeind);
- pinfo = setfield(pinfo, 'base_lvl', base_lvl);
- pinfo = setfield(pinfo, 'base_lvlind', base_lvlind);
- pinfo = setfield(pinfo, 'base_lat', base_lat);
- pinfo = setfield(pinfo, 'base_latind', base_latind);
- pinfo = setfield(pinfo, 'base_lon', base_lon);
- pinfo = setfield(pinfo, 'base_lonind', base_lonind);
- pinfo = setfield(pinfo, 'comp_lvl', comp_lvl);
- pinfo = setfield(pinfo, 'comp_lvlind', comp_lvlind);
+ pinfo.fname = fname;
+ pinfo.base_var = base_var;
+ pinfo.comp_var = comp_var;
+ pinfo.base_time = base_time;
+ pinfo.base_tmeind = base_tmeind;
+ pinfo.base_lvl = base_lvl;
+ pinfo.base_lvlind = base_lvlind;
+ pinfo.base_lat = base_lat;
+ pinfo.base_latind = base_latind;
+ pinfo.base_lon = base_lon;
+ pinfo.base_lonind = base_lonind;
+ pinfo.comp_lvl = comp_lvl;
+ pinfo.comp_lvlind = comp_lvlind;
case 'plotvarvarcorrel'
disp('Getting information for the ''base'' variable.')
- base_var = GetVar(prognostic_vars);
- [base_time, base_tmeind] = GetTime( base_var,dates);
+ base_var = GetVar(pinfo.vars);
+ [base_time, base_tmeind] = GetTime(pinfo.time);
[base_lvl , base_lvlind] = GetLevel( base_var,levels);
[base_lat , base_latind] = GetLatitude( base_var,TmpJ,VelJ);
[base_lon , base_lonind] = GetLongitude(base_var,TmpI,VelI);
disp('Getting information for the ''comparison'' variable.')
- comp_var = GetVar(prognostic_vars, base_var);
- [comp_lvl, comp_lvlind] = GetLevel( comp_var,levels, base_lvl);
- [comp_lat, comp_latind] = GetLatitude( comp_var,TmpJ,VelJ, base_lat);
- [comp_lon, comp_lonind] = GetLongitude(comp_var,TmpI,VelI, base_lon);
+ comp_var = GetVar(pinfo.vars, base_var);
+ [comp_lvl, comp_lvlind] = GetLevel( comp_var, levels, base_lvl);
+ [comp_lat, comp_latind] = GetLatitude( comp_var, TmpJ, VelJ, base_lat);
+ [comp_lon, comp_lonind] = GetLongitude(comp_var, TmpI, VelI, base_lon);
- pinfo = setfield(pinfo, 'model', model);
- pinfo = setfield(pinfo, 'fname', fname);
- pinfo = setfield(pinfo, 'times', dates);
- pinfo = setfield(pinfo, 'base_var', base_var);
- pinfo = setfield(pinfo, 'comp_var', comp_var);
- pinfo = setfield(pinfo, 'base_time', base_time);
- pinfo = setfield(pinfo, 'base_tmeind', base_tmeind);
- pinfo = setfield(pinfo, 'base_lvl', base_lvl);
- pinfo = setfield(pinfo, 'base_lvlind', base_lvlind);
- pinfo = setfield(pinfo, 'base_lat', base_lat);
- pinfo = setfield(pinfo, 'base_latind', base_latind);
- pinfo = setfield(pinfo, 'base_lon', base_lon);
- pinfo = setfield(pinfo, 'base_lonind', base_lonind);
- pinfo = setfield(pinfo, 'comp_lvl', comp_lvl);
- pinfo = setfield(pinfo, 'comp_lvlind', comp_lvlind);
- pinfo = setfield(pinfo, 'comp_lat', comp_lat);
- pinfo = setfield(pinfo, 'comp_latind', comp_latind);
- pinfo = setfield(pinfo, 'comp_lon', comp_lon);
- pinfo = setfield(pinfo, 'comp_lonind', comp_lonind);
+ pinfo.fname = fname;
+ pinfo.base_var = base_var;
+ pinfo.comp_var = comp_var;
+ pinfo.base_time = base_time;
+ pinfo.base_tmeind = base_tmeind;
+ pinfo.base_lvl = base_lvl;
+ pinfo.base_lvlind = base_lvlind;
+ pinfo.base_lat = base_lat;
+ pinfo.base_latind = base_latind;
+ pinfo.base_lon = base_lon;
+ pinfo.base_lonind = base_lonind;
+ pinfo.comp_lvl = comp_lvl;
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list