Extracting Regional Data from ERA5
To me I am not sure why the subset option should not work with a python script. According to me it does work. Here is a sample script.

import cdsapi

c = cdsapi.Client()

c.retrieve(
'reanalysis-era5-single-levels',
{
    'variable':'total_precipitation',
    'product_type':'reanalysis',
    'year':'2010',
    'month':'04',
    'day':'07',
    'area':'60/0/0/100',
    'time':[
        '00:00','01:00','02:00',
        '03:00','04:00','05:00',
        '06:00','07:00','08:00',
        '09:00','10:00','11:00',
        '12:00','13:00','14:00',
        '15:00','16:00','17:00',
        '18:00','19:00','20:00',
        '21:00','22:00','23:00'
    ],
    'format':'netcdf'
},
'precip.nc')


UPDATE
In response to some of the comments raised by @Nemesi - with this CDS Python API 
There is no need to use CDO or NCO to subset the data. The API itself takes care of it for you using the keyword area. Once you specify the bounding box the data that goes get downloaded is the bounding box boundaries you specify. There is also NO programmatic subsetting of global data. Right now the CDS web interface does not allow you to do subsetting(the problem is being fixed as I write according to Copernicus support). But that is not really a problem. Since OP's requirement is one of automating the process the CDS Python Web API fills the purpose. In this context subset is the same as ROI or region of interest. 

Here is what I had to do get this to work. I had to install the cdsapi for python. One can do this using pip (I usually do it - python3.6 setup.py install since I do not have conda). Then create a .cdsapirc(on Linux and other UNIX flavors it should be under the main home directory. Create this file wherever you have the $HOME variable defined) which should be like this - 

url: https://cds.climate.copernicus.eu/api/v2
key: {UID}:{API key}
verify:0


The values of UID and API key you should get when you register at the new CDS web interface site. 

And then just run it. 

Here the key parameter is the  area

and the first two values are starting latitude and longitude followed by the ending latitude and  longitude. When I viewed the downloaded netCDF file using ncdump it did have the right bounding box values. Now you could change it to GRIB format if you so wish but I think the basic format of the script does not change. 

Here is a plot of the subset of the precipitation for a specific instant in a 24 hour period.

Here is the code that plots it (I use matplotlib and cartopy)

from netCDF4 import Dataset,num2date
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.util import add_cyclic_point
import matplotlib as mpl
mpl.rcParams['mathtext.default'] = 'regular'
import matplotlib.pyplot as plt
import numpy as np


file = "precip.nc"
nc_pvFile = Dataset(file,'r')
lats = nc_pvFile.variables['latitude'][:]  # extract/copy the data
lons = nc_pvFile.variables['longitude'][:]
lats = lats[:].squeeze()
lons = lons[:].squeeze()

preciPlot = nc_pvFile.variables['tp'][:]

pp = preciPlot[0,:,:]


ax1 = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
clevs = np.arange(min(pp.flatten()),max(pp.flatten())*1000,1)


shear_fill = ax1.contourf(lons,lats,pp*1000,clevs,
                      transform=ccrs.PlateCarree(),cmap=plt.get_cmap('hsv'),
                                                                      linewidth=(10,),levels=100,extend='both')

ax1.coastlines()
ax1.gridlines()
ax1.set_xticks([0, 10,20,30,40,50,60,70,80,90,100], crs=ccrs.PlateCarree())
ax1.set_yticks([0, 10,20,30,40,50,60], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,
                               number_format='.0f')
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
cbar = plt.colorbar(shear_fill, orientation='horizontal')
plt.title('Total precipitation', fontsize=16)
plt.savefig('precip_era.png')