No commit message

No commit message

da/cosmo/observationoperator.py

@@ -127,15 +127,15 @@ class ObservationOperator(object):
 
         co2_bg = np.empty(self.forecast_nmembers)
 
+        logging.info('Multiplying emissions with parameters for lag %d' % (lag))
         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)
+#            logging.debug('Background CO2 params are (%s)' % co2_bg)
             if dthh=='00':
                 ct(dt.strftime('%Y%m%d'), co2_bg)
 
@@ -145,7 +145,10 @@ class ObservationOperator(object):
         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))
+            if os.path.exists("/scratch/snx3000/parsenov/ctdas/non_opt_"+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/old_non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
+            else:
+                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_ch+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
         os.chdir(dacycle['dir.da_run'])

da/cosmo/pipeline.py

@@ -44,6 +44,7 @@ def ensemble_smoother_pipeline(dacycle, platform, dasystem, samples, statevector
     
     invert(dacycle, statevector, optimizer)
     
+    sys.exit()
     advance(dacycle, samples, statevector, obsoperator)
         
     save_and_submit(dacycle, statevector)    

da/cosmo/statevector.py

@@ -292,24 +292,11 @@ 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 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])
+            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)
+
             newmember = EnsembleMember(member)
- #           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 = (np.hstack((np.dot(C, rands[0:11]),np.dot(C, rands[11:22]), rands_bg)) + newmean).ravel()
             newmember.param_values[10] = 0.
             newmember.param_values[21] = 0.
             newmember.param_values[newmember.param_values<0.] = 0.

template.py

@@ -77,7 +77,7 @@ ensemble_smoother_pipeline(dacycle, platform, dasystem, samples, statevector, ob
 ################### All done, extra stuff can be added next, such as analysis
 ##########################################################################################
 
-logging.info(header + "Starting analysis" + footer) 
+logging.info(header + "All done. Cheers" + footer) 
 sys.exit(0)
 
 save_weekly_avg_1x1_data(dacycle, statevector)

template.rc

@@ -113,7 +113,7 @@ da.obsoperator.rc      : ${da.obsoperator.home}/tm5-ctdas-ei-zoom.rc
 ! The number of ensemble members used in the experiment. Valid entries are integers > 2
 !
 
-da.optimizer.nmembers  : 4
+da.optimizer.nmembers  : 21
 
 ! Finally, info on the archive task, if any. Archive tasks are run after each cycle to ensure that the results of each cycle are
 ! preserved, even if you run on scratch space or a temporary disk. Since an experiment can take multiple weeks to complete, moving