OBCS_SEAICE_NEUMANN flag, more store directives, and cleanup

[mjlosch]

What changes does this PR introduce?

new feature and some additional tweaks of TAF store directives and cleanup

What is the current behaviour?

• no Neumann boundary condition for sea ice variables
• defining ALLOW_OBCS_BALANCE and ALLOW_OBCS_STEVENS leads to (hidden) recomputations and also to some recomputation warnings (untested configuration)

What is the new behaviour

• Neumann boundary conditions for sea ice variables, activate by defining new flag OBCS_SEAICE_NEUMANN, defining this flag will override any other input from files for sea ice variables.
• fewer (hidden) recomputations and no recomputation warnings when both ALLOW_OBCS_BALANCE and ALLOW_OBCS_STEVENS and OBCS_SEAICE_NEUMANN are defined
• new runtime flag for pkg/cost: cost_mask_file, which can be useful for flexible cost function definitions

Does this PR introduce a breaking change?

No

Other information:

These are changes that were required for a specific configuration. I think they are useful for everyone.

Suggested addition to tag-index

o pkg/obcs:

• Neumann boundary conditions for sea ice variables, activate by defining new flag OBCS_SEAICE_NEUMANN
• fewer (hidden) recomputations and no recomputation warnings when both ALLOW_OBCS_BALANCE and ALLOW_OBCS_STEVENS and OBCS_SEAICE_NEUMANN are defined
  o pkg/cost: new runtime flag for pkg/cost: cost_mask_file, which can be useful for flexible cost function definitions

[jm-c]

should this be "= uFld(i,Jobc-1,bi,bj)" instead ?

[jm-c]

I took a look at code changes and have few comments:

1. There are many TAF store dir changes, and few changes to testreport AD output (all 4 global_ocean.90x40x15 tests + global_ocean_ebm + tutorial_global_oce_optim), all small differences and some in forward grad-check and some in AD grad-check + AD monitor. Will need to run this on villon to check if I need to update few ref-output. Also need to run with Tapenade since I notices few changes related to Tapenade code.
2. there is a new CPP "OBCS_SEAICE_NEUMANN" in SEAICE_OPTIONS.h but I see few "#ifdef SEAICE_OBCS_NEUMANN" in obcs_apply_uvice.F & obcs_apply_seaice.F.
3. Also,since it does not involve new arrays in common block, would it be difficult to have a run-time switch instead of a CPP option ?
4. At some point, it would be good to get a minimal doc of pkg/obcs implemetation in the manual (currently 8.3.1.6 and 8.3.1.5.10 are empty), specially for seaice since it does not follow the same sequence (e.g., obcs_apply_seaice called at the very end, after all pkg/seaice update but obcs_apply_uvice called before any update of pkg/seaice fields ).
5. Not sure it would be better, but could we skip the changes in obcs_apply_uvice.F and set OB[N,S,E,W]u,vice instead ?
6. very minor, for "cost config" report in STDOUT, I like what pkg/ctrl reports:

control vector configuration >>> START <<<
control vector configuration >>> END <<<

and could use the same convention here.

[jm-c]

Detail: I don't see anymore any change related to "divisions by zero in taf-generated ad-code of ggl90_calc.F", so may be the PR description could be updated.

[mjlosch]

@jm-c thanks for having a look and for you comments

1. There are many TAF store dir changes, and few changes to testreport AD output (all 4 global_ocean.90x40x15 tests + global_ocean_ebm + tutorial_global_oce_optim), all small differences and some in forward grad-check and some in AD grad-check + AD monitor. Will need to run this on villon to check if I need to update few ref-output. Also need to run with Tapenade since I notices few changes related to Tapenade code.

This is not the cleanest PR, because it mixes the Neumann BCs for seaice with other stuff. The "other" stuff is more reorganisation of store directives to avoid extra computations where they are not necessary (e.g. 85ec739. This reorders the AD code and I can imagine that there are small changes due to reading tapes instead of recomputing (truncation). Not sure where the FW differences come from. I will have another look to make sure that everything is OK with that.
update: On my MacBook with gfortran, where results are not perfect w.r.t. the reference to begin with, I have only these differences (for the three experiments listed above):

< Y Y Y Y 14>12< 7 13 12 13 13 14 16 16 16 16 16 14 14 14 16 16 16 pass  global_ocean.90x40x15.bottomdrag (this branch)
---
> Y Y Y Y 14>12< 7 13 12 13 13 14 16 13 16 16 16 14 14 14 13 16 16 pass  global_ocean.90x40x15.bottomdrag (master)

i.e. the agreement with the reference even improves a little. This somewhat supports my feeling that the differences is a truncation error issue when replacing recomputation by tapes.

2. there is a new CPP "OBCS_SEAICE_NEUMANN" in SEAICE_OPTIONS.h but I see few "#ifdef SEAICE_OBCS_NEUMANN" in obcs_apply_uvice.F & obcs_apply_seaice.F.

That's embarrassing and will be fixed. It's a bad idea not to test new code ... will think of a verification experiment where I can add this functionality.

3. Also,since it does not involve new arrays in common block, would it be difficult to have a run-time switch instead of a CPP option ?

That should be possible. I didn't do it because I thought that extra "if" would make the code less efficient. Do think that this type of code would be OK? Would the runtime flag SEAICEuseNeumannBC be part of data.obcs or data.seaice

          IF ( Jobc.NE.OB_indexNone ) THEN
           IF ( SEAICEuseNeumannBC ) THEN
C     Neumann boundary conditions
            HEFF (i,Jobc,bi,bj) = HEFF (i,Jobc-1,bi,bj)
            AREA (i,Jobc,bi,bj) = AREA (i,Jobc-1,bi,bj)
            HSNOW(i,Jobc,bi,bj) = HSNOW(i,Jobc-1,bi,bj)
#  ifdef SEAICE_VARIABLE_SALINITY
            HSALT(i,Jobc,bi,bj) = HSALT(i,Jobc-1,bi,bj)
#  endif
           ELSE
            HEFF (i,Jobc,bi,bj) = OBNh (i,bi,bj)
            AREA (i,Jobc,bi,bj) = OBNa (i,bi,bj)
            HSNOW(i,Jobc,bi,bj) = OBNsn(i,bi,bj)
#  ifdef SEAICE_VARIABLE_SALINITY
            HSALT(i,Jobc,bi,bj) = OBNsl(i,bi,bj)
#  endif
           ENDIF
          ENDIF

4. At some point, it would be good to get a minimal doc of pkg/obcs implemetation in the manual (currently 8.3.1.6 and 8.3.1.5.10 are empty), specially for seaice since it does not follow the same sequence (e.g., obcs_apply_seaice called at the very end, after all pkg/seaice update but obcs_apply_uvice called before any update of pkg/seaice fields ).

That would be great. There are so many sea ice options that I did not implement, so that I am not sure that could describe them properly, but we could start doing so in this PR. I'll give it a shot.

5. Not sure it would be better, but could we skip the changes in obcs_apply_uvice.F and set OB[N,S,E,W]u,vice instead ?

Would be more efficient, see 3. point

6. very minor, for "cost config" report in STDOUT, I like what pkg/ctrl reports:

control vector configuration >>> START <<<
control vector configuration >>> END <<<

and could use the same convention here.

Yes, we could modify cost_readparms.F in this direction, will try.

[jm-c]

@mjlosch I have not looked at your latest changes but few comments regarding your answers:

1. since I wrote my last set of comments, I did run a Tapenade testreport (AD) and everything is good.

2. the run-time switch vs CPP option: Since there is already a "IF ( I/Jobc.NE.OB_indexNone ) THEN" inside the loop, with I/Jobc dependency on loop index, it might not make a big difference to have an other "IF / THEN / ELSE" inside the same loop. So the run-time option might be a good thing.

3. regaring point 5:

   Not sure it would be better, but could we skip the changes in obcs_apply_uvice.F and set OB[N,S,E,W]u,vice instead ?

   this could be left for later, because I am not sure it would be an improvement and also how this will work with future seaice + obcs development ? not clear to me.

I will take a look at the new docs (exiting !).

[mjlosch]

in d8ac131 I applied some sort of hybrid: Along with replacing the CPP flag with a runtime flag for Neumann BCs, I copied the corresponding uIce/vIce to the OB?u/vice within the routine obcs_apply_uvice.F. It's not ideal, but this block could be moved easily, but where? It would be logical to put it into obcs_calc.F, but after call obcs_prescribe_read, right?

[jm-c]

@mjlosch regarding this:

It would be logical to put it into obcs_calc.F, but after call obcs_prescribe_read, right?

Should this be after the OBCS_CONTROL ? or at the same level as OBCS_STEVENS ? (and in this case just after ?).

[mjlosch]

neither ctrl_getobcs nor obcs_calc_stevens does anything to sea ice, so we are free to put in anywhere.

I guess putting this at the end after obcs_calc_stevens makes the most sense to me, if we agree that it is supposed to override any other boundary condition. I should create a regression test for this anyway before I start moving stuff around.

[mjlosch]

I now added a forward test of the sea ice Neumann boundary conditions to seaice_obcs.tides, which let me immediately catch two errors in the process, commit aa27cd6

Everything after commit 5c100af is preparation for a TAF-AD test for seaice+obcs+Neumann (+and aEVP). This new test required a few adjustments. Especially 031eb42 may appear as an overkill, but there were new RECOMPUTATION warnings in EVP and I used this to fix the EVP AD code. In the test it works nicely.

I had to turn off tides in AD because they lead to massive recomputation warnings which I couldn't fix quickly. That's something on the to-do list.

Both reference solutions are generated on my MacBook and will require updates from the reference machine.