[ENH] add eta to rbcs for nonlinear free surface

[jklymak]

What changes does this PR introduce?

Draft for comment: obcs allows setting etaH, but I wanted an rbcs-type sponge, so I added this to rbcs.

This is potentially valuable for fjords and other embayments where we want to force the tide via sea surface height in a large receiving basin. I do this by adjusting the sponge in the receiving basin so that the surface wave is free in the fjord itself.

Does this PR introduce a breaking change?

No

Other information:

I have a simple version of a toy model that works with this at https://github.com/jklymak/FjordForce. That example could be simplified, its just a simulation I had that had a reasonable fjord with large receiving basin set up.

• Needs documentation
• Needs tests?

[jklymak]

One thing I'm not clear on - for obcs they have a calculation for OBCS_APPLY_SURF_DR. I wasn't clear if this was needed here, or because I was directly manipulating etaH it would get calculated in CALC_SURF_DR.

[mjlosch]

docs compile (html), code compiles in the two verification experiments that use rbcs, but I get this (without the -devel option):

Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16  .  .  .  . pass  exp4
Y Y Y Y>--< 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  .  .  .  . N/O   exp4.nlfs
Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16  .  .  .  .  .  .  .  . pass  exp4.stevens
Y Y Y Y>12<16 16 16 16 16 16 16 16 16 16 16 16 16 16 14 16  .  .  .  .  .  .  .  . pass  tutorial_reentrant_channel

because results are NaNs already in timestep 0:

(PID.TID 0000.0001) %MON surfExpan_theta_mean         =                   NaN
(PID.TID 0000.0001) %MON surfExpan_salt_mean          =                   NaN

with the -devel option, the model does not run exp4.nlfs. I am guessing that there's a problem in the initialisation phase.

[mjlosch]

Maybe only print if the debug level is larger than default, like in e.g. model/src/external_fields_load.F?

[mjlosch]

something went wrong with the indentation.

Since there are so many stops in this file, maybe we could import the "errCount=errCount+1" method here, too? (see, e.g., model/src/config_check.F)

[jklymak]

Thanks @mjlosch for taking a look. For my testing, did you run exp4.nlfs as-is, or did you modify at all?

[mjlosch]

I just checked out the branch and ran ./testreport -t 'exp4.nlfs tutorial_reentrant_channel with and without the -devel option. I did this on my MacBook with gfortran. If you cannot reproduce that, I'll have a closer look.

[jklymak]

Thanks - I can have a look in a few days.

[mjlosch]

Here's the problem: in rbcs_update_etah.F there's a division by tauRelaxEta without checking if it is zero or if useRBCeta is true or not. I suggest to add an IF ( useRBCeta ) THEN in that routine (or set rec_tauRlx=0. if tauRelaxEta==0.)

[jklymak]

Thanks @mjlosch this passes that test report. I added an if statement to rbc_apply_eta.F.

Happy to write a test for this setup.

[mjlosch]

LGTM, thanks for adding some documenation. See my comment there.

I changed the error handling to print all errors before stopping; this has nothing to do with this PR.

[mjlosch]

I am afraid, that mentioning the shape may be confusing. I don't even know what our convention is. For python-generated fields it is clear, i.e. (nt,ny,nx), but with all gendata.m-Matlab codes, it's the other way around (nx,ny,nt), isn't it?

[jrscott]

Our convention is math notation :) and/or fortran, not python (nor explicitly matlab per se either).
But I would agree with Martin here you don't need to mention the shape here.

[jklymak]

I find it pretty confusing to know what dimension gets the time index, and have to look it up in one of the examples or my own past setups. Where would I find this elsewhere in the docs? I think it's a gap in the docs that dimensions are not specified for most i/o files, though nz, ny, nx are a more natural order. Given the number of folks who mess up their input fields in the mailing lists, I don't think I am alone in getting this wrong sometimes

[jrscott]

see section 3.9 in the doc

We've purposely tried to keep the doc from catering to either matlab or python (or netcdf). In the code, the order would be (nx, ny, nz) so that would be the natural order, as it were.

[jrscott]

If explicitly the placement of the time index is the confusion, we could add this case in the 3.9 doc section.

[jklymak]

OK, removed this, and slightly rephrased the relax*File entry to include relaxEtaFile, and added a link back to sec 3.9 (that is very clear, thanks!).

I'll submit a follow-up PR for section 3.9 based on this discussion.

[jm-c]

@jklymak thanks for sharing your code.

I was supposed to look at this PR a long time ago but a little busy these days.

There are several things that I would need to check (including how this fit into time-stepping with various options) as it affects the continuity equation and might affect tracer conservation.
This would not be an issue if the relaxation was only used near Open-Boundary (part of a "sponge layer") but might not always be the case.

The other thing is that we are considering an alternative implementation of Open-Boundary Conditions that would allow naturally to prescribe the SSH at the OB and account for the pressure gradient across the OB.
If this was available and working, it might be used instead of this Eta relaxation, at least in some cases.

So I would suggest we keep this PR open and keep it up-to-date (will merge latest master branch in) and decide later which direction to take.

[jklymak]

Thanks @jm-c - I use this in an a-physical sponge part of the domain where I would also sponge active and passive tracers, so I don't consider tracer conservation important in my use cases. I'm not sure use cases where a-physically forcing eta could conserve anything (definitely not mass!) but perhaps there are such use cases...