Ok, I will test it with GMTSAR but it will take a little while. Need to find & download new files from terrible Copernicus and next set up the scripts.

OK, GMTSAR script finished and results look good but since it was run in parallel mode no messages were printed on screen. What I can say is that at least it didn't screw.

Well, that is a start!

I think what we will find there is that this improvement will actually add time but it should improve convergence. I might possibly find time to do some testing on that. For instance, we have access to the exact solution that surface is trying to fit via greenspline. The only difference will be near the borders as the boundary conditions are different, but we could isolate a sub-region in the middle and do some rms comparisons between the surface solution (with or without -R smarts) and the truth from greenspline.

I ran the tests with this PR and the results for surface/periodic.sh etc are interesting. Much better result with this PR due to more iterations I think. The original PS file in baseline shows bad fits toward the northern boundary but the new ones are much closer to the input circles.

I think this PR has now been completed. To recap, we are breaking backwards compatibility by trying to provide the user with a better result due to more intermediate gridding intervals. It is just too tedious for users to

1. Run surface -Q to figure out what would have been a better -R -I settings
2. Run surface with those settings
3. Use grdcut to exact the desired region

This PR automates this process and makes it the default action. For experts who wish to retain the previous and less accurate solution they can use -Qr which will honor the -R exactly as given in the calculations.

Implementing this feature was relatively simple (we already had the -Q machinery and it was just a question of picking the best improved region and reset to that temporarily, then undo when we write out the result). However, a complication arose with -L where users may supply grids for a limiting upper or lower grid surface. Because these are expected to share the same -R as the desired region, these also need to be extended under the hood, and the extra rows and columns that are added will need to be initialized to NaN since we do not have any constraints for them. So that took a bit of fiddling to get right. I added the new test limits.sh earlier (in master) which for the first time tests the effect of the lower and upper limit for both a constant and a grid, and I do this on our favorite Table 5.11 with -R0/6.3/0/6.3 -I0.1 which yields a 63x63 system which then gets promoted to 64x64 for more intermediate stages. The output of all those tests looks correct, and slightly different.

However, the new convergence scheme means all tests and examples using surface now "fails". Turns out we only have two data sets being used across many tests so we only need to look at two results. I mentioned the Table 5-11 gridding. Here is the diff for one part of example 14:

Obviously it is not clear which solution is "best" in these cases. The new solution iterated a total of 1017 times while the previous (master) version iterated 855 times - so not a huge change for the 63->64 step. Yet, since we know that we keep iterating forever to get closer and closer to convergence I think it is reasonable to assume we want the new result.

However, it is a bit clearer I think with the other case in test/surface/periodic*.sh. Here I create a synthetic grid with a product of sin(lon) times cos(lat/3) so very predictable and smooth data. Below you see the diff and you will notice that near the northern boundary the old solution pulls up while the new solution has converged to a more horizontal trend (matching the synthetic data). Since the synthetic data was randomly sampled (see dots) we do not have lots of data up north hence the continued convergence seem to help quite a bit:

While the master solution would only go through intermediate spacings 10-2-1 for a total of 6642 iterations, the new version has intermediate stages 36-12-4-2-1 for 27610 iterations. Doing the diff grid we see the biggest improvements is along the boundaries where there is little constraints locally and we depend on gradients being propagated in from a distance:

So I think all is well and better. These are the failing scripts:

The following tests FAILED:
	211 - doc/examples/ex14/ex14.sh (Failed)
	213 - doc/examples/ex16/ex16.sh (Failed)
	414 - test/grdcontour/filled_cont.sh (Failed)
	428 - test/grdcontour/view_filled.sh (Failed)
	1131 - test/surface/limits.sh (Failed)
	1132 - test/surface/periodic.sh (Failed)
	1133 - test/surface/periodic_pix.sh (Failed)

so once we are OK with this PR I will need to update some PS files. Note: I called this branch surface-faster but what I meant was surface-better. It is not faster to iterate more, but doing so gives a better result and I think we want GMT to give good result first and foremost. Again, if this slows down Sandwell's gridding too much he can use -Qr or tighten other convergence criteria.