No commit message

No commit message

da/cosmo/icbc4ctdas.py

@@ -130,4 +130,4 @@ def ct(date, co2_bg):
             oco2.units = 'kg kg-1'
             oco2[:] = co2_ens[e,ti,:,:,:]
 
-    ofile.close()
+        ofile.close()

da/cosmo/observationoperator.py

@@ -12,6 +12,7 @@ from datetime import datetime, timedelta
 from dateutil import rrule
 from cdo import *
 from . import site_height
+from da.cosmo.icbc4ctdas import ct
 from itertools import repeat
 from multiprocessing import Pool
 from da.tools.general import to_datetime
@@ -123,18 +124,19 @@ class ObservationOperator(object):
 
  # UNCOMMENT FROM HERE
 
+        co2_bg = np.empty(self.forecast_nmembers)
+
         for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start']+timedelta(hours=24*lag*int(dacycle['time.cycle'])), until=dacycle['time.start']+timedelta(hours=(lag+1)*24*int(dacycle['time.cycle']))):
             logging.info('Multiplying emissions with parameters for lag %d, date %s' % (lag, dt.strftime('%Y%m%d%H')))
             for ens in range(0,self.forecast_nmembers):
                 dthh = dt.strftime('%H')
+                co2_bg[ens] = members[ens].param_values[-1]
                 ens = str(ens).zfill(3)
                 cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['restartmap.dir'],"parameters_gpp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
                 cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['restartmap.dir'],"parameters_resp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
-#            co2_bg=str((members[ens].param_values[-1]))
-            co2_bg=str((members[:].param_values[-1]))
+            logging.info('Background CO2 params are (%s)' % co2_bg)
             if dthh=='00':
-                icbc4ctdas(dt.strftime('%Y%m%d'), co2_bg)
-                    #cdo.setunit("'kg m-2 s-1' -expr,BG_"+ens+"=CO2_BG*"+co2_bg+" -merge "+os.path.join(dacycle['da.bg.input'], 'ct_%s.nc' % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bg.input'], 'ensemble', "ct_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
+                ct(dt.strftime('%Y%m%d'), co2_bg)
 
             cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H')))
             cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H')))