getvar() from many wrf output files?

First and foremost, thanks for providing this fantastic tool.

I'm using wrf.getvar() to open a time series for several variables (HFX, LH, some of the diagnostic variables) that are stored in WRF-written netCDF files. Each file contains a single temporal value (for a 30-minute period). It's a 4-month WRF run, so there are many of these files (> 5000).

For shorter WRF runs with fewer files I've passed getvar() a list of netCDF4.Datset objects.

Now I'm hitting a limit for number of open files (OSError: [Errno 24] Too many open files).

Is there a "best practice" for reading a single variable from lots and lots of netCDF files? It seems that xarray (1) isn't yet supported for getvar() and (2) may not work well anyway because xarray.open_mfdataset seems to want to read every variable from each WRF file and is thus very slow.

My WRF files are netCDF4 (not netCDF4-classic) which seems to rule out netcdf4.MFDataset().

Is my best bet to use something like ncrcat to make a temporary netCDF file containing only the variable I want? This could work but would, I guess, require some digging to supply all the WRF output variables needed for some of the wrf.getvar() diagnostic variables.

The issue names getvar(), netCDF4.Dataset lists, xarray.open_mfdataset, MFDataset, and ncrcat, but no repository file or test. Start by inspecting getvar() handling of Dataset lists and reproducing the file-descriptor failure; done would be a documented supported workflow or a scoped change verified against thousands of files.