[Dart-dev] [4167] DART/trunk/diagnostics: The obs_diag and plotting routines now support

nancy at ucar.edu nancy at ucar.edu
Wed Nov 25 14:36:59 MST 2009


Revision: 4167
Author:   thoar
Date:     2009-11-25 14:36:59 -0700 (Wed, 25 Nov 2009)
Log Message:
-----------
The obs_diag and plotting routines now support 

Modified Paths:
--------------
    DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
    DART/trunk/diagnostics/matlab/plot_profile.m
    DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
    DART/trunk/diagnostics/threed_sphere/obs_diag.f90

-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m	2009-11-25 21:25:19 UTC (rev 4166)
+++ DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m	2009-11-25 21:36:59 UTC (rev 4167)
@@ -61,7 +61,7 @@
 
 diminfo                    = nc_getdiminfo(fname,'region');
 plotdat.nregions           = diminfo.Length;
-region_names               = nc_varget(fname,'region_names') 
+region_names               = nc_varget(fname,'region_names');
 plotdat.region_names       = deblank(region_names);
 
 % Coordinate between time types and dates
@@ -114,15 +114,16 @@
 
    % get appropriate vertical coordinate variable
 
-   guessdims = nc_var_dims(fname, plotdat.guessvar);
-   analydims = nc_var_dims(fname, plotdat.analyvar);
+   guessdims = nc_var_dims(  fname, plotdat.guessvar);
+   analydims = nc_var_dims(  fname, plotdat.analyvar);
+   varinfo   = nc_getvarinfo(fname, plotdat.analyvar);
 
    if ( findstr('surface',guessdims{2}) > 0 )
-      disp(sprintf('%s is a surface field.',plotdat.guessvar))
-      error('Cannot display a surface field this way.')
+      fprintf('%s is a surface field.\n',plotdat.guessvar)
+      fprintf('Cannot display a surface field this way.\n')
    elseif ( findstr('undef',guessdims{2}) > 0 )
-      disp(sprintf('%s has no vertical definition.',plotdat.guessvar))
-      error('Cannot display this field this way.')
+      fprintf('%s has no vertical definition.\n',plotdat.guessvar)
+      fprintf('Cannot display this field this way.\n')
    end
 
    [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
@@ -149,26 +150,36 @@
    
    guess = nc_varget(fname, plotdat.guessvar);  
    analy = nc_varget(fname, plotdat.analyvar); 
+   n = size(analy);
+  
+   % singleton dimensions are auto-squeezed - which is unfortunate.
+   % We want these things to be 3D. [copy-level-region]
+   % Sometimes there is one region, sometimes one level, ...
+   % To complicate matters, the stupid 'ones' function does not allow
+   % the last dimension to be unity ... so you have double the size
+   % of the array ...
 
-   % Check for one region ... if the last dimension is a singleton 
-   % dimension, it is auto-squeezed  - which is bad.
-   % We want these things to be 3D.
-
-   n = size(guess);
-   if ( length(n) < 3 )
-       bob = NaN*ones(n(1),n(2),2);
-       ted = NaN*ones(n(1),n(2),2);
+   if ( plotdat.nregions == 1 )
+       bob = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
+       ted = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
        bob(:,:,1) = guess;
        ted(:,:,1) = analy;
        guess = bob; clear bob
        analy = ted; clear ted
+   elseif ( plotdat.nlevels == 1 )
+       bob = NaN*ones(varinfo.Size);
+       ted = NaN*ones(varinfo.Size);
+       bob(:,1,:) = guess;
+       ted(:,1,:) = analy;
+       guess = bob; clear bob
+       analy = ted; clear ted
    end
    
    % check to see if there is anything to plot
    nposs = sum(guess(plotdat.Npossindex,:,:));
 
    if ( sum(nposs(:)) < 1 )
-      disp(sprintf('No obs for %s...  skipping', plotdat.varnames{ivar}))
+      fprintf('No obs for %s...  skipping\n', plotdat.varnames{ivar})
       continue
    end
 
@@ -276,6 +287,7 @@
 
    axlims = [plotdat.Xrange plotdat.Yrange];
    axis(axlims)
+
    set(gca,'YDir', plotdat.YDir)
    hold on; plot([0 0],plotdat.Yrange,'k-')
 
@@ -312,7 +324,7 @@
            'XTick',          xticks,'XTicklabel',num2str(nicexticks))
        
    set(get(ax2,'Xlabel'),'String','# of obs (o=poss, +=used)')
-   set(get(ax1,'Xlabel'),'String',plotdat.xlabel)
+   set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
    set(ax1,'Position',get(ax2,'Position'))
    grid
 
@@ -380,7 +392,7 @@
 
 
 
-function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname);
+function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname)
 % Find the vertical dimension and harvest some info
 
 varinfo  = nc_getvarinfo(fname,varname);
@@ -450,13 +462,16 @@
       ymax =  ceil(max(glommed));
    end
 
+   if (ymin == ymax)
+     ymin = ymin - 0.1*ymin;
+     ymax = ymax + 0.1*ymax;
+   end
+
    Yrange = [ymin ymax];
 
-   x = [min([Yrange(1) 0.0]) Yrange(2)];
+   x = sort([min([Yrange(1) 0.0]) Yrange(2)] ,'ascend');
 end
 
-
-
 function Stripes(x,edges)
 % EraseMode: [ {normal} | background | xor | none ]
 
@@ -465,7 +480,6 @@
 hold on;
 for i = 1:2:(length(edges)-1)
   yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ];
-  h = fill(xc,yc,[0.8 0.8 0.8], ...
-  'EraseMode','background','EdgeColor','none');
+  fill(xc,yc,[0.8 0.8 0.8],'EraseMode','background','EdgeColor','none');
 end
 hold off;

Modified: DART/trunk/diagnostics/matlab/plot_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_profile.m	2009-11-25 21:25:19 UTC (rev 4166)
+++ DART/trunk/diagnostics/matlab/plot_profile.m	2009-11-25 21:36:59 UTC (rev 4167)
@@ -119,10 +119,10 @@
 
    if ( findstr('surface',guessdims{2}) > 0 )
       fprintf('%s is a surface field.\n',plotdat.guessvar)
-      error('Cannot display a surface field this way.')
+      fprintf('Cannot display a surface field this way.\n')
    elseif ( findstr('undef',guessdims{2}) > 0 )
       fprintf('%s has no vertical definition.\n',plotdat.guessvar)
-      error('Cannot display this field this way.')
+      fprintf('Cannot display this field this way.\n')
    end
 
    [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
@@ -314,7 +314,7 @@
            'XTick',          xticks,'XTicklabel',num2str(nicexticks))
        
    set(get(ax2,'Xlabel'),'String','# of obs (o=poss, +=used)')
-   set(get(ax1,'Xlabel'),'String',plotdat.xlabel)
+   set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
    set(ax1,'Position',get(ax2,'Position'))
    grid
 
@@ -372,7 +372,7 @@
    end
 end
 
-[~,i,~] = unique(basenames);
+[b,i,j] = unique(basenames);
 
 for j = 1:length(i)
    disp(sprintf('%2d is %s',j,basenames{j}))
@@ -452,9 +452,14 @@
       ymax =  ceil(max(glommed));
    end
 
+   if (ymin == ymax)
+     ymin = ymin - 0.1*ymin;
+     ymax = ymax + 0.1*ymax;
+   end
+
    Yrange = [ymin ymax];
 
-   x = [min([Yrange(1) 0.0]) Yrange(2)];
+   x = sort([min([Yrange(1) 0.0]) Yrange(2)] ,'ascend');
 end
 
 function Stripes(x,edges)

Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m	2009-11-25 21:25:19 UTC (rev 4166)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m	2009-11-25 21:36:59 UTC (rev 4167)
@@ -61,7 +61,7 @@
 
 diminfo                    = nc_getdiminfo(fname,'region');
 plotdat.nregions           = diminfo.Length;
-region_names               = nc_varget(fname,'region_names') 
+region_names               = nc_varget(fname,'region_names');
 plotdat.region_names       = deblank(region_names);
 
 % Coordinate between time types and dates
@@ -114,15 +114,16 @@
 
    % get appropriate vertical coordinate variable
 
-   guessdims = nc_var_dims(fname, plotdat.guessvar);
-   analydims = nc_var_dims(fname, plotdat.analyvar);
+   guessdims = nc_var_dims(  fname, plotdat.guessvar);
+   analydims = nc_var_dims(  fname, plotdat.analyvar);
+   varinfo   = nc_getvarinfo(fname, plotdat.analyvar);
 
    if ( findstr('surface',guessdims{2}) > 0 )
-      disp(sprintf('%s is a surface field.',plotdat.guessvar))
-      error('Cannot display a surface field this way.')
+      fprintf('%s is a surface field.\n',plotdat.guessvar)
+      fprintf('Cannot display a surface field this way.\n')
    elseif ( findstr('undef',guessdims{2}) > 0 )
-      disp(sprintf('%s has no vertical definition.',plotdat.guessvar))
-      error('Cannot display this field this way.')
+      fprintf('%s has no vertical definition.\n',plotdat.guessvar)
+      fprintf('Cannot display this field this way.\n')
    end
 
    [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
@@ -149,26 +150,36 @@
    
    guess = nc_varget(fname, plotdat.guessvar);  
    analy = nc_varget(fname, plotdat.analyvar); 
+   n = size(analy);
+  
+   % singleton dimensions are auto-squeezed - which is unfortunate.
+   % We want these things to be 3D. [copy-level-region]
+   % Sometimes there is one region, sometimes one level, ...
+   % To complicate matters, the stupid 'ones' function does not allow
+   % the last dimension to be unity ... so you have double the size
+   % of the array ...
 
-   % Check for one region ... if the last dimension is a singleton 
-   % dimension, it is auto-squeezed  - which is bad.
-   % We want these things to be 3D.
-
-   n = size(guess);
-   if ( length(n) < 3 )
-       bob = NaN*ones(n(1),n(2),2);
-       ted = NaN*ones(n(1),n(2),2);
+   if ( plotdat.nregions == 1 )
+       bob = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
+       ted = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
        bob(:,:,1) = guess;
        ted(:,:,1) = analy;
        guess = bob; clear bob
        analy = ted; clear ted
+   elseif ( plotdat.nlevels == 1 )
+       bob = NaN*ones(varinfo.Size);
+       ted = NaN*ones(varinfo.Size);
+       bob(:,1,:) = guess;
+       ted(:,1,:) = analy;
+       guess = bob; clear bob
+       analy = ted; clear ted
    end
    
    % check to see if there is anything to plot
    nposs = sum(guess(plotdat.Npossindex,:,:));
 
    if ( sum(nposs(:)) < 1 )
-      disp(sprintf('No obs for %s...  skipping', plotdat.varnames{ivar}))
+      fprintf('No obs for %s...  skipping\n', plotdat.varnames{ivar})
       continue
    end
 
@@ -276,6 +287,7 @@
 
    axlims = [plotdat.Xrange plotdat.Yrange];
    axis(axlims)
+
    set(gca,'YDir', plotdat.YDir)
    hold on; plot([0 0],plotdat.Yrange,'k-')
 
@@ -312,7 +324,7 @@
            'XTick',          xticks,'XTicklabel',num2str(nicexticks))
        
    set(get(ax2,'Xlabel'),'String','# of obs (o=poss, +=used)')
-   set(get(ax1,'Xlabel'),'String',plotdat.xlabel)
+   set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
    set(ax1,'Position',get(ax2,'Position'))
    grid
 
@@ -380,7 +392,7 @@
 
 
 
-function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname);
+function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname)
 % Find the vertical dimension and harvest some info
 
 varinfo  = nc_getvarinfo(fname,varname);
@@ -450,9 +462,14 @@
       ymax =  ceil(max(glommed));
    end
 
+   if (ymin == ymax)
+     ymin = ymin - 0.1*ymin;
+     ymax = ymax + 0.1*ymax;
+   end
+
    Yrange = [ymin ymax];
 
-   x = [min([Yrange(1) 0.0]) Yrange(2)];
+   x = sort([min([Yrange(1) 0.0]) Yrange(2)] ,'ascend');
 end
 
 function Stripes(x,edges)
@@ -463,7 +480,6 @@
 hold on;
 for i = 1:2:(length(edges)-1)
   yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ];
-  h = fill(xc,yc,[0.8 0.8 0.8], ...
-  'EraseMode','background','EdgeColor','none');
+  fill(xc,yc,[0.8 0.8 0.8],'EraseMode','background','EdgeColor','none');
 end
 hold off;

Modified: DART/trunk/diagnostics/threed_sphere/obs_diag.f90
===================================================================
--- DART/trunk/diagnostics/threed_sphere/obs_diag.f90	2009-11-25 21:25:19 UTC (rev 4166)
+++ DART/trunk/diagnostics/threed_sphere/obs_diag.f90	2009-11-25 21:36:59 UTC (rev 4167)
@@ -1,4 +1,4 @@
-! Data Assimilation Research Testbed -- DART
+!d Data Assimilation Research Testbed -- DART
 ! Copyright 2004-2007, Data Assimilation Research Section
 ! University Corporation for Atmospheric Research
 ! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
@@ -197,10 +197,13 @@
 ! Variables used to accumulate the statistics.
 !-----------------------------------------------------------------------
 
-integer, parameter :: Ncopies = 10
-character(len = stringlength), dimension(Ncopies) :: copy_names = &
-   (/ 'Nposs      ', 'Nused      ', 'NbadQC     ', 'NbadIZ     ', 'NbadUV     ', &
-      'NbadLV     ', 'rmse       ', 'bias       ', 'spread     ', 'totalspread' /)
+integer, parameter :: Ncopies = 21
+character(len = stringlength), dimension(Ncopies) :: copy_names =                &
+   (/ 'Nposs      ', 'Nused      ', 'NbigQC     ', 'NbadIZ     ', 'NbadUV     ', &
+      'NbadLV     ', 'rmse       ', 'bias       ', 'spread     ', 'totalspread', &
+      'NbadDARTQC ', 'observation', 'ens_mean   ',                               &
+      'N_DARTqc_0 ', 'N_DARTqc_1 ', 'N_DARTqc_2 ', 'N_DARTqc_3 ',                &
+      'N_DARTqc_4 ', 'N_DARTqc_5 ', 'N_DARTqc_6 ', 'N_DARTqc_7 ' /)
 
 type TLRV_type
    ! statistics by time-level-region-variable
@@ -211,11 +214,15 @@
    character(len=8) :: string
    integer :: num_times = 0, num_levels = 0, num_regions = 0, num_variables = 0
    integer,  dimension(:,:,:,:), pointer :: Nposs, Nused
-   integer,  dimension(:,:,:,:), pointer :: NbadQC ! # bad (original) QC values
+   integer,  dimension(:,:,:,:), pointer :: NbigQC ! # original QC values >= input_qc_threshold
    integer,  dimension(:,:,:,:), pointer :: NbadIZ ! # bad (ie huge) Innovation Zscore
    integer,  dimension(:,:,:,:), pointer :: NbadUV ! # unmatched U/V wind pairs
    integer,  dimension(:,:,:,:), pointer :: NbadLV ! # obs above/below top/bottom
    real(r8), dimension(:,:,:,:), pointer :: rmse, bias, spread, totspread
+   integer,  dimension(:,:,:,:), pointer :: NbadDartQC ! # bad DART QC values
+   real(r8), dimension(:,:,:,:), pointer :: observation, ens_mean
+   integer,  dimension(:,:,:,:), pointer :: NDartQC_0, NDartQC_1, NDartQC_2, NDartQC_3
+   integer,  dimension(:,:,:,:), pointer :: NDartQC_4, NDartQC_5, NDartQC_6, NDartQC_7
 end type TLRV_type
 
 type LRV_type
@@ -226,11 +233,15 @@
    character(len=8) :: string
    integer :: num_levels = 0, num_regions = 0, num_variables = 0
    integer,  dimension(:,:,:), pointer :: Nposs, Nused
-   integer,  dimension(:,:,:), pointer :: NbadQC ! # bad (original) QC values
+   integer,  dimension(:,:,:), pointer :: NbigQC ! # bad (original) QC values
    integer,  dimension(:,:,:), pointer :: NbadIZ ! # bad (ie huge) Innovation Zscore
    integer,  dimension(:,:,:), pointer :: NbadUV ! # unmatched U/V wind pairs
    integer,  dimension(:,:,:), pointer :: NbadLV ! # obs above/below top/bottom
    real(r8), dimension(:,:,:), pointer :: rmse, bias, spread, totspread
+   integer,  dimension(:,:,:), pointer :: NbadDartQC ! # bad DART QC values
+   real(r8), dimension(:,:,:), pointer :: observation, ens_mean
+   integer,  dimension(:,:,:), pointer :: NDartQC_0, NDartQC_1, NDartQC_2, NDartQC_3
+   integer,  dimension(:,:,:), pointer :: NDartQC_4, NDartQC_5, NDartQC_6, NDartQC_7
 end type LRV_type
 
 type(TLRV_type) :: analy,    guess
@@ -392,37 +403,49 @@
 allocate(obs_seq_filenames(Nepochs*4))
 obs_seq_filenames = 'null'
 
-allocate(guess%rmse(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%bias(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%spread(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-      guess%totspread(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%Nposs( Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%Nused( Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%NbadQC(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%NbadIZ(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%NbadUV(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         guess%NbadLV(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%rmse(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%bias(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%spread(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-      analy%totspread(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%Nposs( Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%Nused( Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%NbadQC(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%NbadIZ(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%NbadUV(Nepochs, Nlevels, Nregions, num_obs_kinds), &
-         analy%NbadLV(Nepochs, Nlevels, Nregions, num_obs_kinds))
+allocate(guess%rmse(       Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%bias(       Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%spread(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%totspread(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%observation(Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%ens_mean(   Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%Nposs(      Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%Nused(      Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NbigQC(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NbadIZ(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NbadUV(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NbadLV(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_0(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_1(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_2(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_3(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_4(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_5(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_6(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         guess%NDartQC_7(  Nepochs, Nlevels, Nregions, num_obs_kinds)  )
 
-guess%rmse   = 0.0_r8
-guess%bias   = 0.0_r8
-guess%spread = 0.0_r8
-guess%totspread = 0.0_r8
-guess%Nposs  = 0
-guess%Nused  = 0
-guess%NbadQC = 0
-guess%NbadIZ = 0
-guess%NbadUV = 0
-guess%NbadLV = 0
+guess%rmse        = 0.0_r8
+guess%bias        = 0.0_r8
+guess%spread      = 0.0_r8
+guess%totspread   = 0.0_r8
+guess%observation = 0.0_r8
+guess%ens_mean    = 0.0_r8
+guess%Nposs       = 0
+guess%Nused       = 0
+guess%NbigQC      = 0
+guess%NbadIZ      = 0
+guess%NbadUV      = 0
+guess%NbadLV      = 0
+guess%NbadDartQC  = 0
+guess%NDartQC_0   = 0
+guess%NDartQC_1   = 0
+guess%NDartQC_2   = 0
+guess%NDartQC_3   = 0
+guess%NDartQC_4   = 0
+guess%NDartQC_5   = 0
+guess%NDartQC_6   = 0
+guess%NDartQC_7   = 0
 
 guess%string        = 'guess'
 guess%num_times     = Nepochs
@@ -430,71 +453,149 @@
 guess%num_regions   = Nregions
 guess%num_variables = num_obs_kinds
 
-analy%rmse   = 0.0_r8
-analy%bias   = 0.0_r8
-analy%spread = 0.0_r8
-analy%totspread = 0.0_r8
-analy%Nposs  = 0
-analy%Nused  = 0
-analy%NbadQC = 0
-analy%NbadIZ = 0
-analy%NbadUV = 0
-analy%NbadLV = 0
+allocate(analy%rmse(       Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%bias(       Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%spread(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%totspread(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%observation(Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%ens_mean(   Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%Nposs(      Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%Nused(      Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NbigQC(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NbadIZ(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NbadUV(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NbadLV(     Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_0(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_1(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_2(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_3(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_4(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_5(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_6(  Nepochs, Nlevels, Nregions, num_obs_kinds), &
+         analy%NDartQC_7(  Nepochs, Nlevels, Nregions, num_obs_kinds)  )
 
+analy%rmse        = 0.0_r8
+analy%bias        = 0.0_r8
+analy%spread      = 0.0_r8
+analy%totspread   = 0.0_r8
+analy%observation = 0.0_r8
+analy%ens_mean    = 0.0_r8
+analy%Nposs       = 0
+analy%Nused       = 0
+analy%NbigQC      = 0
+analy%NbadIZ      = 0
+analy%NbadUV      = 0
+analy%NbadLV      = 0
+analy%NbadDartQC  = 0
+analy%NDartQC_0   = 0
+analy%NDartQC_1   = 0
+analy%NDartQC_2   = 0
+analy%NDartQC_3   = 0
+analy%NDartQC_4   = 0
+analy%NDartQC_5   = 0
+analy%NDartQC_6   = 0
+analy%NDartQC_7   = 0
+
 analy%string        = 'analy'
 analy%num_times     = Nepochs
 analy%num_levels    = Nlevels
 analy%num_regions   = Nregions
 analy%num_variables = num_obs_kinds
 
-allocate(guessAVG%rmse(  Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%bias(  Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%spread(Nlevels, Nregions, num_obs_kinds), &
-      guessAVG%totspread(Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%Nposs( Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%Nused( Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%NbadQC(Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%NbadIZ(Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%NbadUV(Nlevels, Nregions, num_obs_kinds), &
-         guessAVG%NbadLV(Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%rmse(  Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%bias(  Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%spread(Nlevels, Nregions, num_obs_kinds), &
-      analyAVG%totspread(Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%Nposs( Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%Nused( Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%NbadQC(Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%NbadIZ(Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%NbadUV(Nlevels, Nregions, num_obs_kinds), &
-         analyAVG%NbadLV(Nlevels, Nregions, num_obs_kinds))
+allocate(guessAVG%rmse(       Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%bias(       Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%spread(     Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%totspread(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%observation(Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%ens_mean(   Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%Nposs(      Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%Nused(      Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NbigQC(     Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NbadIZ(     Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NbadUV(     Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NbadLV(     Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NbadDartQC( Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_0(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_1(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_2(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_3(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_4(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_5(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_6(  Nlevels, Nregions, num_obs_kinds), &
+         guessAVG%NDartQC_7(  Nlevels, Nregions, num_obs_kinds)  )
 
-guessAVG%rmse   = 0.0_r8
-guessAVG%bias   = 0.0_r8
-guessAVG%spread = 0.0_r8
-guessAVG%totspread = 0.0_r8
-guessAVG%Nposs  = 0
-guessAVG%Nused  = 0
-guessAVG%NbadQC = 0
-guessAVG%NbadIZ = 0
-guessAVG%NbadUV = 0
-guessAVG%NbadLV = 0
+guessAVG%rmse        = 0.0_r8
+guessAVG%bias        = 0.0_r8
+guessAVG%spread      = 0.0_r8
+guessAVG%totspread   = 0.0_r8
+guessAVG%observation = 0.0_r8
+guessAVG%ens_mean    = 0.0_r8
+guessAVG%Nposs       = 0
+guessAVG%Nused       = 0
+guessAVG%NbigQC      = 0
+guessAVG%NbadIZ      = 0
+guessAVG%NbadUV      = 0
+guessAVG%NbadLV      = 0
+guessAVG%NbadDartQC  = 0
+guessAVG%NDartQC_0   = 0
+guessAVG%NDartQC_1   = 0
+guessAVG%NDartQC_2   = 0
+guessAVG%NDartQC_3   = 0
+guessAVG%NDartQC_4   = 0
+guessAVG%NDartQC_5   = 0
+guessAVG%NDartQC_6   = 0
+guessAVG%NDartQC_7   = 0
 
 guessAVG%string        = 'VPguess'
 guessAVG%num_levels    = Nlevels
 guessAVG%num_regions   = Nregions
 guessAVG%num_variables = num_obs_kinds
 
-analyAVG%rmse   = 0.0_r8
-analyAVG%bias   = 0.0_r8
-analyAVG%spread = 0.0_r8
-analyAVG%totspread = 0.0_r8
-analyAVG%Nposs  = 0
-analyAVG%Nused  = 0
-analyAVG%NbadQC = 0
-analyAVG%NbadIZ = 0
-analyAVG%NbadUV = 0
-analyAVG%NbadLV = 0
+allocate(analyAVG%rmse(       Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%bias(       Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%spread(     Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%totspread(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%observation(Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%ens_mean(   Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%Nposs(      Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%Nused(      Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NbigQC(     Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NbadIZ(     Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NbadUV(     Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NbadLV(     Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NbadDartQC( Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_0(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_1(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_2(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_3(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_4(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_5(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_6(  Nlevels, Nregions, num_obs_kinds), &
+         analyAVG%NDartQC_7(  Nlevels, Nregions, num_obs_kinds)  )
 
+analyAVG%rmse        = 0.0_r8
+analyAVG%bias        = 0.0_r8
+analyAVG%spread      = 0.0_r8
+analyAVG%totspread   = 0.0_r8
+analyAVG%observation = 0.0_r8
+analyAVG%ens_mean    = 0.0_r8
+analyAVG%Nposs       = 0
+analyAVG%Nused       = 0
+analyAVG%NbigQC      = 0
+analyAVG%NbadIZ      = 0
+analyAVG%NbadUV      = 0
+analyAVG%NbadLV      = 0
+analyAVG%NbadDartQC  = 0
+analyAVG%NDartQC_0   = 0
+analyAVG%NDartQC_1   = 0
+analyAVG%NDartQC_2   = 0
+analyAVG%NDartQC_3   = 0
+analyAVG%NDartQC_4   = 0
+analyAVG%NDartQC_5   = 0
+analyAVG%NDartQC_6   = 0
+analyAVG%NDartQC_7   = 0
+
 analyAVG%string        = 'VPanaly'
 analyAVG%num_levels    = Nlevels
 analyAVG%num_regions   = Nregions
@@ -976,12 +1077,51 @@
 
             if( qc_index > 0 ) then
                if (qc(qc_index) > input_qc_threshold ) then
-               call IPE(guess%NbadQC(iepoch,level_index,iregion,flavor), 1)
-               call IPE(analy%NbadQC(iepoch,level_index,iregion,flavor), 1)
+               call IPE(guess%NbigQC(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NbigQC(iepoch,level_index,iregion,flavor), 1)
                endif
             endif
 
             !-----------------------------------------------------------
+            ! Count DART QC values 
+            ! FIXME ... should these be different for prior/posterior?
+            !-----------------------------------------------------------
+
+            if (        qc_integer == 0 ) then
+               call IPE(guess%NDartQC_0(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_0(iepoch,level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer == 1 ) then
+               call IPE(guess%NDartQC_1(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_1(iepoch,level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer == 2 ) then
+               call IPE(guess%NDartQC_2(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_2(iepoch,level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer == 3 ) then
+               call IPE(guess%NDartQC_3(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_3(iepoch,level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer == 4 ) then
+               call IPE(guess%NDartQC_4(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_4(iepoch,level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer == 5 ) then
+               call IPE(guess%NDartQC_5(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_5(iepoch,level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer == 6 ) then
+               call IPE(guess%NDartQC_6(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_6(iepoch,level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer == 7 ) then
+               call IPE(guess%NDartQC_7(iepoch,level_index,iregion,flavor), 1)
+               call IPE(analy%NDartQC_7(iepoch,level_index,iregion,flavor), 1)
+
+            endif
+
+            !-----------------------------------------------------------
             ! Count 'large' innovations
             !-----------------------------------------------------------
 
@@ -1065,11 +1205,50 @@
 
             if( qc_index > 0 ) then
                if (qc(qc_index) > input_qc_threshold ) then
-               call IPE(guessAVG%NbadQC(level_index,iregion,flavor), 1)
-               call IPE(analyAVG%NbadQC(level_index,iregion,flavor), 1)
+               call IPE(guessAVG%NbigQC(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NbigQC(level_index,iregion,flavor), 1)
                endif
             endif
 
+            ! Count DART QC values 
+            ! FIXME ... should these be different for prior/posterior?
+
+            if (        qc_integer    == 0 ) then
+               call IPE(guessAVG%NDartQC_0(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_0(level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer    == 1 ) then
+               call IPE(guessAVG%NDartQC_1(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_1(level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer    == 2 ) then
+               call IPE(guessAVG%NDartQC_2(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_2(level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer    == 3 ) then
+               call IPE(guessAVG%NDartQC_3(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_3(level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer    == 4 ) then
+               call IPE(guessAVG%NDartQC_4(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_4(level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer    == 5 ) then
+               call IPE(guessAVG%NDartQC_5(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_5(level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer    == 6 ) then
+               call IPE(guessAVG%NDartQC_6(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_6(level_index,iregion,flavor), 1)
+
+            elseif (    qc_integer    == 7 ) then
+               call IPE(guessAVG%NDartQC_7(level_index,iregion,flavor), 1)
+               call IPE(analyAVG%NDartQC_7(level_index,iregion,flavor), 1)
+
+            endif
+
+            ! Count 'large' innovations
+
             if(pr_zscore > rat_cri )  then
                call IPE(guessAVG%NbadIZ(level_index,iregion,flavor), 1)
             endif
@@ -1161,52 +1340,73 @@
 do iregion= 1,Nregions
 do ilev   = 1,Nlevels
 do iepoch = 1,Nepochs
-   if (  guess%Nused( iepoch, ilev, iregion, ivar) == 0) then
-         guess%rmse(  iepoch, ilev, iregion, ivar) = MISSING_R4
-         guess%bias(  iepoch, ilev, iregion, ivar) = MISSING_R4
-         guess%spread(iepoch, ilev, iregion, ivar) = MISSING_R4
-      guess%totspread(iepoch, ilev, iregion, ivar) = MISSING_R4
+
+   if (  guess%Nused(      iepoch, ilev, iregion, ivar) == 0) then
+         guess%observation(iepoch, ilev, iregion, ivar) = MISSING_R4
+         guess%ens_mean(   iepoch, ilev, iregion, ivar) = MISSING_R4
+         guess%bias(       iepoch, ilev, iregion, ivar) = MISSING_R4
+         guess%rmse(       iepoch, ilev, iregion, ivar) = MISSING_R4
+         guess%spread(     iepoch, ilev, iregion, ivar) = MISSING_R4
+         guess%totspread(  iepoch, ilev, iregion, ivar) = MISSING_R4
    else
-         guess%rmse(  iepoch, ilev, iregion, ivar) = &
-    sqrt(guess%rmse(  iepoch, ilev, iregion, ivar) / &
-         guess%Nused( iepoch, ilev, iregion, ivar) )
+         guess%observation(iepoch, ilev, iregion, ivar) = &
+         guess%observation(iepoch, ilev, iregion, ivar) / &
+         guess%Nused(      iepoch, ilev, iregion, ivar)  
 
-         guess%bias(  iepoch, ilev, iregion, ivar) = &
-         guess%bias(  iepoch, ilev, iregion, ivar) / &
-         guess%Nused( iepoch, ilev, iregion, ivar)  
+         guess%ens_mean(   iepoch, ilev, iregion, ivar) = &
+         guess%ens_mean(   iepoch, ilev, iregion, ivar) / &
+         guess%Nused(      iepoch, ilev, iregion, ivar)  
 
-         guess%spread(iepoch, ilev, iregion, ivar) = &
-    sqrt(guess%spread(iepoch, ilev, iregion, ivar) / &
-         guess%Nused( iepoch, ilev, iregion, ivar) )
+         guess%bias(       iepoch, ilev, iregion, ivar) = &
+         guess%bias(       iepoch, ilev, iregion, ivar) / &
+         guess%Nused(      iepoch, ilev, iregion, ivar)  
 
-         guess%totspread(iepoch, ilev, iregion, ivar) = &
-    sqrt(guess%totspread(iepoch, ilev, iregion, ivar) / &
-         guess%Nused(    iepoch, ilev, iregion, ivar) )
+         guess%rmse(       iepoch, ilev, iregion, ivar) = &
+    sqrt(guess%rmse(       iepoch, ilev, iregion, ivar) / &
+         guess%Nused(      iepoch, ilev, iregion, ivar) )
 
+         guess%spread(     iepoch, ilev, iregion, ivar) = &
+    sqrt(guess%spread(     iepoch, ilev, iregion, ivar) / &
+         guess%Nused(      iepoch, ilev, iregion, ivar) )
+
+         guess%totspread(  iepoch, ilev, iregion, ivar) = &
+    sqrt(guess%totspread(  iepoch, ilev, iregion, ivar) / &
+         guess%Nused(      iepoch, ilev, iregion, ivar) )
+
    endif
 
-   if (  analy%Nused( iepoch, ilev, iregion, ivar) == 0) then
-         analy%rmse(  iepoch, ilev, iregion, ivar) = MISSING_R4
-         analy%bias(  iepoch, ilev, iregion, ivar) = MISSING_R4
-         analy%spread(iepoch, ilev, iregion, ivar) = MISSING_R4
-      analy%totspread(iepoch, ilev, iregion, ivar) = MISSING_R4
+   if (  analy%Nused(      iepoch, ilev, iregion, ivar) == 0) then
+         analy%observation(iepoch, ilev, iregion, ivar) = MISSING_R4
+         analy%ens_mean(   iepoch, ilev, iregion, ivar) = MISSING_R4
+         analy%bias(       iepoch, ilev, iregion, ivar) = MISSING_R4
+         analy%rmse(       iepoch, ilev, iregion, ivar) = MISSING_R4
+         analy%spread(     iepoch, ilev, iregion, ivar) = MISSING_R4
+         analy%totspread(  iepoch, ilev, iregion, ivar) = MISSING_R4
    else
-         analy%rmse(  iepoch, ilev, iregion, ivar) = &
-    sqrt(analy%rmse(  iepoch, ilev, iregion, ivar) / &
-         analy%Nused( iepoch, ilev, iregion, ivar) )
+         analy%observation(iepoch, ilev, iregion, ivar) = &
+         analy%observation(iepoch, ilev, iregion, ivar) / &
+         analy%Nused(      iepoch, ilev, iregion, ivar)  
 
-         analy%bias(  iepoch, ilev, iregion, ivar) = &
-         analy%bias(  iepoch, ilev, iregion, ivar) / &
-         analy%Nused( iepoch, ilev, iregion, ivar)  
+         analy%ens_mean( iepoch, ilev, iregion, ivar) = &
+         analy%ens_mean( iepoch, ilev, iregion, ivar) / &
+         analy%Nused(      iepoch, ilev, iregion, ivar)  
 
-         analy%spread(iepoch, ilev, iregion, ivar) = &
-    sqrt(analy%spread(iepoch, ilev, iregion, ivar) / &
-         analy%Nused( iepoch, ilev, iregion, ivar) )
+         analy%bias(       iepoch, ilev, iregion, ivar) = &
+         analy%bias(       iepoch, ilev, iregion, ivar) / &
+         analy%Nused(      iepoch, ilev, iregion, ivar)  
 
-         analy%totspread(iepoch, ilev, iregion, ivar) = &
-    sqrt(analy%totspread(iepoch, ilev, iregion, ivar) / &
-         analy%Nused(    iepoch, ilev, iregion, ivar) )
+         analy%rmse(       iepoch, ilev, iregion, ivar) = &
+    sqrt(analy%rmse(       iepoch, ilev, iregion, ivar) / &
+         analy%Nused(      iepoch, ilev, iregion, ivar) )
 
+         analy%spread(     iepoch, ilev, iregion, ivar) = &
+    sqrt(analy%spread(     iepoch, ilev, iregion, ivar) / &
+         analy%Nused(      iepoch, ilev, iregion, ivar) )
+
+         analy%totspread(  iepoch, ilev, iregion, ivar) = &
+    sqrt(analy%totspread(  iepoch, ilev, iregion, ivar) / &
+         analy%Nused(      iepoch, ilev, iregion, ivar) )
+
    endif
 enddo
 enddo
@@ -1238,54 +1438,74 @@
 do iregion=1, Nregions
 do ilev=1, Nlevels
 
-   if (    guessAVG%Nused( ilev, iregion, ivar) == 0) then
-           guessAVG%rmse(  ilev, iregion, ivar) = MISSING_R4
-           guessAVG%bias(  ilev, iregion, ivar) = MISSING_R4
-           guessAVG%spread(ilev, iregion, ivar) = MISSING_R4
-        guessAVG%totspread(ilev, iregion, ivar) = MISSING_R4
+   if (    guessAVG%Nused(      ilev, iregion, ivar) == 0) then
+           guessAVG%observation(ilev, iregion, ivar) = MISSING_R4
+           guessAVG%ens_mean(   ilev, iregion, ivar) = MISSING_R4
+           guessAVG%bias(       ilev, iregion, ivar) = MISSING_R4
+           guessAVG%rmse(       ilev, iregion, ivar) = MISSING_R4
+           guessAVG%spread(     ilev, iregion, ivar) = MISSING_R4
+           guessAVG%totspread(  ilev, iregion, ivar) = MISSING_R4
 
    else
-           guessAVG%rmse(  ilev, iregion, ivar) = &
-      sqrt(guessAVG%rmse(  ilev, iregion, ivar) / &
-           guessAVG%Nused( ilev, iregion, ivar) )
+           guessAVG%observation(ilev, iregion, ivar) = &
+           guessAVG%observation(ilev, iregion, ivar) / &
+           guessAVG%Nused(      ilev, iregion, ivar)
 
-           guessAVG%bias(  ilev, iregion, ivar) = &
-           guessAVG%bias(  ilev, iregion, ivar) / &
-           guessAVG%Nused( ilev, iregion, ivar)
+           guessAVG%ens_mean( ilev, iregion, ivar) = &
+           guessAVG%ens_mean( ilev, iregion, ivar) / &
+           guessAVG%Nused(      ilev, iregion, ivar)
 
-           guessAVG%spread(ilev, iregion, ivar) = &
-      sqrt(guessAVG%spread(ilev, iregion, ivar) / &
-           guessAVG%Nused( ilev, iregion, ivar) )
+           guessAVG%bias(       ilev, iregion, ivar) = &
+           guessAVG%bias(       ilev, iregion, ivar) / &
+           guessAVG%Nused(      ilev, iregion, ivar)
 
-           guessAVG%totspread(ilev, iregion, ivar) = &
-      sqrt(guessAVG%totspread(ilev, iregion, ivar) / &
-           guessAVG%Nused(    ilev, iregion, ivar) )
+           guessAVG%rmse(       ilev, iregion, ivar) = &
+      sqrt(guessAVG%rmse(       ilev, iregion, ivar) / &
+           guessAVG%Nused(      ilev, iregion, ivar) )
 
+           guessAVG%spread(     ilev, iregion, ivar) = &
+      sqrt(guessAVG%spread(     ilev, iregion, ivar) / &
+           guessAVG%Nused(      ilev, iregion, ivar) )
+
+           guessAVG%totspread(  ilev, iregion, ivar) = &
+      sqrt(guessAVG%totspread(  ilev, iregion, ivar) / &
+           guessAVG%Nused(      ilev, iregion, ivar) )
+
    endif
 
-   if (    analyAVG%Nused( ilev, iregion, ivar) == 0) then
-           analyAVG%rmse(  ilev, iregion, ivar) = MISSING_R4
-           analyAVG%bias(  ilev, iregion, ivar) = MISSING_R4
-           analyAVG%spread(ilev, iregion, ivar) = MISSING_R4
-        analyAVG%totspread(ilev, iregion, ivar) = MISSING_R4
+   if (    analyAVG%Nused(      ilev, iregion, ivar) == 0) then
+           analyAVG%observation(ilev, iregion, ivar) = MISSING_R4
+           analyAVG%ens_mean(   ilev, iregion, ivar) = MISSING_R4
+           analyAVG%bias(       ilev, iregion, ivar) = MISSING_R4
+           analyAVG%rmse(       ilev, iregion, ivar) = MISSING_R4
+           analyAVG%spread(     ilev, iregion, ivar) = MISSING_R4
+           analyAVG%totspread(  ilev, iregion, ivar) = MISSING_R4
 

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list