This PR addresses a few bugs and improves documentation:

1. Let grdflexure handle the combination of -H -W: If the load has variable load density then we must use the density at each node to compute the equivalent height for an average load density. Also, this must be done before we determine if any part of the load (now) sticks above sea level for the additional boost due to the greater density contrast with air.
2. The check for NaNs in the input grid gave an error right away, hence we never even got to the section where the NaNs were replaced by zeros. Now fixed.
3. Since this module as well as many others (gravfft, grdgravmag, grdvector) need to check if the two grids are truly co-registered I have added a new function gmt_grd_domains_match which consolidates the (many) checks into one function. These modules now use this function to simplify those checks.
4. The effects of -H and -W on the load heights were never discussed - now the documentation explains what is going on.

No changes to any tests.