Update outer boundary implementation to support open boundary conditions

[apcraig]

PR checklist

• Short (1 sentence) summary of your PR:
  Update outer boundary implementation to support open boundary conditions
• Developer(s):
  apcraig
• Suggest PR reviewers from list in the column to the right.
• Please copy the PR test results link or provide a summary of testing completed below.
  All tests are bit-for-bit except the halochk unittest which was refactored and passes. There is an issue in the evp1d implementation where on some compilers when the debug flag is on, the code will trap on a floating point error. See more details below. https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#85e4ed8a966197ba8479928cbd93ba8e0597df1a, https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#a8862c3bc36f31d3f5fa07d180f45e3a8e8979b6
• How much do the PR code changes differ from the unmodified code?
  • bit for bit
  • different at roundoff level
  • more substantial
• Does this PR create or have dependencies on Icepack or any other models?
  • Yes
  • No
• Does this PR update the Icepack submodule? If so, the Icepack submodule must point to a hash on Icepack's main branch.
  • Yes
  • No
• Does this PR add any new test cases?
  • Yes
  • No
• Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/. A test build of the technical docs will be performed as part of the PR testing.)
  • Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • No
• Please document the changes in detail, including why the changes are made. This will become part of the PR commit log.

Refactor CICE to support open boundary conditions where values may be imposed on the outer boundary. Prior to this update, the HaloUpdate was zeroing out the halo before updating the halo. This meant any values set on the halo were lost anytime a HaloUpdate was executed. Several other issues upstream and downstream were found as this was implemented. In particular, the global index set on the outer boundary for open and closed boundary conditions was "special", a zero global index was used internally as a flag, and both the scatter_global and haloUpdate was covering up some issues with uninitialized data.

Update the i_glob and j_glob indexing for open and closed boundary conditions on the outer boundary. For open, closed, and tripole (south) boundaries, now extrapolate index values on the outer boundary. This might mean indices are less than 0 or greater than n[x,y]_global+2*nghost. In the original implementation, these outer boundary indices were either set to some internal values so they would be set with the scatter_global method or to zero which was treated as a special value. The zero special value feature has been removed and the scatter_global implementation has been refactored to fill all gridpoints with valid global indices and leave others unset. The only special case remaining in the global indexing is for tripole grids on the north boundary where the j global index is set to a negative value.

• Refactor HaloUpdate in ice_boundary in serial and mpi to avoid zeroing out the halo. The ew and ns directions are handled separately. The new implementation will always zero the internal halo points. It will also zero the outer halo if cyclic or tripole is specified or if a fillvalue is passed into the haloUpdate. Otherwise, the outer halo is left unchanged.
  • Added ewBoundaryType and nsBoundaryType to halo datatype to keep track of the boundary conditions for each halo datatype.
  • Zero out the outer boundary during HaloUpdate only under certain conditions including cyclic or tripole bcs or if a fillValue is explicitly passed
• Add a few fillvalue=c0 arguments to HaloUpdate calls during initialization to force the outer boundary values to zero, mostly during initialization of grid data and similar cases.
• Add several HaloExtrapolate calls to initialize boundary values for grid fields. This replaces those fields being set by the incorrect global indices and the scatter.
• Explicitly zero out all allocated variables at initialization
• Explicitly zero out restart fields before reading to be extra careful
• Explicitly zero out some arrays in the evp1d before use
• Replace the G_HTN and G_HTE initialization with a call to gather_global_ext on HTN and HTE respectively after HTN and HTE are set and the halos updated.
• Move the closed boundary condition abort to initialization in ice_domain.F90
• Refactor halochk unit test, add explicit fill tests

Additional Features

• Update ice_HaloExtrapolate to support more than 1 ghost cell
• Add support for kmt_type = none with in subroutine rectgrid for rectangular grids. This allows a user to setup a rectangular grid with no land.
• Update the restart_ext implementation in ice_read_write. There was a bug if restart_ext=.false. was passed, nothing would happen. The update also allows some extra restart_ext logic in io_netcdf/ice_restart.F90 to be cleaned up.
• Carry out some code cleanup in init_domain_distribution and init_grid2
• Rename iglb/jglb to isrc/jsrc consistent in scatter_global_ext to be consistent with other scatter methods.

All tests are bit-for-bit except the halochk unittest which was refactored and passes. There is an issue in the evp1d implementation where on some compilers when the debug flag is on, the code will trap on a floating point error. If the floating point error isn't trapped, the cases run to completion and are bit-for-bit with the current trunk. This suggests there is something like a 0/0 calculation being trapped in the evp1d which has no impact on the solution, is only caught by a subset of compilers, and probably requires an if test to avoid. I spent some time trying to track down the error without any success.

[dabail10]

This is a significant amount of work and the code looks sensible. I'm confused though why everything should be bfb? Are we just not doing the right tests to capture this behavior?

[TillRasmussen]

This is necessary because the G_ARRAY of gather_global_ext is intent(inout).
Are there places where it is intent in?

[apcraig]

This may not be necessary in practice if we can be sure that the gather_global_ext sets all gridcells in the G_ arrays. But the cost of zeroing the arrays at initialization is also very small and probably reasonable to do. The hope is that zeroing these arrays makes the implementation more robust. Given we have a small technical issue in the evp1d implementation where some compilers trap a floating point error, this seems like a sensible thing to do for now.

[apcraig]

A slightly bigger worry is resetting all the G_ arrays to zero in convert_1d_2d_dyn below which happens each timestep. Again, it may improve robustness in the short term. But at some point, we may decide it's unnecessary and want to remove that code again.

[TillRasmussen]

Many of these conversions are not really needed or only needed for history output which should be done smarter than it currently is (only convert when needed).

[apcraig]

This is a significant amount of work and the code looks sensible. I'm confused though why everything should be bfb? Are we just not doing the right tests to capture this behavior?

It's bfb because open boundary conditions don't work except when there is land or no ice on the boundary. And all of our tests with open boundary conditions meet that criteria (global grids at j=1, box grids, etc). So, open boundary conditions weren't working and we weren't testing them. But now we need them to work. As we create cases with sea ice on an open boundary, we will probably have more debugging to do. This initial goal was just to refactor stuff to get ready for open boundary conditions and check whether answers changed. I hope to start work on another PR that would implement zero-gradient, constant-gradient, or other "technical" boundary conditions on open boundaries and create box test cases for them. That could be done concurrently with the open boundary condition work where specified fluxes/states are imposed on the boundary based on model/observations. The open boundary conditions create a lot of complexity though in terms of when and how to impose them in the timestepping.

[TillRasmussen]

I have looked through and the changes seems reasonably. The changes I have noted
Ice_haloupdate. Only fill halo if requested
ice_scatter has been simplified
ice_blocks only work global indices if cyclic.
ice_domain remove ocean count and abort if closed boundaries. I am not sure why on the latter.

In addition some semantics has been changed.

With regards to the 1d evp solver I will take a look at it along with the conversions of 1d to 2d data and vice versa. I should probably take a look at the halo updates of 1dEVP as well now that they have changed.

[anton-seaice]

Suggested change

	logical (kind=log_kind) :: lrestart_ext ! local value of reatart_ext
	logical (kind=log_kind) :: lrestart_ext ! local value of restart_ext