docs: def generate_seasonal_averages(qout_file, seasonal_average_file,
                                 num_cpus=multiprocessing.cpu_count()):
      """
      This function loops through a CF compliant rapid streamflow
      file to produce a netCDF file with a seasonal average for
      365 days a year
      """
      with RAPIDDataset(qout_file) as qout_nc_file:
          print("Generating seasonal average file ...")
          seasonal_avg_nc = Dataset(seasonal_average_file, 'w')
  
          seasonal_avg_nc.createDimension('rivid', qout_nc_file.size_river_id)
          seasonal_avg_nc.createDimension('day_of_year', 365)
  
          time_series_var = seasonal_avg_nc.createVariable('rivid', 'i4',
                                                           ('rivid',))
          time_series_var.long_name = (
              'unique identifier for each river reach')
  
          average_flow_var = \
              seasonal_avg_nc.createVariable('average_flow', 'f8',
                                             ('rivid', 'day_of_year'))
          average_flow_var.long_name = 'seasonal average streamflow'
          average_flow_var.units = 'm3/s'
  
          std_dev_flow_var = \
              seasonal_avg_nc.createVariable('std_dev_flow', 'f8',
                                             ('rivid', 'day_of_year'))
          std_dev_flow_var.long_name = 'seasonal std. dev. streamflow'
          std_dev_flow_var.units = 'm3/s'
  
          std_dev_flow_var = \
              seasonal_avg_nc.createVariable('max_flow', 'f8',
                                             ('rivid', 'day_of_year'))
          std_dev_flow_var.long_name = 'seasonal max streamflow'
          std_dev_flow_var.units = 'm3/s'
  
          std_dev_flow_var = \
              seasonal_avg_nc.createVariable('min_flow', 'f8',
                                             ('rivid', 'day_of_year'))
          std_dev_flow_var.long_name = 'seasonal min streamflow'
          std_dev_flow_var.units = 'm3/s'
  
          lat_var = seasonal_avg_nc.createVariable('lat', 'f8', ('rivid',),
                                                   fill_value=-9999.0)
  
          lon_var = seasonal_avg_nc.createVariable('lon', 'f8', ('rivid',),
                                                   fill_value=-9999.0)
          add_latlon_metadata(lat_var, lon_var)
  
          seasonal_avg_nc.variables['lat'][:] = \
              qout_nc_file.qout_nc.variables['lat'][:]
          seasonal_avg_nc.variables['lon'][:] = \
              qout_nc_file.qout_nc.variables['lon'][:]
  
          river_id_list = qout_nc_file.get_river_id_array()
          seasonal_avg_nc.variables['rivid'][:] = river_id_list
          seasonal_avg_nc.close()
  
      # generate multiprocessing jobs
      mp_lock = multiprocessing.Manager().Lock()  # pylint: disable=no-member
      job_combinations = []
      for day_of_year in range(1, 366):
          job_combinations.append((qout_file,
                                   seasonal_average_file,
                                   day_of_year,
                                   mp_lock
                                   ))
  
      pool = multiprocessing.Pool(num_cpus)
      pool.map(generate_single_seasonal_average,
               job_combinations)
      pool.close()
      pool.join()