diff --git a/graphics.py b/graphics.py
index b5e4ee19977d12fbba8992f01ba92c9fe8dc4598..2a600d36e2509202d84397603064ef6676e37405 100644
--- a/graphics.py
+++ b/graphics.py
@@ -508,7 +508,10 @@ def plot_free_ensemble(cbl,show='spread',ax=None,cmap='viridis'):
 
     if ax is not None:
         # Make plots
-        c = ax.pcolormesh(times/3600.,zt,bkgd,cmap=cmap,norm=colors.CenteredNorm())
+        if show=='cov' or show=='corr':
+            c = ax.pcolormesh(times/3600.,zt,bkgd,cmap=cmap,norm=colors.CenteredNorm())
+        else:
+            c = ax.pcolormesh(times/3600.,zt,bkgd,cmap=cmap)
         ax.contour(times/3600.,zt,theta.mean(axis=0),
                     np.linspace(cbl.theta_0,cbl.theta_0+cbl.gamma*zmax,ncont),
                     colors=contour_color,
@@ -666,18 +669,34 @@ def plot_desroziers(dgs,filename='desroziers_diagnostics.png',labels=None,logsca
         fig.set_size_inches(4,3)
         ax = fig.add_subplot(1,1,1)
 
+    # Check if all experiments have the same obs error variance
+    all_varo_equal = False
+    varo =[]
+    for dg in dgs:
+        varo.append(dg.varo)
+    varo = np.array(varo)
+    if np.all(np.isclose(varo, varo[0,0])):
+        varo = varo[0,0]
+        all_varo_equal = True
+
+
     colors = p.rcParams['axes.prop_cycle'].by_key()['color'][:len(dgs)]
     for i,dg,color in zip( range(len(dgs)) , dgs , colors):
+        # Select time interval to consider
+        idx = np.arange(0,dg.ominb.shape[0]) # whole run
+        #idx = np.arange(0,24) # first two hours
         if labels is not None:
             # Two formulations of the same thing; they give the same numbers!
             #ax1.scatter( dg.dobdob_ta-dg.mean_ominb**2 , dg.varo+dg.varb_ta , s=5, label = labels[i])
-            ax.scatter( np.mean(dg.sigma_ominb**2) , np.mean(dg.varo+dg.varb_ta) , s=13, label = labels[i],zorder=10, c=color, edgecolors='black')
-            ax.scatter( dg.sigma_ominb**2 , dg.varo+dg.varb_ta , s=10, alpha = 0.3, c=color, linewidths=0)
+            ax.scatter( np.mean(dg.ominb[idx,:].std(axis=0)**2) , np.mean(varo+(dg.sigmab[idx,:]**2).mean(axis=0)) , s=13, label = labels[i],zorder=10, c=color, edgecolors='black')
+            ax.scatter( dg.ominb[idx,:].std(axis=0)**2 , varo+(dg.sigmab[idx,:]**2).mean(axis=0) , s=10, alpha = 0.3, c=color, linewidths=0)
         else:
             # Two formulations of the same thing; they give the same numbers!
             #ax1.scatter( dg.dobdob_ta-dg.mean_ominb**2 , dg.varo+dg.varb_ta , s=5)
-            ax.scatter( np.mean(dg.sigma_ominb**2) , np.mean(dg.varo+dg.varb_ta) , s=13, zorder = 10)
-            ax.scatter( dg.sigma_ominb**2 , dg.varo+dg.varb_ta , s=10, alpha = 0.3) # 
+            ax.scatter( np.mean(dg.ominb[idx,:].std(axis=0)**2) , np.mean(varo+(dg.sigmab[idx,:]**2).mean(axis=0)) , s=13, zorder = 10)
+            ax.scatter( dg.ominb[idx,:].std(axis=0)**2 , varo+(dg.sigmab[idx,:]**2).mean(axis=0) , s=10, alpha = 0.3) # 
+        if all_varo_equal:
+            ax.axhline( y=varo, dashes=[3,1], color='lightgrey')
         #ax.set_title("Consistency of innovations",fontsize=8)
         ax.set_xlabel(r'${\langle}\sigma_d^2{\rangle}$ (K)') # '${\langle}d_{ob}^2{\rangle}_t$'
         ax.set_ylabel(r'$\sigma_{y^o}^2+{\langle}\sigma_{y^b}^2{\rangle}$ (K)')