Stokes velocity fills dry cell uphill

[stef]

I have tried ROMS+SWAN coupling with wetting/drying and the model blew up.

jcwarner already had told me to use the COAWST version for ROMS+SWAN coupling, but I wasn't able to compile it (see below for the error. Also, I'm not sure where to report errors.). So I kept working with the Rutgers version.

It turned out that a dry cell kept exporting water until the depth became negative. I looked for differences between the Rutgers and COAWST versions of step2d_LF_AM3.h, and saw that the COAWST version applies the "export inhibitor" not only to the momentum equation (as Rutgers does), but also to the free surface equation, and there only for the wind-wave-induced velocities. I assume that if this is not done, the Stokes velocities continue to export water from a dry cell.

However, while trying to figure this out, I noticed that Stokes velocities are able to *import* water to a dry cell, even if that cell is located at a much higher elevation. I have one initially dry cell at h=-5 that constantly imports water from an adjecent wet cell (where zeta roughly 0). At some point the high and dry cell turns wet, and discharges the water to the surrounding cells (probably that's just a normal gravity wave, since it's higher than it's surroundings.)

Did I interpret this correctly, and if so, could it be solved by allowing the Stokes velocities only to contribute to the total transport velocity if u/vbar are not identical to zero, or perhaps if u/vbar are not disallowed to export? Not sure if this makes sense.

Errors from COAWST svn revision 1402:

I use gfortran and the following directives:

```
#define ROMS_MODEL
#define SWAN_MODEL

#define UV_LDRAG
#define ANA_INITIAL
#define MASKING
#define CURVGRID
#define WET_DRY
#define ATM_PRESS
#define NEARSHORE_MELLOR08
#define UV_ADV
#define WET_DRY_LIMIT

#define MCT_LIB
```

Building errors with
```
main2d.f90:170:50:

& MOD(iic(ng)-1,CoupleSteps(Iwaves,ng)).eq.0) THEN
1
Error: Symbol ‘iwaves’ at (1) has no IMPLICIT type; did you mean ‘iname’?
main2d.f90:172:50:

CALL ocn2wav_coupling (ng, tile)
1
Error: Missing actual argument for argument ‘tile’ at (1)
```

In main2d.F, mod_coupler is not used because COAWST_MODEL is defined. I couldn't find this switch in the coawst manual, so I simply commented it out in globaldefs.h

Building errors with
```
main2d.f90:173:50:

CALL ocn2wav_coupling (ng, tile)
1
Error: Missing actual argument for argument ‘tile’ at (1)
```

I guess this is because in mct_roms_swan.h, the subroutine is defined with 3 arguments like
```
SUBROUTINE ocn2wav_coupling (ng, iw, tile)
```
but I haven't figured out yet what `iw` means.

[jcwarner]

Thank you for the post. A couple of things:
- looks like you are compiling for just a 2D application? because the error was in main2d. I have not tried to run a 2d application in a long time. i really suggest that you do SOLVE3D. Many of the physics we have added are for 3D applications, like turbulence closures, sed transport, etc.
- i also notice that you are using the mellor approach. i recommend that you switch to use vortex force. we have not seen any new updates to the mellor approach in a long time, and the community is leaning towards the dissipation (vf) method.
- suggest you do not just turn things off in global_defs.h that is a sensitive routine, with things listed in certain orders. if you go to a 3d applications, and compile INLET_TEST, that should compile and run.
- for the stokes transport uphill into a dry cell: that does not sound correct. I do believe we have approaches to prevent that, but if it is happening then i need to dig in deep and correct it. maybe you could try with 3d, VF, and see if that problem still exists.

-john

[stef]

Thanks for your reply!

Yes, I should have mentioned it's a 2d application. I'm trying to simulate storm surges. At some point I want to have at least a bottom boundary layer, but having 2d simulations for reference and comparison is interesting. (also I just noticed that I have the WET_DRY_LIMIT directive listed above, which I believe does nothing in 2d anyways, but I have to check)

Yes, I used the Mellor approach although you already had advised me not to do it in a previous thread. That's only because I had trouble compiling the COAWST version.

I'll follow your advice and configure a COAWST 3d application with the vortex method next week, and report back about the uphill transport.

[stef]

I tested the spurious wetting of high elevation grid points now with the 3d version of COAWST. No spurious wetting occurred in my test, neither using WEC_VF nor WEC_MELLOR.

Browsing through the source code and Warner et al. (2013) I currently think the following (unverified/untested hypothesis):

In a SOLVE3D setup, the radiation stress is computed once per baroclinic timestep. u/vbar_stokes stay constant during the barotropic timesteps. They are vertical integrals of u/v_stokes, so the "binary state" u/vmask_wet used in the radiation stress routines for u/v_stokes carry over to the barotropic velocity. This effectively means that u/vbar_stokes do not participate in the wetting (I mean at the wet/dry boundaries).

As described in Warner et al. (2013), in a SOLVE3D setup, u/vmask_wet is computed as the average over the barotropic timesteps (at the end of the barotropic stepping). This is what I mean by "binary state" (0 or 1) vs. "multi-valued state" (0 or 1 or -1 or 2) as during the barotropic stepping.

In a 2d setup, the "binary state" of u/vmask_wet is not used at all. The "multi-valued state" is used to inhibit drying of dry cells below DCRIT (same procedure as in SOLVE3D) in the momentum equation. In particular, u/vbar are set to zero at "critical" dry/wet boundaries (i.e. where water is exported from already dry cells). This carries over to the next free surface equation, but there, critical export by u/vbar_stokes *alone* is not necessarily inhibited (contrary to the SOLVE3D case).

Regarding the adaptions of the free surface stepping in COAST vs. Rutgers (described in the post above), I'm now confused as to why they are there. Critical export should ideally be inhibited based on the total flow, but u/vbar components at the free-surface stepping had already been set to zero in the momentum stepping in the previous iteration of step2d. So it's not possible to compute the total flow at the free-surface stepping any more. Considering u/vbar_stokes in isolation at the free surface stepping can lead to cases where u/vbar_stokes has the "correct" sign (importing into a dry cell), but is possibly added to a zero u/vbar that was artificially set zero because of gross exporting. The hypothesis is that this leads to spurious wetting.

I'm not sure if anyone is actually interested in 2D simulations with radiation stress (here on the ROMS forum), but if so I appreciate any feedback. Anyways I'm going to do more tests. I'll try to compute the total ("Lagrangian") barotropic transport and apply the inhibitor based on it.