[Dart-dev] [6143] DART/branches/development: The oned obs_diag program can now skip the spinup.
nancy at ucar.edu
nancy at ucar.edu
Thu May 16 16:24:26 MDT 2013
Revision: 6143
Author: thoar
Date: 2013-05-16 16:24:25 -0600 (Thu, 16 May 2013)
Log Message:
-----------
The oned obs_diag program can now skip the spinup.
I removed the useless 'rat_cri' and 'input_qc_threshold' because
the DART QC values are available to do the same thing.
The matlab scripts and input namelists had to be modified
to not have 'rat_cri' and 'input_qc_threshold' in them.
The input namelists also had some other variables removed,
like 'iskip_days' (replaced by init_skip_days, augmented by init_skip_seconds)
and 'obs_select' ... which was horrible to begin with. 'obs_select' required
you to specify the INTEGER depicting the obs TYPE you wanted to explore.
the two_experiments_evolution.m and plot_rank_histogram.m functions also support
the new 1D obs_diag_output.nc. The oned/obs_diag.html is updated but does not
have any example graphics. Mebbe later ...
Modified Paths:
--------------
DART/branches/development/diagnostics/matlab/plot_bias_xxx_profile.m
DART/branches/development/diagnostics/matlab/plot_evolution.m
DART/branches/development/diagnostics/matlab/plot_profile.m
DART/branches/development/diagnostics/matlab/plot_rank_histogram.m
DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m
DART/branches/development/diagnostics/matlab/plot_rmse_xxx_profile.m
DART/branches/development/diagnostics/matlab/two_experiments_evolution.m
DART/branches/development/diagnostics/oned/obs_diag.f90
DART/branches/development/diagnostics/oned/obs_diag.html
DART/branches/development/diagnostics/oned/obs_diag.nml
DART/branches/development/models/9var/work/input.nml
DART/branches/development/models/CESM/matlab/plot_evolution.m
DART/branches/development/models/CESM/matlab/plot_rmse_xxx_evolution.m
DART/branches/development/models/CESM/matlab/plot_spread_xxx_evolution.m
DART/branches/development/models/CESM/matlab/two_experiments_evolution.m
DART/branches/development/models/POP/matlab/plot_evolution.m
DART/branches/development/models/POP/matlab/plot_rmse_xxx_evolution.m
DART/branches/development/models/POP/matlab/plot_spread_xxx_evolution.m
DART/branches/development/models/POP/matlab/two_experiments_evolution.m
DART/branches/development/models/clm/matlab/plot_evolution.m
DART/branches/development/models/clm/matlab/plot_rmse_xxx_evolution.m
DART/branches/development/models/dynamo/work/input.nml
DART/branches/development/models/forced_lorenz_96/work/input.nml
DART/branches/development/models/ikeda/work/input.nml
DART/branches/development/models/lorenz_04/work/input.nml
DART/branches/development/models/lorenz_63/work/input.nml
DART/branches/development/models/lorenz_84/work/input.nml
DART/branches/development/models/lorenz_96/work/input.nml
DART/branches/development/models/lorenz_96/work/input.workshop.nml
DART/branches/development/models/lorenz_96_2scale/work/input.nml
DART/branches/development/models/null_model/work/input.nml
DART/branches/development/models/simple_advection/work/input.nml
DART/branches/development/models/template/work/input.nml
-------------- next part --------------
Modified: DART/branches/development/diagnostics/matlab/plot_bias_xxx_profile.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_bias_xxx_profile.m 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/matlab/plot_bias_xxx_profile.m 2013-05-16 22:24:25 UTC (rev 6143)
@@ -41,7 +41,6 @@
plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
-plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
Modified: DART/branches/development/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_evolution.m 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/matlab/plot_evolution.m 2013-05-16 22:24:25 UTC (rev 6143)
@@ -85,17 +85,15 @@
plotdat.region_names = deblank(plotdat.region_names');
end
-dimensionality = local_nc_attget(fname, nc_global, 'LocationRank');
-plotdat.binseparation = local_nc_attget(fname, nc_global, 'bin_separation');
-plotdat.binwidth = local_nc_attget(fname, nc_global, 'bin_width');
-time_to_skip = local_nc_attget(fname, nc_global, 'time_to_skip');
-plotdat.rat_cri = local_nc_attget(fname, nc_global, 'rat_cri');
-plotdat.input_qc_threshold = local_nc_attget(fname, nc_global, 'input_qc_threshold');
-plotdat.lonlim1 = local_nc_attget(fname, nc_global, 'lonlim1');
-plotdat.lonlim2 = local_nc_attget(fname, nc_global, 'lonlim2');
+dimensionality = nc_attget(fname, nc_global, 'LocationRank');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
+plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
+time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
plotdat.latlim1 = local_nc_attget(fname, nc_global, 'latlim1');
plotdat.latlim2 = local_nc_attget(fname, nc_global, 'latlim2');
-plotdat.biasconv = local_nc_attget(fname, nc_global, 'bias_convention');
% Coordinate between time types and dates
Modified: DART/branches/development/diagnostics/matlab/plot_profile.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_profile.m 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/matlab/plot_profile.m 2013-05-16 22:24:25 UTC (rev 6143)
@@ -41,7 +41,6 @@
plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
-plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
Modified: DART/branches/development/diagnostics/matlab/plot_rank_histogram.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_rank_histogram.m 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/matlab/plot_rank_histogram.m 2013-05-16 22:24:25 UTC (rev 6143)
@@ -84,28 +84,32 @@
plotdat.fname = fname;
plotdat.timecenters = nc_varget(fname,'time');
plotdat.timeedges = nc_varget(fname,'time_bounds');
-plotdat.mlevel = nc_varget(fname,'mlevel');
-plotdat.plevel = nc_varget(fname,'plevel');
-plotdat.plevel_edges = nc_varget(fname,'plevel_edges');
-plotdat.hlevel = nc_varget(fname,'hlevel');
-plotdat.hlevel_edges = nc_varget(fname,'hlevel_edges');
+plotdat.mlevel = local_nc_varget(fname,'mlevel');
+plotdat.plevel = local_nc_varget(fname,'plevel');
+plotdat.plevel_edges = local_nc_varget(fname,'plevel_edges');
+plotdat.hlevel = local_nc_varget(fname,'hlevel');
+plotdat.hlevel_edges = local_nc_varget(fname,'hlevel_edges');
plotdat.Nrhbins = length(nc_varget(fname,'rank_bins'));
plotdat.ncopies = length(nc_varget(fname,'copy'));
plotdat.nregions = length(nc_varget(fname,'region'));
plotdat.region_names = strtrim(nc_varget(fname,'region_names'));
+% Matlab does a crazy thing with the strings for 1 region
+if (plotdat.nregions == 1 && (size(plotdat.region_names,2) == 1) )
+ plotdat.region_names = deblank(plotdat.region_names');
+end
+
+dimensionality = nc_attget(fname, nc_global, 'LocationRank');
plotdat.timeseparation = nc_attget(fname, nc_global, 'bin_separation');
plotdat.timewidth = nc_attget(fname, nc_global, 'bin_width');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
-plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
-plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
+plotdat.outlierstring = nc_attget(fname, nc_global, 'outliers_in_histogram');
+plotdat.QCsusedstring = nc_attget(fname, nc_global, 'DART_QCs_in_histogram');
plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
-plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
-plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
-plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.outlierstring = nc_attget(fname, nc_global, 'outliers_in_histogram');
-plotdat.QCsusedstring = nc_attget(fname, nc_global, 'DART_QCs_in_histogram');
+plotdat.latlim1 = local_nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = local_nc_attget(fname, nc_global, 'latlim2');
% Make sure the time index makes sense.
@@ -114,12 +118,6 @@
length(plotdat.timecenters), timeindex)
end
-% Matlab does a crazy thing with the strings for 1 region
-
-if (plotdat.nregions == 1 && (size(plotdat.region_names,2) == 1) )
- plotdat.region_names = deblank(plotdat.region_names');
-end
-
% Coordinate between time types and dates
calendar = nc_attget(fname,'time','calendar');
@@ -198,10 +196,13 @@
guessdims = nc_var_dims(fname, plotdat.guessvar);
rhistdims = nc_var_dims(fname, plotdat.rhistvar);
- if ( findstr('surface',guessdims{3}) > 0 )
+ if ( dimensionality == 1 ) % observations on a unit circle, no level
+ plotdat.level = 1;
+ plotdat.level_units = [];
+ elseif ( strfind(guessdims{3},'surface') > 0 )
plotdat.level = 1;
plotdat.level_units = 'surface';
- elseif ( findstr('undef',guessdims{3}) > 0 )
+ elseif ( strfind(guessdims{3},'undef') > 0 )
plotdat.level = 1;
plotdat.level_units = 'undefined';
else
@@ -270,10 +271,14 @@
plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
- plotdat.title = sprintf('%s @ %d %s', ...
- plotdat.myvarname, ...
- plotdat.level(ilevel), ...
- plotdat.level_units);
+ if ( isempty(plotdat.level_units) )
+ plotdat.title = plotdat.myvarname;
+ else
+ plotdat.title = sprintf('%s @ %d %s', ...
+ plotdat.myvarname, ...
+ plotdat.level(ilevel), ...
+ plotdat.level_units);
+ end
plotdat.rank_hist = squeeze(rhist(plotdat.timeindex, :, ilevel,iregion));
@@ -377,7 +382,7 @@
j = 0;
for i = 1:length(x.allvarnames)
- indx = findstr('time',x.allvardims{i});
+ indx = strfind(x.allvardims{i},'time');
if (indx > 0)
j = j + 1;
@@ -386,7 +391,7 @@
end
end
-[b,i,j] = unique(basenames);
+[~,i,j] = unique(basenames);
y = cell(length(i),1);
ydims = cell(length(i),1);
for k = 1:length(i)
@@ -398,22 +403,22 @@
function s = ReturnBase(string1)
-inds = findstr('_guess',string1);
+inds = strfind(string1,'_guess');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_analy',string1);
+inds = strfind(string1,'_analy');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_VPguess',string1);
+inds = strfind(string1,'_VPguess');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_VPanaly',string1);
+inds = strfind(string1,'_VPanaly');
if (inds > 0 )
s = string1(1:inds-1);
end
@@ -430,3 +435,38 @@
else
s = rankvar;
end
+
+
+
+function value = local_nc_varget(fname,varname)
+%% If the variable exists in the file, return the contents of the variable.
+% if the variable does not exist, return empty value instead of error-ing
+% out.
+
+[variable_present, varid] = nc_var_exists(fname,varname);
+if (variable_present)
+ ncid = netcdf.open(fname);
+ value = netcdf.getVar(ncid, varid);
+else
+ value = [];
+end
+
+
+
+function value = local_nc_attget(fname,varid,varname)
+%% If the (global) attribute exists, return the value.
+% If it does not, do not throw a hissy-fit.
+
+value = [];
+if (varid == nc_global)
+ finfo = ncinfo(fname);
+ for iatt = 1:length(finfo.Attributes)
+ if (strcmp(finfo.Attributes(iatt).Name, deblank(varname)))
+ value = finfo.Attributes(iatt).Value;
+ return
+ end
+ end
+else
+ fprintf('function not supported for local variables, only global atts.\n')
+end
+
Modified: DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m 2013-05-16 22:24:25 UTC (rev 6143)
@@ -88,8 +88,6 @@
plotdat.binseparation = local_nc_attget(fname, nc_global, 'bin_separation');
plotdat.binwidth = local_nc_attget(fname, nc_global, 'bin_width');
time_to_skip = local_nc_attget(fname, nc_global, 'time_to_skip');
-plotdat.rat_cri = local_nc_attget(fname, nc_global, 'rat_cri');
-plotdat.input_qc_threshold = local_nc_attget(fname, nc_global, 'input_qc_threshold');
plotdat.lonlim1 = local_nc_attget(fname, nc_global, 'lonlim1');
plotdat.lonlim2 = local_nc_attget(fname, nc_global, 'lonlim2');
plotdat.latlim1 = local_nc_attget(fname, nc_global, 'latlim1');
Modified: DART/branches/development/diagnostics/matlab/plot_rmse_xxx_profile.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_rmse_xxx_profile.m 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/matlab/plot_rmse_xxx_profile.m 2013-05-16 22:24:25 UTC (rev 6143)
@@ -41,7 +41,6 @@
plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
-plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
Modified: DART/branches/development/diagnostics/matlab/two_experiments_evolution.m
===================================================================
--- DART/branches/development/diagnostics/matlab/two_experiments_evolution.m 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/matlab/two_experiments_evolution.m 2013-05-16 22:24:25 UTC (rev 6143)
@@ -153,8 +153,8 @@
commondata{i}.copyindex = get_copy_index(filenames{i},copystring);
commondata{i}.lonlim1 = nc_attget(filenames{i},nc_global,'lonlim1');
commondata{i}.lonlim2 = nc_attget(filenames{i},nc_global,'lonlim2');
- commondata{i}.latlim1 = nc_attget(filenames{i},nc_global,'latlim1');
- commondata{i}.latlim2 = nc_attget(filenames{i},nc_global,'latlim2');
+ commondata{i}.latlim1 = local_nc_attget(filenames{i},nc_global,'latlim1');
+ commondata{i}.latlim2 = local_nc_attget(filenames{i},nc_global,'latlim2');
commondata{i}.region_names = nc_varget(filenames{i},'region_names');
@@ -215,22 +215,21 @@
plotdat.levelindex = levelindex;
plotdat.bincenters = nc_varget(fname,'time');
plotdat.binedges = nc_varget(fname,'time_bounds');
-plotdat.mlevel = nc_varget(fname,'mlevel');
-plotdat.plevel = nc_varget(fname,'plevel');
-plotdat.plevel_edges = nc_varget(fname,'plevel_edges');
-plotdat.hlevel = nc_varget(fname,'hlevel');
-plotdat.hlevel_edges = nc_varget(fname,'hlevel_edges');
+plotdat.mlevel = local_nc_varget(fname,'mlevel');
+plotdat.plevel = local_nc_varget(fname,'plevel');
+plotdat.plevel_edges = local_nc_varget(fname,'plevel_edges');
+plotdat.hlevel = local_nc_varget(fname,'hlevel');
+plotdat.hlevel_edges = local_nc_varget(fname,'hlevel_edges');
plotdat.ncopies = length(nc_varget(fname,'copy'));
+dimensionality = nc_attget(fname, nc_global, 'LocationRank');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
-plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
-plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
-plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
-plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
-plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
+plotdat.latlim1 = local_nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = local_nc_attget(fname, nc_global, 'latlim2');
% Coordinate between time types and dates
@@ -271,7 +270,10 @@
guessdims = nc_var_dims(fname, plotdat.priorvar);
analydims = nc_var_dims(fname, plotdat.postevar);
-if ( findstr('surface',guessdims{3}) > 0 )
+if ( dimensionality == 1 ) % observations on a unit circle, no level
+ plotdat.level = 1;
+ plotdat.level_units = [];
+elseif ( findstr('surface',guessdims{3}) > 0 )
plotdat.level = 1;
plotdat.level_units = 'surface';
plotdat.level_edges = [];
@@ -478,10 +480,14 @@
set(get(ax1,'Ylabel'),'String',ylabel,'Interpreter','none','FontSize',figdata.fontsize)
-th = title({deblank(plotobj.region_names(plotobj.region,:)), ...
+if ( isempty(plotobj.level_units) )
+ th = title(plotobj.varname)
+else
+ th = title({deblank(plotobj.region_names(plotobj.region,:)), ...
sprintf('%s @ %d %s', plotobj.varname, ...
plotobj.level(plotobj.levelindex), ...
plotobj.level_units)});
+end
set(th,'Interpreter','none','FontSize',figdata.fontsize);
@@ -600,3 +606,39 @@
'fontsize',fontsize, 'orientation',orientation);
clf; orient(gcf, orientation)
+
+
+
+function value = local_nc_varget(fname,varname)
+%% If the variable exists in the file, return the contents of the variable.
+% if the variable does not exist, return empty value instead of error-ing
+% out.
+
+[variable_present, varid] = nc_var_exists(fname,varname);
+if (variable_present)
+ ncid = netcdf.open(fname);
+ value = netcdf.getVar(ncid, varid);
+else
+ value = [];
+end
+
+
+
+function value = local_nc_attget(fname,varid,varname)
+%% If the (global) attribute exists, return the value.
+% If it does not, do not throw a hissy-fit.
+
+value = [];
+if (varid == nc_global)
+ finfo = ncinfo(fname);
+ for iatt = 1:length(finfo.Attributes)
+ if (strcmp(finfo.Attributes(iatt).Name, deblank(varname)))
+ value = finfo.Attributes(iatt).Value;
+ return
+ end
+ end
+else
+ fprintf('function not supported for local variables, only global atts.\n')
+end
+
+
Modified: DART/branches/development/diagnostics/oned/obs_diag.f90
===================================================================
--- DART/branches/development/diagnostics/oned/obs_diag.f90 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/oned/obs_diag.f90 2013-05-16 22:24:25 UTC (rev 6143)
@@ -144,40 +144,40 @@
character(len = stringlength), dimension(MaxTrusted) :: trusted_obs = 'null'
-integer :: bin_width_days = 0 ! width of the assimilation bin - seconds
-integer :: bin_width_seconds = 0 ! width of the assimilation bin - days
integer :: max_num_bins = 9999 ! maximum number of temporal bins to consider
-real(r8):: rat_cri = 4.0
-real(r8):: input_qc_threshold = 4.0 ! input obs QC values >= this not used
+integer :: bin_width_days = -1 ! width of the assimilation bin - seconds
+integer :: bin_width_seconds = -1 ! width of the assimilation bin - days
+integer :: init_skip_days = 0
+integer :: init_skip_seconds = 0
logical :: verbose = .false.
-logical :: outliers_in_histogram = .false.
+logical :: outliers_in_histogram = .true.
logical :: create_rank_histogram = .true.
! index 1 == region 1 == [0.0, 1.0) i.e. Entire domain
! index 2 == region 2 == [0.0, 0.5)
! index 3 == region 3 == [0.5, 1.0)
-integer :: Nregions = 3
+integer :: Nregions = MaxRegions
real(r8), dimension(MaxRegions) :: lonlim1 = (/ 0.0_r8, 0.0_r8, 0.5_r8, -1.0_r8 /)
real(r8), dimension(MaxRegions) :: lonlim2 = (/ 1.0_r8, 0.5_r8, 1.0_r8, -1.0_r8 /)
character(len=6), dimension(MaxRegions) :: reg_names = &
(/ 'whole ','yin ','yang ','bogus '/)
-namelist /obs_diag_nml/ obs_sequence_name, obs_sequence_list, &
- bin_width_days, bin_width_seconds, &
- max_num_bins, rat_cri, input_qc_threshold, &
- Nregions, lonlim1, lonlim2, reg_names, &
- verbose, outliers_in_histogram, &
+namelist /obs_diag_nml/ obs_sequence_name, obs_sequence_list, &
+ bin_width_days, bin_width_seconds, &
+ init_skip_days, init_skip_seconds, max_num_bins, &
+ Nregions, lonlim1, lonlim2, reg_names, &
+ verbose, outliers_in_histogram, &
create_rank_histogram, trusted_obs
!-----------------------------------------------------------------------
! Variables used to accumulate the statistics.
!-----------------------------------------------------------------------
-integer, parameter :: Ncopies = 20
+integer, parameter :: Ncopies = 18
character(len = stringlength), dimension(Ncopies) :: copy_names = &
- (/ 'Nposs ', 'Nused ', 'NbigQC ', 'NbadIZ ', &
+ (/ 'Nposs ', 'Nused ', &
'rmse ', 'bias ', 'spread ', 'totalspread', &
'NbadDARTQC ', 'observation', 'ens_mean ', &
'N_DARTqc_0 ', 'N_DARTqc_1 ', 'N_DARTqc_2 ', 'N_DARTqc_3 ', &
@@ -191,8 +191,6 @@
character(len=8) :: string
integer :: num_times = 0, num_regions = 0, num_variables = 0
integer, dimension(:,:,:), pointer :: Nposs, Nused, Ntrusted
- integer, dimension(:,:,:), pointer :: NbigQC ! # original QC values >= input_qc_threshold
- integer, dimension(:,:,:), pointer :: NbadIZ ! # bad (ie huge) Innovation Zscore
real(r8), dimension(:,:,:), pointer :: rmse, bias, spread, totspread
integer, dimension(:,:,:), pointer :: NbadDartQC ! # bad DART QC values
real(r8), dimension(:,:,:), pointer :: observation, ens_mean
@@ -216,7 +214,7 @@
integer :: iregion, iepoch, ivar, ifile, num_obs_in_epoch
real(r8) :: rlocation
-integer :: obsindex, i, iunit, ierr, io
+integer :: obsindex, i, iunit, ierr, io, ireg
integer :: seconds, days, Nepochs, Nfiles
integer :: num_trusted
@@ -237,7 +235,7 @@
type(time_type) :: binwidth, halfbinwidth
type(time_type) :: seqT1, seqTN ! first,last time of obs sequence
type(time_type) :: obsT1, obsTN ! first,last time of all observations
-type(time_type) :: obs_time
+type(time_type) :: obs_time, skip_time
character(len = 129) :: msgstring1, msgstring2
character(len = stringlength) :: str1, str2, str3
@@ -308,6 +306,8 @@
! presume the first/last times in the input file(s) are to be used.
!----------------------------------------------------------------------
+skip_time = set_time(init_skip_seconds, init_skip_days)
+
call DetermineFilenames(obsT1, obsTN, Nepochs, Nfiles) ! fills obs_seq_filenames array
call DefineTimes() ! Sets binwidth, halfbinwidth
@@ -324,8 +324,19 @@
obs_used_in_epoch = 0
!----------------------------------------------------------------------
+! Rectify the region namelist information
!----------------------------------------------------------------------
+ireg = MaxRegions
+Regions: do i = 1,MaxRegions
+ if ((lonlim1(i) < 0.0_r8) .or. (lonlim2(i) < 0.0_r8) ) then
+ exit Regions
+ else
+ ireg = i
+ endif
+enddo Regions
+Nregions = min(Nregions, ireg)
+
if (verbose) then
write(logfileunit,*)
write( * ,*)
@@ -337,10 +348,6 @@
enddo
endif
-prior_mean(1) = 0.0_r8
-prior_spread(1) = 0.0_r8
-posterior_mean(1) = 0.0_r8
-posterior_spread(1) = 0.0_r8
!----------------------------------------------------------------------
! Declares and initializes the guess and analy structures.
@@ -539,6 +546,9 @@
obs_err_var = get_obs_def_error_variance(obs_def)
rlocation = get_location(obs_loc)
+ ! Check to make sure we are past the burn-in
+ if (obs_time < skip_time) cycle ObservationLoop
+
! Check to see if this is a trusted observation
trusted = .false.
if ( num_trusted > 0 ) then
@@ -647,17 +657,6 @@
if ( .not. keeper ) cycle Areas
!-----------------------------------------------------------
- ! Count original QC values 'of interest' ...
- !-----------------------------------------------------------
-
- if ( org_qc_index > 0 ) then
- if (qc(org_qc_index) > input_qc_threshold ) then
- call IPE(guess%NbigQC(iepoch, iregion, flavor), 1)
- call IPE(analy%NbigQC(iepoch, iregion, flavor), 1)
- endif
- endif
-
- !-----------------------------------------------------------
! Count DART QC values
!-----------------------------------------------------------
@@ -670,25 +669,6 @@
call Bin3D(qc_integer, iepoch, iregion, flavor, trusted, obs(1), &
obs_err_var, pr_mean, pr_sprd, po_mean, po_sprd, rank_histogram_bin)
- !-----------------------------------------------------------
- ! Count 'large' innovations
- !-----------------------------------------------------------
-
- if( pr_zscore > rat_cri ) then
- call IPE(guess%NbadIZ(iepoch, iregion, flavor), 1)
- endif
-
- if( po_zscore > rat_cri ) then
- call IPE(analy%NbadIZ(iepoch, iregion, flavor), 1)
- endif
-
- if ( pr_zscore > rat_cri ) then
- if (verbose) then
- write(logfileunit,*)'obsindex ',obsindex,' pr_zscore ', pr_zscore
- write(logfileunit,*)'val,prm,prs,var ',obs(1), pr_mean, pr_sprd, obs_err_var
- endif
- endif
-
enddo Areas
!-----------------------------------------------------------------
@@ -743,8 +723,6 @@
write(*,*) '# observations used : ',sum(obs_used_in_epoch)
write(*,*) 'Count summary over all regions - obs may count for multiple regions:'
write(*,*) '# identity : ',Nidentity
-write(*,*) '# big Innov (ratio) : ',sum(guess%NbadIZ)
-write(*,*) '# big (original) QC : ',sum(guess%NbigQC)
write(*,*) '# bad DART QC prior : ',sum(guess%NbadDartQC)
write(*,*) '# bad DART QC post : ',sum(analy%NbadDartQC)
write(*,*) '# TRUSTED : ',sum(analy%Ntrusted)
@@ -766,14 +744,11 @@
write(*,*) '# poste DART QC 5 : ',sum(analy%NDartQC_5)
write(*,*) '# poste DART QC 6 : ',sum(analy%NDartQC_6)
write(*,*) '# poste DART QC 7 : ',sum(analy%NDartQC_7)
-write(*,*)
write(logfileunit,*)
write(logfileunit,*) '# observations used : ',sum(obs_used_in_epoch)
write(logfileunit,*) 'Count summary over all regions - obs may count for multiple regions:'
write(logfileunit,*) '# identity : ',Nidentity
-write(logfileunit,*) '# big Innov (ratio) : ',sum(guess%NbadIZ)
-write(logfileunit,*) '# big (original) QC : ',sum(guess%NbigQC)
write(logfileunit,*) '# bad DART QC prior : ',sum(guess%NbadDartQC)
write(logfileunit,*) '# bad DART QC post : ',sum(analy%NbadDartQC)
write(logfileunit,*) '# TRUSTED : ',sum(analy%Ntrusted)
@@ -795,7 +770,22 @@
write(logfileunit,*) '# poste DART QC 5 : ',sum(analy%NDartQC_5)
write(logfileunit,*) '# poste DART QC 6 : ',sum(analy%NDartQC_6)
write(logfileunit,*) '# poste DART QC 7 : ',sum(analy%NDartQC_7)
+
+! Print the histogram of innovations as a function of standard deviation.
+write( * ,*)
+write( * ,*) 'Table that reflects the outlier_threshold -- '
+write( * ,*) 'How are the (good) innovations distributed?'
write(logfileunit,*)
+write(logfileunit,*) 'Table that reflects the outlier_threshold -- '
+write(logfileunit,*) 'How are the (good) innovations distributed?'
+do i=0,MaxSigmaBins
+ if(nsigma(i) /= 0) then
+ write( * ,*)'innovations within ',i+1,' stdev = ',nsigma(i)
+ write(logfileunit,*)'innovations within ',i+1,' stdev = ',nsigma(i)
+ endif
+enddo
+write( * ,*)
+write(logfileunit,*)
!----------------------------------------------------------------------
! Open netCDF output file
@@ -830,8 +820,6 @@
guess%ens_mean( Nepochs, Nregions, num_obs_types), &
guess%Nposs( Nepochs, Nregions, num_obs_types), &
guess%Nused( Nepochs, Nregions, num_obs_types), &
- guess%NbigQC( Nepochs, Nregions, num_obs_types), &
- guess%NbadIZ( Nepochs, Nregions, num_obs_types), &
guess%NbadDartQC( Nepochs, Nregions, num_obs_types), &
guess%NDartQC_0( Nepochs, Nregions, num_obs_types), &
guess%NDartQC_1( Nepochs, Nregions, num_obs_types), &
@@ -851,8 +839,6 @@
guess%ens_mean = 0.0_r8
guess%Nposs = 0
guess%Nused = 0
-guess%NbigQC = 0
-guess%NbadIZ = 0
guess%NbadDartQC = 0
guess%NDartQC_0 = 0
guess%NDartQC_1 = 0
@@ -877,8 +863,6 @@
analy%ens_mean( Nepochs, Nregions, num_obs_types), &
analy%Nposs( Nepochs, Nregions, num_obs_types), &
analy%Nused( Nepochs, Nregions, num_obs_types), &
- analy%NbigQC( Nepochs, Nregions, num_obs_types), &
- analy%NbadIZ( Nepochs, Nregions, num_obs_types), &
analy%NbadDartQC( Nepochs, Nregions, num_obs_types), &
analy%NDartQC_0( Nepochs, Nregions, num_obs_types), &
analy%NDartQC_1( Nepochs, Nregions, num_obs_types), &
@@ -898,8 +882,6 @@
analy%ens_mean = 0.0_r8
analy%Nposs = 0
analy%Nused = 0
-analy%NbigQC = 0
-analy%NbadIZ = 0
analy%NbadDartQC = 0
analy%NDartQC_0 = 0
analy%NDartQC_1 = 0
@@ -928,7 +910,7 @@
deallocate(guess%rmse, guess%bias, guess%spread, guess%totspread, &
guess%observation, guess%ens_mean, guess%Nposs, guess%Nused, &
- guess%NbigQC, guess%NbadIZ, guess%NbadDartQC,guess%Ntrusted, &
+ guess%NbadDartQC,guess%Ntrusted, &
guess%NDartQC_0, guess%NDartQC_1, guess%NDartQC_2, guess%NDartQC_3, &
guess%NDartQC_4, guess%NDartQC_5, guess%NDartQC_6, guess%NDartQC_7)
@@ -938,7 +920,7 @@
deallocate(analy%rmse, analy%bias, analy%spread, analy%totspread, &
analy%observation, analy%ens_mean, analy%Nposs, analy%Nused, &
- analy%NbigQC, analy%NbadIZ, analy%NbadDartQC,analy%Ntrusted, &
+ analy%NbadDartQC,analy%Ntrusted, &
analy%NDartQC_0, analy%NDartQC_1, analy%NDartQC_2, analy%NDartQC_3, &
analy%NDartQC_4, analy%NDartQC_5, analy%NDartQC_6, analy%NDartQC_7)
@@ -1152,28 +1134,31 @@
! type(time_type), intent(out) GLOBAL :: halfbinwidth ! half that period
! bin_width_days, bin_width_seconds GLOBAL from namelist
-logical :: error_out = .false.
integer :: nbins
type(time_type) :: test_time
! do some error-checking first
-if ( (bin_width_days < 0) .or. (bin_width_seconds < 0) ) then
+if ( (bin_width_days < 0) .and. (bin_width_seconds >= 0) ) then
+
write(msgstring1,*)'bin_width_[days,seconds] must be non-negative, they are ', &
bin_width_days, bin_width_seconds
- call error_handler(E_WARN,'DefineTimes',msgstring1,source,revision,revdate)
- error_out = .true.
-endif
+ call error_handler(E_ERR,'DefineTimes',msgstring1,source,revision,revdate, &
+ text2='namelist parameter out-of-bounds. Fix and try again.')
-if ( error_out ) call error_handler(E_ERR,'DefineTimes', &
- 'namelist parameter out-of-bounds. Fix and try again.',source,revision,revdate)
+elseif ( (bin_width_days >= 0) .and. (bin_width_seconds < 0) ) then
-! 'space-filling' strategy: bin width and bin separation are same.
-! If the user input is empty, use default Nepochs.
-! If not, ...
+ write(msgstring1,*)'bin_width_[days,seconds] must be non-negative, they are ', &
+ bin_width_days, bin_width_seconds
+ call error_handler(E_ERR,'DefineTimes',msgstring1,source,revision,revdate, &
+ text2='namelist parameter out-of-bounds. Fix and try again.')
-if ((bin_width_days == 0) .and. (bin_width_seconds == 0)) then
+elseif ( (bin_width_days <= 0) .and. (bin_width_seconds <= 0) ) then
+
! This is the 'default' case ... use all possible, up to "max_num_bins".
+ ! 'space-filling' strategy: bin width and bin separation are same.
+ ! Using Nepochs that comes from the number of unique times in the files.
+
binwidth = (obsTN - obsT1) / (Nepochs - 1)
if (Nepochs > max_num_bins) then
write(msgstring1,*)'default calculation results in ',Nepochs,' time bins.'
@@ -1182,6 +1167,7 @@
Nepochs = max_num_bins
obsTN = obsT1 + (Nepochs-1)*binwidth
endif
+
else
! honor the user input
binwidth = set_time(bin_width_seconds, bin_width_days)
@@ -1802,14 +1788,10 @@
call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'bin_width_seconds', seconds ), &
'WriteNetCDF', 'put_att bin_width_seconds '//trim(fname))
call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'time_to_skip', &
- (/ 0, 0, 0, 0, 0, 0/) ), &
+ (/ 0, 0, init_skip_days, 0, 0, init_skip_seconds /) ), &
'WriteNetCDF', 'put_att time_to_skip '//trim(fname))
call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'max_num_bins', Nepochs ), &
'WriteNetCDF', 'put_att max_num_bins '//trim(fname))
- call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'rat_cri', rat_cri ), &
- 'WriteNetCDF', 'put_att rat_cri '//trim(fname))
- call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'input_qc_threshold', input_qc_threshold ), &
- 'WriteNetCDF', 'put_att input_qc_threshold '//trim(fname))
call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'Nregions', Nregions ), &
'WriteNetCDF', 'put_att Nregions '//trim(fname))
call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'lonlim1', lonlim1(1:Nregions) ), &
@@ -2129,24 +2111,22 @@
rchunk(iregion, 1,itime) = vrbl%Nposs( itime,iregion,ivar)
rchunk(iregion, 2,itime) = vrbl%Nused( itime,iregion,ivar)
- rchunk(iregion, 3,itime) = vrbl%NbigQC( itime,iregion,ivar)
- rchunk(iregion, 4,itime) = vrbl%NbadIZ( itime,iregion,ivar)
- rchunk(iregion, 5,itime) = vrbl%rmse( itime,iregion,ivar)
- rchunk(iregion, 6,itime) = vrbl%bias( itime,iregion,ivar)
- rchunk(iregion, 7,itime) = vrbl%spread( itime,iregion,ivar)
- rchunk(iregion, 8,itime) = vrbl%totspread( itime,iregion,ivar)
- rchunk(iregion, 9,itime) = vrbl%NbadDartQC( itime,iregion,ivar)
- rchunk(iregion,10,itime) = vrbl%observation(itime,iregion,ivar)
- rchunk(iregion,11,itime) = vrbl%ens_mean( itime,iregion,ivar)
- rchunk(iregion,12,itime) = vrbl%NDartQC_0( itime,iregion,ivar)
- rchunk(iregion,13,itime) = vrbl%NDartQC_1( itime,iregion,ivar)
- rchunk(iregion,14,itime) = vrbl%NDartQC_2( itime,iregion,ivar)
- rchunk(iregion,15,itime) = vrbl%NDartQC_3( itime,iregion,ivar)
- rchunk(iregion,16,itime) = vrbl%NDartQC_4( itime,iregion,ivar)
- rchunk(iregion,17,itime) = vrbl%NDartQC_5( itime,iregion,ivar)
- rchunk(iregion,18,itime) = vrbl%NDartQC_6( itime,iregion,ivar)
- rchunk(iregion,19,itime) = vrbl%NDartQC_7( itime,iregion,ivar)
- rchunk(iregion,20,itime) = vrbl%Ntrusted( itime,iregion,ivar)
+ rchunk(iregion, 3,itime) = vrbl%rmse( itime,iregion,ivar)
+ rchunk(iregion, 4,itime) = vrbl%bias( itime,iregion,ivar)
+ rchunk(iregion, 5,itime) = vrbl%spread( itime,iregion,ivar)
+ rchunk(iregion, 6,itime) = vrbl%totspread( itime,iregion,ivar)
+ rchunk(iregion, 7,itime) = vrbl%NbadDartQC( itime,iregion,ivar)
+ rchunk(iregion, 8,itime) = vrbl%observation(itime,iregion,ivar)
+ rchunk(iregion, 9,itime) = vrbl%ens_mean( itime,iregion,ivar)
+ rchunk(iregion,10,itime) = vrbl%NDartQC_0( itime,iregion,ivar)
+ rchunk(iregion,11,itime) = vrbl%NDartQC_1( itime,iregion,ivar)
+ rchunk(iregion,12,itime) = vrbl%NDartQC_2( itime,iregion,ivar)
+ rchunk(iregion,13,itime) = vrbl%NDartQC_3( itime,iregion,ivar)
+ rchunk(iregion,14,itime) = vrbl%NDartQC_4( itime,iregion,ivar)
+ rchunk(iregion,15,itime) = vrbl%NDartQC_5( itime,iregion,ivar)
+ rchunk(iregion,16,itime) = vrbl%NDartQC_6( itime,iregion,ivar)
+ rchunk(iregion,17,itime) = vrbl%NDartQC_7( itime,iregion,ivar)
+ rchunk(iregion,18,itime) = vrbl%Ntrusted( itime,iregion,ivar)
enddo
enddo
Modified: DART/branches/development/diagnostics/oned/obs_diag.html
===================================================================
--- DART/branches/development/diagnostics/oned/obs_diag.html 2013-05-16 21:41:58 UTC (rev 6142)
+++ DART/branches/development/diagnostics/oned/obs_diag.html 2013-05-16 22:24:25 UTC (rev 6143)
@@ -2,44 +2,85 @@
"http://www.w3.org/TR/html4/strict.dtd">
<HTML>
<HEAD>
-<TITLE>program obs_diag</TITLE>
+<TITLE>program obs_diag (for 1D observations)</TITLE>
<link rel="stylesheet" type="text/css" href="../../doc/html/doc.css">
</HEAD>
<BODY>
<A NAME="TOP"></A>
-<center>
-<A HREF="#Modules">MODULES</A> /
+<H1>PROGRAM <em class=program>obs_diag</em> (for 1D observations)</H1>
+
+<table border=0 summary="" cellpadding=5>
+<tr>
+ <td valign=middle>
+ <img src="../../doc/html/Dartboard7.png" alt="DART project logo" height=70 />
+ </td>
+ <td>
+ <P>Jump to <a href="../../index.html">DART Documentation Main Index</a><br />
+ <small><small>version information for this file: <br />
+ <!-- version tag follows, do not edit -->
+ $Id$</small></small>
+ </P></td>
+</tr>
+</table>
+
<A HREF="#Namelist">NAMELIST</A> /
+<A HREF="#Modules">INTERFACES</A> /
<A HREF="#FilesUsed">FILES</A> /
<A HREF="#References">REFERENCES</A> /
<A HREF="#Errors">ERRORS</A> /
<A HREF="#FuturePlans">PLANS</A> /
<A HREF="#PrivateComponents">PRIVATE COMPONENTS</A> /
<A HREF="#Legalese">TERMS OF USE</A>
-</center>
-<H1>PROGRAM <em class=program>obs_diag</em></H1>
-<!-- version tag follows, do not edit --><P>$Id$</P>
+<H2>Overview/Usage</H2>
<P>
Main program for observation-space diagnostics for the models with
- 1D locations. The model space is carved up into regions, hence the
- dependency on the location type. There is no discrimination between
- observation types, but the framework is still in the code - the main
- selection block is commented out. The result of the code is a set of
- ASCII files that contain the RMS of the bias and the RMS of the spread
- of the prior (aka 'guess') and posterior (aka 'analysis') estimates
- as a function of time and region. Separate (pairs of)
- files are created for each observation type in the observation
- sequence file.
-<br><br>
- The <em class=program>obs_diag</em> program performs an outlier test on
- observations, which can be more restrictive than the outlier test in the
- filter, but not less
- (see <em class=code>outlier_threshold</em> in <em class=code>filter_nml</em>
- and <em class=code>rat_cri</em> here).
-<br><br>
+ 1D locations. 18 quantities are calculated for each region for each
+ temporal bin specified by user input. The result of the code is a netCDF
+ file that contains the 18 quantities of the prior (aka 'guess') and
+ posterior (aka 'analysis') estimates as a function of time and region as
+ well as all the metadata to create meaningful figures.
+</P>
+<P>
+ Each <em class="file">obs_seq.final</em> file contains an observation
+ sequence that has multiple 'copies' of the observation. One copy is
+ the actual observation, another copy is the prior ensemble mean estimate
+ of the observation, one is the spread of the prior ensemble estimate,
+ one may be the prior estimate from ensemble member 1, ... etc.
+ The only observations for the 1D models are the result of a
+ 'perfect model' experiment, so there is an additional copy called
+ the 'truth' - the noise-free expected observation given the true
+ model state. Since this copy does not, in general, exist for the
+ high-order models, all comparisons are made with the copy labelled
+ 'observation'. It may, in some instances be useful to compare against
+ the 'truth', in which case you will have to hand-edit the code and
+ recompile. Caveat emptor.
+</P>
+<P>
+ Each ensemble member applies a forward operator to the state to compute
+ the "expected" value of an observation. Given multiple estimates of
+ the observation, several quantities can be calculated. It is possible to
+ compute the expected observations from the state vector before
+ assimilating (the "guess", "forecast", or "prior") or after the
+ assimilation (the "analysis", or "posterior").
+</P>
+<P>
+ Even with <em class="file">input.nml</em>:<em class="code">filter_nml:num_output_obs_members</em>
+ set to <em class="code">0</em>; the full [prior,posterior] ensemble mean
+ and [prior,posterior] ensemble spread are preserved in the
+ <em class="file">obs_seq.final</em> file. Consequently, the ensemble
+ means and spreads are used to calculate the diagnostics. If the
+ <em class="file">input.nml</em>:<em class="code">filter_nml:num_output_obs_members</em>
+ is set to <em class="code">80</em> (for example); the first 80 ensemble
+ members prior and posterior "expected" values of the observation are
+ also included. In this case, the <em class="file">obs_seq.final</em>
+ file contains enough information to calculate a rank histograms,
+ verify forecasts, etc.
+ The ensemble means are still used for many other calculations.
+</P>
+<P>
The J release of DART also implements a DART QC flag that provides
information about how the observation was or was not assimilated.
The DART QC flag is intended to provide information about whether the
@@ -70,37 +111,90 @@
</table>
<P>
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list