Index: trunk/tools/validation_sim.py
===================================================================
--- trunk/tools/validation_sim.py	(revision 501)
+++ trunk/tools/validation_sim.py	(revision 502)
@@ -1268,4 +1268,13 @@
 (opts, args) = parser.parse_args()
 
+####### ###### ##### #### ### ## #
+# Number of different statistics according to the temporal coincidence
+#  0: Exact time
+#  1: Simulation values between consecutive observed times
+Nstsim = 2
+
+stdescsim = ['E', 'B']
+Lstdescsim = ['exact time', 'between observational time-steps']
+
 #######    #######
 ## MAIN
@@ -1594,4 +1603,5 @@
 newdim = onewnc.createDimension('stats',5)
 newdim = onewnc.createDimension('tstats',6)
+newdim = onewnc.createDimension('Nstsim, 2)
 
 # Variable dimensions
@@ -1618,4 +1628,9 @@
 basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
 writing_str_nc(newvar, ostatsn, StringLength)
+
+newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
+basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
+writing_str_nc(newvar, Lstdescsim, StringLength)
+
 
 if obskind == 'trajectory':
@@ -1671,8 +1686,8 @@
         oFillValue = None
 
-# Observed values temporally around coincident times
-#    simobsvalues, simobsSvalues, simobsTvalues, simobsTtvalues, trjsim =             \
-#        getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, ovobs,      \
-#        exacttvalues, oFillValue, Ngrid)
+# Observed values temporally exact times
+    Esimobsvalues, EsimobsSvalues, EsimobsTvalues, EsimobsTtvalues, trjsim =         \
+        getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, ovobs,      \
+        exacttvalues, oFillValue, Ngrid)
 
 # Observed values temporally around coincident times
@@ -1682,17 +1697,26 @@
 
 # Re-arranging values...
+    Earrayvals = np.array(Esimobsvalues)
     arrayvals = np.array(simobsvalues)
     if len(valvars[ivar]) > 2:
         const=np.float(valvars[ivar][3])
         if valvars[ivar][2] == 'sumc':
+            EsimobsSvalues = EsimobsSvalues + const
+            Earrayvals[:,0] = Earrayvals[:,0] + const
             simobsSvalues = simobsSvalues + const
             arrayvals[:,0] = arrayvals[:,0] + const
         elif valvars[ivar][2] == 'subc':
+            EsimobsSvalues = EsimobsSvalues - const
+            Earrayvals[:,0] = Earrayvals[:,0] - const
             simobsSvalues = simobsSvalues - const
             arrayvals[:,0] = arrayvals[:,0] - const
         elif valvars[ivar][2] == 'mulc':
+            EsimobsSvalues = EsimobsSvalues * const
+            Earrayvals[:,0] = Earrayvals[:,0] * const
             simobsSvalues = simobsSvalues * const
             arrayvals[:,0] = arrayvals[:,0] * const
         elif valvars[ivar][2] == 'divc':
+            EsimobsSvalues = EsimobsSvalues / const
+            Earrayvals[:,0] = Earrayvals[:,0] / const
             simobsSvalues = simobsSvalues / const
             arrayvals[:,0] = arrayvals[:,0] / const
@@ -1703,10 +1727,16 @@
 
 # statisics sim
-    simstats = np.zeros((5), dtype=np.float)
-    simstats[0] = np.min(arrayvals[:,0])
-    simstats[1] = np.max(arrayvals[:,0])
-    simstats[2] = np.mean(arrayvals[:,0])
-    simstats[3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
-    simstats[4] = np.sqrt(simstats[3] - simstats[2]*simstats[2])
+    simstats = np.zeros((2,5), dtype=np.float)
+    simstats[0,0] = np.min(Earrayvals[:,0])
+    simstats[0,1] = np.max(Earrayvals[:,0])
+    simstats[0,2] = np.mean(Earrayvals[:,0])
+    simstats[0,3] = np.mean(Earrayvals[:,0]*Earrayvals[:,0])
+    simstats[0,4] = np.sqrt(simstats[0,3] - simstats[0,2]*simstats[0,2])
+
+    simstats[1,0] = np.min(arrayvals[:,0])
+    simstats[1,1] = np.max(arrayvals[:,0])
+    simstats[1,2] = np.mean(arrayvals[:,0])
+    simstats[1,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
+    simstats[1,4] = np.sqrt(simstats[1,3] - simstats[1,2]*simstats[1,2])
 
 # statisics obs
@@ -1724,17 +1754,75 @@
  
 # Statistics sim-obs
-    simobsstats = np.zeros((9), dtype=np.float)
+    simobsstats = np.zeros((2,13), dtype=np.float)
+    Ediffvals = np.zeros((Nexactt), dtype=np.float)
     diffvals = np.zeros((Ncoindt), dtype=np.float)
 
+    Ediffvals = Earrayvals[:,0] - obsmask
     diffvals = arrayvals[:,0] - obsmask
 
-    simobsstats[0] = simstats[0] - obsstats[0]
-    simobsstats[1] = np.mean(arrayvals[:,0]*obsmask)
-    simobsstats[2] = np.min(diffvals)
-    simobsstats[3] = np.max(diffvals)
-    simobsstats[4] = np.mean(diffvals)
-    simobsstats[5] = np.mean(np.abs(diffvals))
-    simobsstats[6] = np.sqrt(np.mean(diffvals*diffvals))
-    simobsstats[7], simobsstats[8] = sts.mstats.pearsonr(arrayvals[:,0], arrayvals[:,1])
+    Ediff2vals = Ediffvals*Ediffvals
+    Esumdiff = Ediffvals.sum()
+    Esumdiff2 = Ediff2vals.sum()
+
+    diff2vals = diffvals*diffvals
+    sumdiff = diffvals.sum()
+    sumdiff2 = diff2vals.sum()
+
+    simobsstats[0,0] = simstats[0,0] - obsstats[0]
+    simobsstats[0,1] = np.mean(Earrayvals[:,0]*obsmask)
+    simobsstats[0,2] = Ediffvals.min()
+    simobsstats[0,3] = Ediffvals.max()
+    simobsstats[0,4] = Ediffvals.mean()
+    simobsstats[0,5] = np.abs(Ediffvals).mean()
+    simobsstats[0,6] = np.sqrt(Ediff2vals.mean())
+    simobsstats[0,7], simobsstats[0,8] = sts.mstats.pearsonr(Earrayvals[:,0],        \
+      obsmask)
+
+    simobsstats[1,0] = simstats[1,0] - obsstats[0]
+    simobsstats[1,1] = np.mean(arrayvals[:,0]*obsmask)
+    simobsstats[1,2] = diffvals.min()
+    simobsstats[1,3] = diffvals.max()
+    simobsstats[1,4] = diffvals.mean()
+    simobsstats[1,5] = np.abs(diffvals).mean()
+    simobsstats[1,6] = np.sqrt(diff2vals.mean())
+    simobsstats[1,7], simobsstats[1,8] = sts.mstats.pearsonr(arrayvals[:,0],         \
+      obsmask)
+# From:
+#Willmott, C. J. 1981. On the validation of models. Physical Geography, 2, 184â194
+#Willmott, C. J. (1984). On the evaluation of model performance in physical geography. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds., 443-460
+#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R. Legates, J. O'Donnell, and C. M. Rowe (1985), Statistics for the Evaluation and Comparison of Models, J. Geophys. Res., 90(C5), 8995-9005
+#Legates, D. R., and G. J. McCabe Jr. (1999), Evaluating the Use of "Goodness-of-Fit" Measures in Hydrologic and Hydroclimatic Model Validation, Water Resour. Res., 35(1), 233â241 
+#
+# Deviation of residuals (SDR)
+    simobsstats[0,9] = np.mean(np.sqrt(np.abs((Ediffvals-simobsstats[0,0])*(Ediffvals-\
+      obsstats[0]))))
+    simobsstats[1,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[1,0])*(diffvals- \
+      obsstats[0]))))
+# Index of Efficiency (IOE)
+    obsbias2series = (obsmask - obsstats[0])*(obsmask - obsstats[0])
+    sumobsbias2series = obsbias2series.sum()
+
+    simobsstats[1,10] = 1. - Esumdiff2/(sumobsbias2series)
+    simobsstats[1,10] = 1. - sumdiff2/(sumobsbias2series)
+# Index of Agreement (IOA)
+    Esimbias2series = Earrayvals[:,0] - obsstats[0]
+    simbias2series = arrayvals[:,0] - obsstats[0]
+
+    obsbias2series = obsmask - obsstats[0]
+
+    Eabssimbias2series = np.abs(Esimbias2series)
+    abssimbias2series = np.abs(simbias2series)
+    absobsbias2series = obbias2series.abs()
+
+    Eabssimobsbias2series = (Eabssimbias2series+absobsbias2series)*(Eabssimbias2series +\
+      absobsbias2series)
+    abssimobsbias2series = (abssimbias2series+absobsbias2series)*(abssimbias2series +\
+      absobsbias2series)
+
+    simobsstats[0,11] = 1. - Esumdiff2/(Eabssimobsbias2series.sum())
+    simobsstats[1,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
+# Fractional Mean Bias (FMB)
+    simobsstats[0,12] = (simstats[0,0]-obsstats[0])/(0.5*(simstats[0,0]+obsstats[0]))
+    simobsstats[1,12] = (simstats[1,0]-obsstats[0])/(0.5*(simstats[1,0]+obsstats[0]))
 
 # Statistics around sim values
@@ -1793,7 +1881,8 @@
 
 # sim Statistics
+    for ist in range(Nstsim):
     if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
-        newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('stats'),   \
-          fill_value=fillValueF)
+        newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('Nstsim',   \
+          'stats'), fill_value=fillValueF)
         descvar = 'simulated statisitcs: ' + valvars[ivar][0]
         basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, ovobs.getncattr('units'))
@@ -1811,5 +1900,5 @@
 # sim-obs Statistics
     if not searchInlist(onewnc.variables,varsimobs + 'st'):
-        newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('gstats'),            \
+        newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('Nstsim', 'gstats'),  \
           fill_value=fillValueF)
         descvar = 'simulated-observed statisitcs: ' + varsimobs
