Allow variable sea level

[rjleveque]

WIP - don't merge yet. I still plan to write a test that illustrates this, and we should discuss some design choices. Also I just realized there are several other routines that use sea_level where something different should be done when this is variable.

This adds a new subroutine src/2d/shallow/set_set_level.f90 that takes a grid patch and returns an array of the same shape that has the value of sea level to use in each cell, which is only used when filling a cell with fluid up to this level. This is done in setting the initial conditions (which I just realized needs to be modified) and also when a new level of AMR is introduced to initialize cells that cover a dry grid cell on the coarser level if the shoreline cuts through that dry cell. After refinement, some of the finer cells may be wet and are filled to the level determined by sea_level. With this modification, this value can vary spatially (and with time).

The use case this was developed for was to handle a coastal region where

1. Subsidence (or uplift) occurs due to an earthquake, as specified in the dtopo file,
2. The tsunami generated by the quake does not arrive until several hours after the quake, and so refinement to very fine grids (e.g. 1/3 arcsecond) does not need to be done until a couple hours into the simulation. This can happen if the location is in large shallow bay, e.g. Aberdeen, WA.

When the ground subsides, the water in the bay subsides too, so if refinement happens later the value of h should be chosen to bring the water up to a sea level that is lower than the original sea level by an amount given by interpolating the dtopo file to this location. This is what the set_sea_level function provided here computes. (Note this assumes no water has flowed in the intervening time near the coast, which is not generally true, but is a better assumption than filling these cells to the original sea level, which might immediately flood all subsided regions with water that had no opportunity to get there.)

For other applications, the user could provide their own version of set_sea_level. One possible use case is initializing a river, where the initial water level varies with the elevation of the river, or to fill onshore lakes near shore that have different water levels than the ocean's sea level.

Question:
Is sea_level the right term for these parameters, especially given the other possible use cases?
Would initial_surface_level or something be clearer?

Still needs fixing:

• qinit.f90 to set properly at initial time
• qinit_module.f90 also uses sea_level in initializing from a qinit file
• friction_module.f90 use sea_level but only in a test case, so it's ok?

[mandli]

This would also be useful for overland flooding situations and rivers (storm surge or just plain flood events). Also @aks2203 had implemented a more generic eta_init for the multilayer code as sea_level was one value and we needed multiple, might be an easy way to name this something other than sea level.

[rjleveque]

@mandli: I'm not sure why travis is failing, seems to be a netcdf problem...

======================================================================
ERROR: Test reading and writing of storm formats
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/virtualenv/python2.7.14/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/home/travis/build/clawpack/geoclaw/tests/test_storm.py", line 143, in test_storm_IO
    raise e
RuntimeError: NetCDF: Not a valid ID

[mandli]

Better solution for this incoming from @rjleveque.