Update halo and interior restoring

[apcraig]

PR checklist

• Short (1 sentence) summary of your PR:
  Update halo and interior restoring
• Developer(s):
  isievers, 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.
  results are bit-for-bit except boxrestore tests have changed answers.
• How much do the PR code changes differ from the unmodified code?
  • bit for bit except boxrestore tests which have changed answers because the implementation changed, https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#5e4aed141c8b3b55bfd694727be06f0f8f2445b5
  • 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.

This is a significant refactor of the restoring implementation. Separated the exterior halo boundary update and the interior restoring capability. Added namelist set_boundary_flds to set exterior halo values. Added namelist restore_data, restore_flds, restore_mask, restore_width, and restore_timescale which work when restore_ice is on to computed interior restoring terms. Calls to initialize and execute the boundary and interior restoring were updated in the code.

In summary, several model defined fields can be restored (controlled by set_boundary_flds and restore_flds) in ice_restoring.F90. The restoring dataset is specified with restore_data. Exterior halo values are set with set_boundary_flds. Interior restoring is turned on with restore_ice and the restoring strength is specified with restore_mask, restore_width, and restore_timescale. The implementation provides some basic options and a framework for adding new fields, new masks, and new datasets by application users as needed.

The current restore_data options are

• defined, internally defined state at all time
• initial, uses the initial model state at all time
• restartfiles, reads in extended restart files at each timestep, filenames are somewhat hardcoded into the implementation

The fields (set_boundary_flds, restore_flds) currently supported are

• aicen, restore ice concentration in all categories
• state, restore aicen, vicen, vsnon, and trcrn
• trcrn, restore all tracers in all categories
• vicen, restore ice thickness in all categories
• vsnon, restore snow thickness in all categories
• velocity, restore uvel and vvel

The interior restoring strength is set as follows. The interior mask (restore_mask) specifies a mask that is applied to interior points based on restore_width (the number of gridcells around the exterior to restore) and restore_timescale (the restoring time in days) as follows

• all, sets all interior points with full restoring value (restore_timescale)
• none, no interior restoring is carried out
• constant, restores restore_width gridpoints on the exterior of the domain with full restoring value (restore_timescale)
• linear, restores restore_width gridpoints on the exterior of the domain with a linearly decreasing restoring value moving from the exterior (set to restore_timescale) to the interior.

The boundary and interior restoring uses the same internal arrays to hold the restoring data. The boundary updates always just set the values on the exterior halo to the restoring values. The interior restoring updates the interior fields with a forcing strength proportional to the model timestep (dt) over the restoring time (restore_timescale). It's important not to do interior restoring more than once per model timestep for each field. The halo updates are done via subroutine ice_restoring_halo while the interior restoring is done via calls to subroutine ice_restoring_interior.

A new timer was added, timer_restore, to track restoring time. Documentation was updated, and several tests were added or updated. There is a new script, configuration/scripts/tests/full2nest.sh, that provides a basis for extracting a smaller (nest) regional grid from a (full) run and generating the extended restart files to force the boundary and interior for the nest testcase. The full2nest.sh script describes how to setup that test.

All results are bit-for-bit except the boxrestore test case which was modified to produce new results.

There are still several issues to investigate

• Can the exterior halo values for grid and technical fields (like the advection metrics) be computed directly and exactly for a regional grid rather than by a "zero_gradient" type implementation? Up to now, the only time these values were needed on the halo was in cyclic cases and those halo values can be set by the cyclic bc exactly.
• Can/should we be able to mix "zero_gradient" and "linear_extrapolation" for regional bcs in a model run depending on fields (maybe velocity should be extrapolated and concentration should not?)
• We should implement restoring for aice, vice, vsno, tracers (not defined by category but then distributed across categories) and maybe subsets of tracer fields (instead of all or none).
• Will the regional grids work better with 2 or 3 ghost cells? We need to validate the model with ghost cells greater than 1 and test whether extending the halo and specifying the halo values will improve the regional solutions.
• Are the restoring calls in the best place in the code? For instance, can we implement interior restoring in dynamics subcycling (ndtd or ndte) while preserving the restoring forcing strength somehow?
• Should we get rid of either restore_ocn or restore_bgc?

Several other minor modifications were made

• bound_state was separated from ice_state into it's own module due to some circular dependencies created by adding ice_restoring_halo calls to bound_state to update the exterior halo when required. Also added an optional argument that turns off any restoring in bound_state, needed during initialization.
• deprecated restore_bgc feature
• maxval_spval in runtime_diagnostics was updated to better reflect the implementation
• the restoring calculation for restore_ocn and restore_bgc was maintained but updated to reduce flops
• added some memory allocation checks
• updated some abort_ice calls
• cleaned up indentation in ice_forcing_bgc.F90
• updated restoring calls in all drivers including unit test drivers as needed
• removed a couple of unused variables
• updated the Boreas port

[apcraig]

Here are some full/nest results. The full case is an 80x80 box with a centered gaussian ice distribution, east wind forcing, and cyclic boundary conditions. The nest case is a 12x12 box just northeast of the center of the 80x80 box, east wind forcing, and zero_gradient boundary conditions. The nest case initial condition is the restart from the first timestep of the full case, and the nest case is setting the bc (on the outside halo) for the state (aicen,vicen,vsnon,trcrn) and the velocity (uvel,vvel) at each timestep from the full restart files. There is no interior restoring in either case.

timestep 1, full / nest:

timestep 1, difference:

day 4, full / nest:

day 4, difference:

There are still several questions including why the gaussian blob of ice doesn't seem to be advecting east properly in the nest case. The velocities quickly spin up to free drift with aice=~1 (at all gridcells). The velocities and aice match well in the full and nest cases, but the thickness doesn't.

[eclare108213]

Wow, this is a lot to digest. I've made an initial pass through the changes, and for the most part this looks okay overall.

Is there documentation for the regional/nest test, yet?

Looking for that, I came across this relatively recent addition to the code, originally from NRL: https://cice-consortium-cice.readthedocs.io/en/main/developer_guide/dg_assim.html
It seems like that capability ought to be combined with this new restoring design, or at least documented in the same place. Pinging @NickSzapiro-NOAA since this DA adition appears to be used in one or more NOAA models (see #1060).

[eclare108213]

Is upwind transport included in the testing for these boundary conditions?

[apcraig]

I just tested the upwind advection and the full/nest results are MUCH BETTER! Differences in full-nest vice of 0.003 in upwind vs 0.12 (40x) in remap. I think this has to do with the haloupdates required for the upwind advection remap metrics. We need to improve the zero_gradient extrapolation currently used in the remap advection metrics (see first "todo" bullet point in PR).

[eclare108213]

I don't understand "upwind advection remap metrics". ?

There seems to be something wrong with the incremental remapping code that has little to do with regional modeling. See #1072

[apcraig]

Sorry for the confusion. I've been speculating that the haloUpdates needed to compute "dpx, dpy, mc, mc, my, tc, tx, and ty" in the remap advection are leading to problems in the remap advection. Those fields are what I'm calling the remap metrics. Those fields need to be computed on the halo. When the bcs are cyclic or closed, the metrics are computed correctly and easily. With the open bc, we are using "zero_gradient" to compute those metrics which is clearly not exact. One of my points in the PR (Can the exterior halo values for grid and technical fields (like the advection metrics) be computed directly and exactly for a regional grid rather than by a "zero_gradient" type implementation?) is whether we can compute the remap metrics in the halo directly. We never had to before, so we never tried. Maybe we need to look into that. The fact that the upwind advection behaves so much better than the remap advection for the full/nest test makes a stronger case that this is an issue specifically in the remap advection. But at this point, the cause is still speculation. My intuition says the remap advection issue on the regional grid is bigger than #1072. Again, more investigation is needed. I think we should merge this PR and then continue to tackle these other issues and requirements.

[NickSzapiro-NOAA]

I agree there can be nice overlap with restoring, boundary conditions, and initialization

"adjust_aice" adjusts the domain interior to some prescribed aice on restart.
I'm not sure how to combine that with options that update every timestep, but the underlying strategies for creating ice states could be shared.

If I'm following, some boundary condition options can extend the interior solution to the domain haloes.
Separate flavors like adjust_aice (or wrappers to add_new_ice) could also set haloes or restoring targets from external aice/hi/hs/...

[apcraig]

Thanks, lots of good feedback. Yes, we should consolidate the "adjust_aice" feature and the ice restoring. Let me try to do that in this PR. I will also remove the restore_bgc option. More soon.

[apcraig]

Had a look at the code again. First, I am fixing several of the issues noted above. I removed the restore_bgc code, I renamed crestore to frestore and made that scalar consistent through the model and I removed the "abs" part of the crestore calculation. I will push the changes to the PR soon.

Second, I had a look at the "adjust_aice" restart da implementation and compared it to the restoring. In theory, these should be similar and overlap. In practice, I'm not sure we're ready to do that yet. The "adjust_aice" is a very particular implementation that has a bunch of hard-wired "close enough" parameters and other "application specific" details. Plus, it is implemented not as a restoring, but as a initial condition reset. Gently restoring might need a slightly different approach. While I believe the "adjust_aice" could be generalized further and possibly shared with the ice restoring option, I think we should defer for the moment.

Third, the ability to restore aice and vice (concentration and thickness) is high priority. Right now, we support restoring to the "ncat" values, but not to a total value.

• The "ice_restoring_data_define" implementation provides one way to convert aice and vice (hbar) to the ncat layers, but it does so by placing all the ice in one category and then computing the appropriate enthalpy for the to be restored to state.
• "adjust_aice" modifies the initial ice state by removing ice from the first category first and moving up categories or adding ice just to the first category. In all cases, it adjusts the thickness, enthalpy, salinity, etc as it goes. It has some tricks when all ice is removed, etc. Again, this is a state reset, not restoring.
• The "SIS2" restoring method initially assigns the restoring ice to the "primary category", like "ice_restoring_ data_define" then redistributes to thinner categories. The aice and vice categories are restored then it specifically talks about checking and adjusting salinity and enthalpy after restoring. I don't think we've considered this yet.

So, again, a science question mostly, how should this really be done? We can assume that the model ice state is always self consistent with respect to salinity and enthalpy. We seem to make some effort to set the restoring state to be internally consistent wrt salinity and enthalpy (maybe). But right now, we do not explicitly adjust the state after restoring to make sure the updated state in internally consistent. Or maybe that happens later in the timestep anyway? And when we restore aice and vice, should we be restoring the salinity and enthalpy tracers as well or should we restore aice and vice then adjust the salinity and enthalpy based on the new model state? If we are restoring aicen, vicen, vsnon, and trcrn based on restart data, should all fields be restored with (or without) adjustments or should only a subset of fields be restored with (or without) adjustments?

Do we want to implement a "few" aice/vice restoring options. Maybe one like "ice_restoring_data_define" to start?

This PR is mostly about infrastructure, I think there are still alot of science questions and lots of implementation details still to work out that will require future PRs.

[NickSzapiro-NOAA]

Thanks @apcraig !

On the infrastructure, I've been wondering if it's nice for some of these choices to build from internal or external states to eventually live as options in Icepack.

To help the options be modular, combinable, and testable, say with the MOSAiC test case

[apcraig]

Moving the restoring into Icepack might not be ideal. All the spatially varying restoring information would have to be passed into Icepack with each call for each gridcell, so I think it would look like it does now with just a call into an icepack restoring subroutine just at the point the restoring term is computed. I think most of the restoring code would still exist outside Icepack in the end.

[eclare108213]

And when we restore aice and vice, should we be restoring the salinity and enthalpy tracers as well or should we restore aice and vice then adjust the salinity and enthalpy based on the new model state? If we are restoring aicen, vicen, vsnon, and trcrn based on restart data, should all fields be restored with (or without) adjustments or should only a subset of fields be restored with (or without) adjustments?

Whether a given field needs to be restored following changes to aice and vice depends on the field. Most tracers (e.g. salt, enthalpy, bgc) are "densities", and so if the tracer is not updated, then lopping off a chunk of ice or adding more ice retains the same "density" of whatever it is, and that is applied to the total area or volume after adjustment. That means the total quantity of tracer is not conserved, but in this case that's probably okay. I think we should follow the lead of other groups that are already doing restoring or data assimilation, because they've likely already dealt with finicky thermo convergence issues.

Do we want to implement a "few" aice/vice restoring options. Maybe one like "ice_restoring_data_define" to start?

Possibly. You've already made the infrastructure fairly flexible, so let's start with just what's needed by the developers in the regional modeling WG, then add other options later if needed.

This PR is mostly about infrastructure, I think there are still alot of science questions and lots of implementation details still to work out that will require future PRs.

totally agree

[TillRasmussen]

I agree that we should push the infrastructure changes when we are happy with these.

For the restoring_
We only modify enthalpy and salinity if new ice is added to a ice category or completely removed during restoring. As Elizabeth mentioned this add/remove total tracer volume but assimilation do not necessarily go hand in hand with conservation of properties. Especially in forced (by atmosphere) we dont really know where the bias/uncertainty originates from. It can be precipitation/evaporation (salinity) or heat.

[DmitryDukhovskoy]

And when we restore aice and vice, should we be restoring the salinity and enthalpy tracers as well or should we restore aice and vice then adjust the salinity and enthalpy based on the new model state? If we are restoring aicen, vicen, vsnon, and trcrn based on restart data, should all fields be restored with (or without) adjustments or should only a subset of fields be restored with (or without) adjustments?

Whether a given field needs to be restored following changes to aice and vice depends on the field. Most tracers (e.g. salt, enthalpy, bgc) are "densities", and so if the tracer is not updated, then lopping off a chunk of ice or adding more ice retains the same "density" of whatever it is, and that is applied to the total area or volume after adjustment. That means the total quantity of tracer is not conserved, but in this case that's probably okay. I think we should follow the lead of other groups that are already doing restoring or data assimilation, because they've likely already dealt with finicky thermo convergence issues.

FYI: Ice restoration (relaxation) has been recently added to SIS2.
https://github.com/NOAA-GFDL/SIS2/blob/dev/gfdl/src/SIS_sponge.F90
Regarding restoring salinity and enthalpy, in our code, S and enthalpy are updated in newly added ice or if enthalpy is too "warm" (to prevent rapid melting of restored ice) or bulk S < threshold (3 psu)
e.g.:
https://github.com/NOAA-GFDL/SIS2/blob/dev/gfdl/src/SIS_sponge.F90#L468

[eclare108213]

@DmitryDukhovskoy, thanks for the detail. It makes sense to only update the enthalpy and salinity when ice first arrives or disappears completely, unless a correction is needed to stay within the bounds needed by the model.

@apcraig, are you planning to make additional changes, or is this PR ready for final review (other than addressing the few remaining comments above)?

[apcraig]

I think this is ready for final review. Hopefully I've addressed the reviews thus far.

[eclare108213]

There are a few comments/questions waiting for a response.

[eclare108213]

We know what "interior restoring" means, but someone outside our core group might think it means "interior of the ice", which is partly true but not entirely correct. Every instance of "interior restoring" doesn't need to be changed, but at least once in each file where it occurs, please use "interior domain restoring" or something similar. E.g. I'd change lines 475 and 478 here, and maybe lines 608-613 in cice_index.rst.

[apcraig]

I see the point. Not sure I love "interior domain restoring", maybe we should just call it "restoring"? The problem is that we're transitioning from a single restoring feature that did boundary and interior "restoring" to two distinct features (almost) that are setting exterior boundary conditions and setting interior restoring (nudging?). I'm happy to go through the code and documentation and use different language everywhere as needed if we can decide on what terms we want to use for these two features. I think I've generally used "boundary conditions" and "interior restoring". I will define "interior restoring" formally in the code at line 478 above, but happy to do more.

[eclare108213]

Thanks for all the updates. This looks good to me, with an understanding that the infrastructure will likely continue to evolve as regional modeling needs are better defined.