Small fix related to total domain error

pycif/plugins/controlvects/standard/build_hcorr.py

@@ -62,6 +62,7 @@ def build_hcorrelations(zlat, zlon, lsm,
     try:
         evalues, evectors = \
             read_hcorr(nlon, nlat, sigma_sea, sigma_land, dir_dump, hresoldim=hresoldim)
+        check.verbose("Reading hcorr from file")
         
     
     # Else build correlations
@@ -105,23 +106,25 @@ def build_hcorrelations(zlat, zlon, lsm,
         corr = np.exp(old_div(-dx, sigma))
         corr[sigma <= 0] = 0
 
-        # Scale covariance in accordance with global total error
 
-        # Prior error averaged over time steps
-        xerr = np.mean(np.reshape(cntrlv.std, (tracer.ndates,
-                                int(len(cntrlv.std)/tracer.ndates))), axis=0)
-        
-        # cov = np.outer(cntrlv.std, cntrlv.std) * corr
-        cov = np.outer(xerr, xerr) * corr
-        covsum = cov*np.outer(cntrlv.area_reg, cntrlv.area_reg)
-        toterr = np.sqrt(covsum.sum())*3600.*24.*365./1.e9/float(tracer.numscale)
-        errscalar = glob_err/toterr
-        verbose("Total error scaled by "+ str(errscalar))
-        verbose("covsum "+ str(covsum.sum()))
-
-        # ESO TODO: try scaling std instead 
-        # corr = corr*errscalar**2
-        cntrlv.std = cntrlv.std*errscalar**2
+        if glob_err is not None:
+            # Scale covariance in accordance with global total error
+
+            # Prior error averaged over time steps
+            xerr = np.mean(np.reshape(cntrlv.std, (tracer.ndates,
+                                    int(len(cntrlv.std)/tracer.ndates))), axis=0)
+
+            # cov = np.outer(cntrlv.std, cntrlv.std) * corr
+            cov = np.outer(xerr, xerr) * corr
+            covsum = cov*np.outer(cntrlv.area_reg, cntrlv.area_reg)
+            toterr = np.sqrt(covsum.sum())*3600.*24.*365./1.e9/float(tracer.numscale)
+            errscalar = glob_err/toterr
+            verbose("Total error scaled by "+ str(errscalar))
+            cntrlv.std = cntrlv.std*errscalar
+
+            # ESO: this is done in flexinvert / correl.f90
+            corr[corr < 0.0001] = 0.0
+                           
         
         # Component analysis
         evalues, evectors = np.linalg.eigh(corr)
@@ -143,8 +146,8 @@ def build_hcorrelations(zlat, zlon, lsm,
     
     # Truncating values < evalmin
     # mask = evalues >= evalmin
+    
     # ESO: This is how it is done in flexinvert
-
     mask = evalues >= evalmin*evalues.max()
     check.verbose("Truncating eigenvalues at " +  str(evalmin*evalues.max()))
     
@@ -189,7 +192,7 @@ def read_hcorr(nlon, nlat, sigma_sea, sigma_land, dir_dump, hresoldim=None):
         nlon, nlat (ints): dimensions of the domain
         sigma_land, sigma_sea (floats): horizontal correlation distances
         dir_dump (str): where the horizontal correlations have been stored
-        hresoldim: if not None, indicates numer of regions used
+        hresoldim: if not None, indicates number of regions used
     
     """
 

pycif/plugins/controlvects/standard/init_bprod.py

@@ -91,6 +91,7 @@ def init_bprod(cntrlv, options={}, **kwargs):
 
             elif hasattr(tracer, 'hcorrelations') and tracer.hresol == 'regions':
                 # Scale by domain total error (given as Tg/y)
+                glob_err=None
                 if hasattr(tracer, 'glob_err'):
                     glob_err = float(tracer.glob_err)
                     if glob_err < 0: