The code now tests for the presence of the optimizer file, which is False for...

The code now tests for the presence of the optimizer file, which is False for the purely forward runs. If not present, some steps related to the forecast variables are skipped.

da/analysis/expand_mixingratios.py

@@ -37,6 +37,7 @@ def WriteMixingRatios(DaCycle):
     from string import join
     import shutil
     import logging
+    import netCDF4
 
        
     dirname     = 'data_molefractions'
@@ -57,7 +58,7 @@ def WriteMixingRatios(DaCycle):
 
     # Step (1): Get the posterior sample output data file for this cycle
 
-    infile   = os.path.join(DaCycle['dir.output'],'sampleinfo_%s.nc'% DaCycle['time.sample.stamp'])
+    infile   = os.path.join(DaCycle['dir.output'],'sampleinfo_%s__newstyle.nc'% DaCycle['time.sample.stamp'])
 
     ncf_in   = io.CT_Read(infile,'read')
 
@@ -65,8 +66,9 @@ def WriteMixingRatios(DaCycle):
     obs_val   = ncf_in.GetVariable('observed')
     simulated = ncf_in.GetVariable('modelsamples')
     infilename= ncf_in.GetVariable('inputfilename')
+    infiles   = netCDF4.chartostring(infilename).tolist()
 
-    infiles   = [join(s.compressed(),'') for s in infilename]
+    #infiles   = [join(s.compressed(),'') for s in infilename]
 
     ncf_in.close()
 
@@ -74,25 +76,33 @@ def WriteMixingRatios(DaCycle):
 
     infile   = os.path.join(DaCycle['dir.output'],'optimizer.%s.nc'% startdate.strftime('%Y%m%d') )
 
-    ncf_fc_in = io.CT_Read(infile,'read')
+    if os.path.exists(infile): 
+        optimized_present = True
+    else:
+        optimized_present = False
 
-    fc_obs_num   = ncf_fc_in.GetVariable('obspack_num')
-    fc_obs_val   = ncf_fc_in.GetVariable('observed')
-    fc_simulated = ncf_fc_in.GetVariable('modelsamplesmean_prior')
-    fc_simulated_ens = ncf_fc_in.GetVariable('modelsamplesdeviations_prior')
-    fc_flag      = ncf_fc_in.GetVariable('flag')
-    fc_r         = ncf_fc_in.GetVariable('modeldatamismatchvariance').diagonal()
-    fc_hphtr     = ncf_fc_in.GetVariable('totalmolefractionvariance').diagonal()
-    filesitecode = ncf_fc_in.GetVariable('sitecode')
+    if optimized_present:
 
-    fc_sitecodes = [join(s.compressed(),'') for s in filesitecode]
+        ncf_fc_in = io.CT_Read(infile,'read')
 
-    ncf_fc_in.close()
+        fc_obs_num   = ncf_fc_in.GetVariable('obspack_num')
+        fc_obs_val   = ncf_fc_in.GetVariable('observed')
+        fc_simulated = ncf_fc_in.GetVariable('modelsamplesmean_prior')
+        fc_simulated_ens = ncf_fc_in.GetVariable('modelsamplesdeviations_prior')
+        fc_flag      = ncf_fc_in.GetVariable('flag')
+        fc_r         = ncf_fc_in.GetVariable('modeldatamismatchvariance').diagonal()
+        fc_hphtr     = ncf_fc_in.GetVariable('totalmolefractionvariance').diagonal()
+        filesitecode = ncf_fc_in.GetVariable('sitecode')
 
-    # Expand the list of input files with those available from the forecast list
+        fc_sitecodes  = netCDF4.chartostring(filesitecode).tolist()
+        #fc_sitecodes = [join(s.compressed(),'') for s in filesitecode]
 
-    infiles_rootdir = os.path.split(infiles[0])[0]
-    infiles.extend(os.path.join(infiles_rootdir,f+'.nc') for f in fc_sitecodes)
+        ncf_fc_in.close()
+
+        # Expand the list of input files with those available from the forecast list
+
+        infiles_rootdir = os.path.split(infiles[0])[0]
+        infiles.extend(os.path.join(infiles_rootdir,f+'.nc') for f in fc_sitecodes)
 
 
     #Step (2): For each observation timeseries we now have data for, open it and fill with data
@@ -240,7 +250,7 @@ def WriteMixingRatios(DaCycle):
         # Get all obs_nums related to this file, determine their indices in the local arrays
 
         selected_obs_nums = [num for infile,num in zip(infiles,obs_num) if infile == orig_file]
-        selected_fc_obs_nums = [num for sitecode,num in zip(fc_sitecodes,fc_obs_num) if sitecode in orig_file]
+
 
         # Optimized data 1st: For each index, get the data and add to the file in the proper file index location
 
@@ -263,28 +273,32 @@ def WriteMixingRatios(DaCycle):
 
         # Now forecast data too: For each index, get the data and add to the file in the proper file index location
 
-        for num in selected_fc_obs_nums:
+        if optimized_present:
 
-            model_index = fc_obs_num.tolist().index(num)
-            file_index  = file_obs_nums.tolist().index(num)
+            selected_fc_obs_nums = [num for sitecode,num in zip(fc_sitecodes,fc_obs_num) if sitecode in orig_file]
+
+            for num in selected_fc_obs_nums:
+
+                model_index = fc_obs_num.tolist().index(num)
+                file_index  = file_obs_nums.tolist().index(num)
 
-            var = ncf_out.variables['modeldatamismatch']
-            var[file_index] = np.sqrt(fc_r[model_index])
+                var = ncf_out.variables['modeldatamismatch']
+                var[file_index] = np.sqrt(fc_r[model_index])
 
-            var = ncf_out.variables['modelsamplesmean_forecast']
-            var[file_index] = fc_simulated[model_index]
+                var = ncf_out.variables['modelsamplesmean_forecast']
+                var[file_index] = fc_simulated[model_index]
 
-            var = ncf_out.variables['modelsamplesstandarddeviation_forecast']
-            var[file_index] = fc_simulated_ens[model_index,1:].std()
+                var = ncf_out.variables['modelsamplesstandarddeviation_forecast']
+                var[file_index] = fc_simulated_ens[model_index,1:].std()
 
-            var = ncf_out.variables['modelsamplesensemble_forecast']
-            var[file_index] = fc_simulated_ens[model_index,:]
+                var = ncf_out.variables['modelsamplesensemble_forecast']
+                var[file_index] = fc_simulated_ens[model_index,:]
 
-            var = ncf_out.variables['totalmolefractionvariance_forecast']
-            var[file_index] = fc_hphtr[model_index]
+                var = ncf_out.variables['totalmolefractionvariance_forecast']
+                var[file_index] = fc_hphtr[model_index]
 
-            var = ncf_out.variables['flag_forecast']
-            var[file_index] = fc_flag[model_index]
+                var = ncf_out.variables['flag_forecast']
+                var[file_index] = fc_flag[model_index]
 
 
         # close the file