Let grdconvert write coregistered grids into a single netCDF file

**Description of the desired feature**

While **grdinterpolate** and **greenspline** can create 3-D netCDF data cubes _v_(_x,y,z_), we have no way to create a single netCDF file with multiple 2-D grids _z1_(_x,y_), z2(_x,y_), etc., where these grids all shar the same (_x,y_) lattice. GMT can already read these as lots of Earth and Ocean science data are distributed with multiple layers (temperature, salinity in oceanography, P and S velocities in seismology, etc). It would be convenient for GMT users to simply place related grids in a single netCDF file instead of learning how to use CDO, NCO and other tools (which can do much more but requires installation, learning, etc.). Hence, I propose a simple extension to **grdconvert**:

**gmt grdconvert -L** _grid1 grid2 ... gridn_ **-G**_outfile.nc_

where the different data sharing the same **-R -I -r** are simply stored in outfile.nc. Other options like **-R** for selecting subsets) remain.

A related issue that may become important to address as well is the naming of single netCDF layers. Many GMT grid creators (e.g., **surface**, **sphinterpolate**) simply calls the layer "z", while others (**grdinterpolate**, **greenspline**) has the ability to name the layer via the (relatively new) **-D+v**_name_ modifier. Clearly, we cannot place layers with the same names in **grdconvert**. So in addition to the above capability, perhaps

1. Make a simple way to name an output grid variable [z] from all modules that create grids, maybe as a modifier to the grid name (e.g., `-Gbathymetry.grd+vdepth`). If this was implemented we could deprecate the recent **+v**_name_.
2. Let **grdconvert -L** take arguments that provide the desired names (e.g., `grdconvert -Lvs,vp vs.grd vp.grd -Gseis_vel.nc`)

If anyone has other thoughts on this, please comment, @GenericMappingTools/core and @GenericMappingTools/gmt-contributors.