Time bounds returned after an operation with resample-method

Problem description

For datamining with xarray there is always the following issue with the resampling-method.
If i resample e.g. a timeseries with hourly values to monthly values, the netcdf-standards tell us to put into the result file information about:

1. the bounds for each timestep over which the aggregation was taken (for each month the beginning and the end of the month)
2. the method which was used for aggregation decoded by the variable attribute 'cell_method' (e.g. 'time: mean').

The recent implementation should be improved which is proven by the following data example.

Data example

I have a dataset with hourly values over a period of 5 month.

<xarray.Dataset>
Dimensions:       (bnds: 2, time: 3672)
Coordinates:
    rlon          float32 22.06
    rlat          float32 5.06
  * time          (time) datetime64[ns] 2006-05-01 2006-05-01T01:00:00 ...
Dimensions without coordinates: bnds
Data variables:
    rotated_pole  int32 1
    time_bnds     (time, bnds) float64 1.304e+07 1.305e+07 1.305e+07 ...
    TOT_PREC      (time) float64 nan nan nan nan nan nan nan nan nan nan nan ...
Attributes:

Doing a resample process using the mean operator gives

In [36]: frs
Out[36]: 
<xarray.Dataset>
Dimensions:       (bnds: 2, time: 5)
Coordinates:
  * time          (time) datetime64[ns] 2006-05-31 2006-06-30 2006-07-31 ...
Dimensions without coordinates: bnds
Data variables:
    rotated_pole  (time) float64 1.0 1.0 1.0 1.0 1.0
    time_bnds     (time, bnds) float64 1.438e+07 1.438e+07 1.702e+07 ...
    TOT_PREC      (time) float64 12.0 nan nan nan nan

Here the time_bnds is still in the file but the content is very strange:

In [37]: frs["time_bnds"]
Out[37]: 
<xarray.DataArray 'time_bnds' (time: 5, bnds: 2)>
array([[  1.438020e+07,   1.438380e+07],
       [  1.701540e+07,   1.701900e+07],
       [  1.965060e+07,   1.965420e+07],
       [  2.232900e+07,   2.233260e+07],
       [ -6.330338e+10,  -6.330338e+10]])
Coordinates:
  * time     (time) datetime64[ns] 2006-05-31 2006-06-30 2006-07-31 ...
Dimensions without coordinates: bnds

So, he still knows that time_bnds is related to the coordinate time. However, the values are not correct. The first time_bnds entry should be [1.5.2006 00:00,31.5.2006 23:00]. That is definitely not the case, i.e. the numbers here are related to the original file (seconds since 2005-12-01), but they do not match to my expection. 1.438020e+07 equals "Dienstag, 16. Mai 2006, 10:30:00" and 1.438380e+07 equals "Dienstag, 16. Mai 2006, 11:30:00".
Moreover, the xarray's do not consider to change the unit of the time_bnds according the unit of the variable 'time' if data is written to netcdf. Output of the program ncdump reveals that time was changed to days since but time_bnds seems to be still coded in "seconds since".

ncdump -v time_bnds try.nc
netcdf try {
dimensions:
	time = 5 ;
	bnds = 2 ;
variables:
	double rotated_pole(time) ;
		rotated_pole:_FillValue = NaN ;
	double time_bnds(time, bnds) ;
		time_bnds:_FillValue = NaN ;
	double TOT_PREC(time) ;
		TOT_PREC:_FillValue = NaN ;
	int64 time(time) ;
		time:units = "days since 2006-05-31 00:00:00" ;
		time:calendar = "proleptic_gregorian" ;
data:

 time_bnds =
  14380200, 14383800,
  17015400, 17019000,
  19650600, 19654200,
  22329000, 22332600,
  -63303379200, -63303379200 ;
}

Is there a recommendation what to do?

No source file, test, or entry point is named. Start by reproducing the hourly-to-monthly resampling example, then inspect how resampling produces time_bnds and how the result is written to NetCDF. Done means monthly bounds cover the aggregation periods, their units match time on output, and the aggregation method is represented as requested.