Added up and downwind stations; use datetime instead of index for times

da/ffdas/observationoperator.py

@@ -57,7 +57,6 @@ d13C             = -9.
 Nav              = 6.023e23                                                       # molecules/mole Avogadro's number
 tau_14CO2        = 5700 * 365*24*3600.                                            # s              half life time of 14CO2
 lambda_14CO2     = np.log(2)/tau_14CO2                                            # s-1            radioactive decay
-tau_14CO2        = 5700 * 365*24*3600.                                            # s              half life time of 14CO2
 
 MOL2UMOL = 1.0e6
 KMOL2UMOL= 1.0e9
@@ -180,7 +179,8 @@ class STILTObservationOperator(ObservationOperator):
             if line[0] == '#': continue
             else:
                 ct_, filename, lat, lon, height, species_name, species_mass, recalc_factor, *_ = line.split(',')
-                sitename = filename.split('_')[1]
+                two_names = any(x in filename for x in ['UW', 'DW'])
+                sitename = ('_'.join(filename.split('_')[1:2 + two_names]))
                 ct = int(ct_) - 1 # Set the counter
                 self.obsnc.append(filename)
                 self.sitenames.append(sitename)
@@ -206,7 +206,6 @@ class STILTObservationOperator(ObservationOperator):
         self.ops_sector = []
         self.ops_id = []
         self.temporal_var = []
-        ct = 0
         for line in lines:
             if line[0] == '#': continue
             else:
@@ -215,7 +214,6 @@ class STILTObservationOperator(ObservationOperator):
                 if id != 0:
                     self.ops_sector.append(ct)
                     self.ops_id.append(id)
-                ct += 1
                 self.temporal_var.append(temporal_var)
 
         #set time control variables for this cycle
@@ -244,6 +242,18 @@ class STILTObservationOperator(ObservationOperator):
             self.datelist.append(dumdate)
             dumdate=dumdate+dt
 
+
+        self.get_c14_time_profile()
+
+    def get_c14_time_profile():
+        """Function that loads the nuclear power temporal variation"""
+        # Get the time profiles of the nuclear plants
+        temp_distribution_file = self.model_settings['nuclear_timeprofiledir']
+        workbook      = xlrd.open_workbook(temp_distribution_file)
+        worksheet     = workbook.sheet_by_name('Sheet1')
+        time_profile_nuc = np.array(worksheet.col_values(19)
+        self.nuc_time_profile = time_profile_nuc
+
     ### !!!! Only cache if made dependent on datepoint/datetime !!!!
     def get_time_index_nc(self, time=None):
         """Function that gets the time index from the flux files
@@ -332,12 +342,12 @@ class STILTObservationOperator(ObservationOperator):
         return background_conc
     
     @cached
-    def get_background_orig(self, i_species, i_datepoint):
+    def get_background_orig(self, i_species, datepoint):
         """Function that finds the background concentration, non-time dependent and hard-coded.
         Input:
             i_species: int: the index of the species for which the background should be found
             i_loc    : int: the index of the locatin for which the background should be found
-            datepoint: datetime.datetime: the datetime of the background concentration
+            datepoint: datetime for which the background concentration should be found.
         Returns:
             float: background concentration
             """
@@ -346,12 +356,12 @@ class STILTObservationOperator(ObservationOperator):
         return backgrounds[i_species]
 
     @cached
-    def get_c14_concentration(self, i_loc, i_datepoint,  gpp, ter, ff_flux, background):
+    def get_c14_concentration(self, i_loc, datepoint,  gpp, ter, ff_flux, background):
         """Function that gets the c14 concentration based on the emissions by the biosphere, 
             fossil fuels and nuclear power. The constants are defined at the top of this script.
         Input: 
             i_loc: int: the index of the location for which the c14 concentration should be calculated
-            i_datepoint: int: the index of the datepoint for which the c14 should be calculated
+            datepoint: datepoint: datetime for which the c14 concentration should be calculated
             gpp: float: The gross primary production in umol/s
             ter: float: The total ecosystem respiration in umol/s
             ff_flux: float: The fossil fuel flux in umol/s
@@ -359,7 +369,6 @@ class STILTObservationOperator(ObservationOperator):
         Returns: 
             float: The C14 in ppm
             float: Delta(14C) """
-        datepoint = self.datelist[i_datepoint]
         # First, get the footprint
         site = self.sitenames[i_loc]
         foot = self.get_foot(i_loc, datepoint)
@@ -367,16 +376,12 @@ class STILTObservationOperator(ObservationOperator):
         indices = self.get_time_indices(datepoint)
         
         # Get the fluxes from the nuclear plants
-        # Get the time profiles of the nuclear plants
-        temp_distribution_file = self.model_settings['nuclear_timeprofiledir']
-        workbook      = xlrd.open_workbook(temp_distribution_file)
-        worksheet     = workbook.sheet_by_name('Sheet1')
-        time_profile_nuc = np.array(worksheet.col_values(19,start_rowx=indices.start, end_rowx = indices.stop ))
+        nuclear_time_profile = self.nuc_time_profile[indices]
         # Get the fluxes and multiply them by the footprint and time profile
         # Note that the unit is already in Delta notation
         file = self.model_settings['nuclear_fluxdir']
         nuc_flux = self.get_nc_variable(file, 'E_14CO2_nuc')
-        nuc_flux = (time_profile_nuc[:, None, None] * nuc_flux[None,:] * foot).sum()
+        nuc_flux = (nuclear_time_profile[:, None, None] * nuc_flux[None,:] * foot).sum()
 
         # now get the disequilibrium flux
         file = self.model_settings['biosphere_fluxdir']
@@ -392,11 +397,11 @@ class STILTObservationOperator(ObservationOperator):
         co2_total = gpp + ter + ff_flux + background
 
 
-        print('ff', co2_14_ff * 1E12)
-        print('bg', co2_14_bg * 1E12)
-        print('gpp', alpha_14CO2_gpp * R_14CO2_bg * gpp * 1E12)
-        print('ter', R_14CO2_ter * ter)
-        print('nuc', nuc_flux)
+#        print('ff', co2_14_ff * 1E12)
+#        print('bg', co2_14_bg * 1E12)
+#        print('gpp', alpha_14CO2_gpp * R_14CO2_bg * gpp * 1E12)
+#        print('ter', R_14CO2_ter * ter)
+#        print('nuc', nuc_flux)
         # Conversion 14CO2 in ppm to DELTA ----
         Rm           = (co2_14_total * MASS_14CO2) / (co2_total * MASS_CO2)   # ppm to mass ratio (kg 14CO2 / kgCO2)
         A14          = Rm * lambda_14CO2 * 1e3*Nav / MASS_14CO2 # kg14CO2/kgCO2 * 1/s * molecules/kmole / (kg14CO2/kmole14CO2) = molecules 14CO2 / kgCO2 / s 
@@ -451,17 +456,16 @@ class STILTObservationOperator(ObservationOperator):
         if total: return gpp_increase + ter_increase
         else: return gpp_increase, ter_increase
 
-    @cached
-    def get_temporal_profiles(self, i_datepoint, i_station,  i_member):
+    def get_temporal_profiles(self, datepoint, i_station,  i_member):
         """ Function that reads in the temporal profiles for the current timestep
         Input:
-            i_datepoint: int: the index of the datepoint for which the c14 should be calculated
+            datepoint: datepoint: datetime for which the temporal profile should be found
             i_station: int: the index of the location for which the c14 concentration should be calculated
             i_member: int: the index of the member for which the simulation is run
         Returns:
             np.array (2): the time profiles for all categories. Indices: time, category"""
         #read temporal profiles for the times within the footprint
-        time_indices = self.get_time_indices(self.datelist[i_datepoint])
+        time_indices = self.get_time_indices(datepoint)
         station_name = self.sitenames[i_station]
         temporal_prior = os.path.join(self.model_settings['datadir'],'prior_temporal_{0}_{1:03}.nc'.format(station_name, i_member))
         temporal_prior = io.ct_read(temporal_prior, 'read')
@@ -490,10 +494,10 @@ class STILTObservationOperator(ObservationOperator):
         return emis
 
     @cached
-    def run_STILT(self, i_datepoint,i_member,i_loc, i_species, do_pseudo=0):
+    def run_STILT(self, datepoint, i_member, i_loc, i_species, do_pseudo=0):
         """This function reads in STILT footprints and returns hourly concentration data
         Input:
-            i_datepoint: int: index of the datepoint from self.datelist
+            datepoint: datepoint: datetime: the datepoint for which the concentrations should be calculated
             i_member   : int: index of the ensemblemember
             i_loc      : int: index of the location
             i_species  : int: index of the species
@@ -501,9 +505,7 @@ class STILTObservationOperator(ObservationOperator):
         Returns:
             float: the concentration increase at the respective location due to the emissions from the respective species"""
         # get the date:
-        datepoint = self.datelist[i_datepoint]
-
-        temp_profiles = self.get_temporal_profiles(i_datepoint, i_loc,  i_member)
+        temp_profiles = self.get_temporal_profiles(datepoint, i_loc,  i_member)
         spatial_emissions = self.get_spatial_emissions(i_member, i_species)
         #find correct background concentrationa
         footprint = self.get_foot(i_loc, datepoint)
@@ -545,11 +547,11 @@ class STILTObservationOperator(ObservationOperator):
                     logging.debug('... and species {}'.format(self.spname[i_species]))
                     for i_datepoint, datepoint in enumerate(self.datelist):
                         logging.debug('working on time {}'.format(datepoint))
-                        background_concentration = self.get_background_orig(i_species, i_datepoint)
-                        conc_increase_stilt = self.run_STILT(i_datepoint, i_member, i_loc, i_species)
+                        background_concentration = self.get_background_orig(i_species, datepoint)
+                        conc_increase_stilt = self.run_STILT(datepoint, i_member, i_loc, i_species)
                         if self.spname[i_species].upper() == 'CO2':
                             gpp, ter = self.get_biosphere_concentration(i_loc, datepoint)
-                            logging.debug('C14 (delta) = {}'.format(self.get_c14_concentration(i_loc, i_datepoint, gpp, ter, conc_increase_stilt, background_concentration)))
+#                            logging.debug('C14 (delta) = {}'.format(self.get_c14_concentration(i_loc, datepoint, gpp, ter, conc_increase_stilt, background_concentration)))
                         logging.debug('{0:.3} for FF and {1:.3} for bio for {2}'.format(conc_increase_stilt, gpp+ter, self.spname[i_species]))
                         conc_STILT.append(conc_increase_stilt + gpp + ter + background_concentration)
             totser[i_member] = np.array(conc_STILT)