diff --git a/da/cosmo/base_optimizer.py b/da/cosmo/base_optimizer.py
index 595299a35549c440ac48f73621e94445f13f7101..a2788037722267484d60cbbb250b569688dcaac4 100755
--- a/da/cosmo/base_optimizer.py
+++ b/da/cosmo/base_optimizer.py
@@ -50,7 +50,7 @@ class Optimizer(object):
self.nlag = dims[0]
self.nmembers = dims[1]
# self.nparams = dims[2]
- self.nparams = 22
+ self.nparams = 23
self.nobs = dims[3]
self.create_matrices()
diff --git a/da/cosmo/observationoperator.py b/da/cosmo/observationoperator.py
index a29805ddb7bd8442e33ded832f969bfd8e0e43e8..d6f21d4723d5400c4bacc51e7b815897985ad6af 100755
--- a/da/cosmo/observationoperator.py
+++ b/da/cosmo/observationoperator.py
@@ -124,44 +124,42 @@ 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')))
- logging.info('Background CO2 params are (%s)' % co2_bg)
- if dthh=='00':
- 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')))
-# if dthh=='00' or dthh=='03' or dthh=='06' or dthh=='09' or dthh=='12' or dthh=='15' or dthh=='18' or dthh=='21':
- # cdo.merge(input = os.path.join(dacycle['da.bg.input'], 'ensemble', "ct_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bg.input'], 'ensemble', "ct_%s.nc" % dt.strftime('%Y%m%d%H')))
-
- os.chdir(dacycle['da.obsoperator.home'])
-
- if os.path.exists("/scratch/snx3000/parsenov/ctdas/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
- os.rename("/scratch/snx3000/parsenov/ctdas/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), "/scratch/snx3000/parsenov/ctdas/non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
-
- os.system('python run_chain.py ctdas '+absolute_start_time+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
- os.chdir(dacycle['dir.da_run'])
-
- args = [
- (dacycle, hstart, hstop, self.forecast_nmembers)
- for dacycle, (hstart, hstop), self.forecast_nmembers
- in zip(repeat(dacycle),
- [(starth+168*lag,endh+168*lag-1)],
- repeat(self.forecast_nmembers))
- ]
-
- with Pool(3) as pool:
- pool.starmap(self.extract_model_data, args)
-
+# 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')))
+ # logging.info('Background CO2 params are (%s)' % co2_bg)
+ # if dthh=='00':
+ # 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')))
+#
+ # os.chdir(dacycle['da.obsoperator.home'])
+#
+ # if os.path.exists("/scratch/snx3000/parsenov/ctdas/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
+ # os.rename("/scratch/snx3000/parsenov/ctdas/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), "/scratch/snx3000/parsenov/ctdas/non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
+#
+ # os.system('python run_chain.py ctdas '+absolute_start_time+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
+ # os.chdir(dacycle['dir.da_run'])
+#
+ # args = [
+ # (dacycle, hstart, hstop, self.forecast_nmembers)
+ # for dacycle, (hstart, hstop), self.forecast_nmembers
+ # in zip(repeat(dacycle),
+ # [(starth+168*lag,endh+168*lag-1)],
+ # repeat(self.forecast_nmembers))
+ # ]
+#
+ # with Pool(3) as pool:
+ # pool.starmap(self.extract_model_data, args)
+#
for i in range(0,self.forecast_nmembers):
idx = str(i).zfill(3)
cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
@@ -224,7 +222,7 @@ class ObservationOperator(object):
for s,ss in enumerate(sites):
site_height.main(cosmo_out, str(ens), ss, time_stamp)
- cdo.intlevel("860", input = cosmo_out+"CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = cosmo_out+"/modelled_"+ens+"_lhw_"+time_stamp+".nc")
+ cdo.intlevel("860", input = cosmo_out+"CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc")
cdo.intlevel("797", input = cosmo_out+"CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc")
cdo.intlevel("3580", input = cosmo_out+"CO2_60lev_"+ens+"_jfj_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc")
cdo.intlevel("1205", input = cosmo_out+"CO2_60lev_"+ens+"_ssl_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc")
diff --git a/da/cosmo/statevector.py b/da/cosmo/statevector.py
index f472ad345496fd16f0d3f40f758641785a7a5794..230c5d48650966e647b0aa5fa5bbfa696b97be21 100644
--- a/da/cosmo/statevector.py
+++ b/da/cosmo/statevector.py
@@ -292,18 +292,24 @@ class StateVector(object):
# Create members 1:nmembers and add to ensemble_members list
for member in range(1, self.nmembers):
- rands = np.random.uniform(low=-1., high=1., size=self.nparams-1)
- rands_bg = np.random.uniform(low=-0.05, high=0.05, size=1)
-# if (member==1):
- # rands = np.array([1.836451,1.91572,1.719371,1.947495,1.788597,1.109507,0.9401712,2,0.003201536,1.2276,0])-1.
- # elif (member==2):
- # rands = np.array([2,0.04773442,1.809479,1.136452,2,0.8753765,1.488298,2,1.908293,0.02784704,0])-1.
- # elif (member==3):
- # rands = np.array([0.9136019,1.386204,2,0.1250395,1.815127,2,0.3002012,1.666945,0,0,0])-1.
-
+# rands = np.random.uniform(low=-1., high=1., size=self.nparams-1)
+ # rands_bg = np.random.uniform(low=-0.05, high=0.05, size=1)
+ if member==1 and lag==0:
+ rands = np.array([0.9251788, 0.9034393, 0.336998, 0.9244774, -0.05096465, 1.039527, 0.8610526, -0.6545241, -0.151757, 0.6358872, -1, 0.4625102, 0.1186881, 0.6585811, 1.080793, 0.7087865, 0.1603414, 0.7908162, 0.9887264, -0.5482104, 0.8389163, -1, 0.02753718 ])
+ elif member==2 and lag==0:
+ rands = np.array([0.6557752, -0.1970569, 0.04818814, -0.688437, -0.7006997, -0.1784677, 0.2049817, 0.2842229, -0.656664, -0.4288206, -1, -0.4610326, 0.7538341, -0.4886873, -0.715229, 0.3835333, -0.8013363, 0.1057041, 0.4560323, 0.4176753, -0.2784889, -1, -0.04719913])
+ elif member==3 and lag==0:
+ rands = np.array([0.01135811, -0.01033817, 0.8016006, 0.2941527, -0.3657097, 0.1622472, -0.14451, -0.9826015, -0.3264703, 0.5378729, -1, 0.936269, 0.2612394, 0.8938854, 0.3610448, 1.104215, -0.1021504, 0.2102389, 0.803875, -0.4403061, -0.8015866, -1, -0.04969293])
+
+ if member==1 and lag==1:
+ rands = np.array([0.5848618, -0.263966, -0.8484242, -0.4811731, 0.008701986, -0.5143792, 0.196903, 0.9151311, -0.3057676, 1.025769, -1, 0.6145338, -0.109208, 0.9366637, -0.3343047, -0.5159711, -0.7831933, -0.114835, -0.1131854, -0.9170142, 0.01517715, -1, 0.04161017])
+ if member==2 and lag==1:
+ rands = np.array([0.2462473, -0.4267351, -0.7604715, 0.3000261, 0.690653, 0.5368793, 0.617731, 0.1198399, -0.5313886, 0.1087431, -1, 0.6121695, -0.02330998, -0.03271778, 0.2461925, -0.00954848, 0.1422206, -0.3688108, -0.357688, -0.474921, -0.8561779, -1, 0.02497981])
+ if member==3 and lag==1:
+ rands = np.array([0.6437422, -0.6527606, 0.4303354, 0.5737818, 0.4010006, 0.8635789, -0.6543487, 0.4480304, 0.4105155, 0.7875805, -1, -0.7137064, 0.1760372, 0.7394882, -0.4832366, 0.6765335, 0.4377979, -0.09937549, -0.75212, -0.9841197, 0.0008806029, -1, 0.02755764])
newmember = EnsembleMember(member)
-# newmember.param_values = np.concatenate([np.dot(C, rands[0:self.nparams/2]),np.dot(C, rands[self.nparams/2:self.nparams])]) + newmean
- newmember.param_values = (np.hstack((np.dot(C, rands[0:11]),np.dot(C, rands[11:22]), rands_bg)) + newmean).ravel()
+ # newmember.param_values = (np.hstack((np.dot(C, rands[0:11]),np.dot(C, rands[11:22]), rands_bg)) + newmean).ravel()
+ newmember.param_values = rands+1.
newmember.param_values[10] = 0.
newmember.param_values[21] = 0.
newmember.param_values[newmember.param_values<0.] = 0.