Update halo and interior restoring#1110
Conversation
…uld not be restored
…dyrestore1 Update to #ce4dd0b
apcraig
requested review from
TillRasmussen,
anton-seaice,
dabail10,
eclare108213 and
isievers
|
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
left a comment
eclare108213
left a comment
There was a problem hiding this comment.
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).
| call ice_timer_start(timer_bound) | ||
| call ice_HaloExtrapolate(uee, distrb_info, & | ||
| ew_boundary_type, ns_boundary_type) | ||
| call ice_HaloExtrapolate(vnn, distrb_info, & |
There was a problem hiding this comment.
Is upwind transport included in the testing for these boundary conditions?
There was a problem hiding this comment.
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).
| crestore = c1 ! use data instantaneously | ||
| else | ||
| trest = real(trestore,kind=dbl_kind) * secday ! seconds | ||
| crestore = max(abs(dt/(trestore*secday)),c1) |
There was a problem hiding this comment.
Why is abs (absolute value) needed? Should this be int() instead?
| if (trestore == c0) then | ||
| crestore = c1 ! use data instantaneously | ||
| else | ||
| crestore = max(abs(dt/(trestore*secday)),c1) |
There was a problem hiding this comment.
same comment as for restore_ocn
There was a problem hiding this comment.
I recommend removing restore_bgc competely
| "``restore_data``", "defined", "set restoring data to internally defined values", "unknown" | ||
| "", "``initial``", "set restoring data to initial conditions", "" | ||
| "", "``restartfiles``", "set restoring data from restart files", "" | ||
| "``restore_flds``", "", "array of strings specifying interior boundary fields restoring", "" |
There was a problem hiding this comment.
should the word 'boundary' be removed from this descriptor?
| ! PRIVATE: | ||
|
|
||
| real (dbl_kind) :: & | ||
| crestore ! restoring value, dt/trestore |
There was a problem hiding this comment.
The actual calculation of crestore is more complicated than dt/trestore. Maybe just say "! restoring value (days)" ?
| that define which fields to set. The boundary update | ||
| is carried out with a call to *ice_restoring_halo* method in the file **ice_restoring.F90**. | ||
| The exterior boundary data is set with the ``restore_data`` namelist option and the halo | ||
| is always set to the boundary data, not restored, despite some of the naming conventions. |
There was a problem hiding this comment.
So trestore = 0? Or is it imposed via logic?
| to ensure the restoring forcing is consistent with the namelist settings. | ||
|
|
||
| The interior restoring strength is specified by the ``restore_mask``, ``restore_width``, and | ||
| ``restore_timescale`` namelist. These specified where and how strong to restore the interior |
There was a problem hiding this comment.
change 'specified' to 'specify'
| this_block ! block info for current block | ||
|
|
||
| real (dbl_kind) :: & | ||
| fval, & ! fval computed from restore_timescale |
There was a problem hiding this comment.
"! restoring parameter computed from restore_timescale" (?). Is this the same thing as crestore elsewhere? If so, I'd prefer to use one or the other rather than both.
| indxj(:) = 0 | ||
| !----------------------------------------------------------------- | ||
| ! initial area and thickness in ice-occupied restoring cells | ||
| !----------------------------------------------------------------- |
There was a problem hiding this comment.
Add a comment here that the hardwired values in this subroutine are just placeholders and the user is expected to provide the data definitions?
| if (restore_timescale == c0) then | ||
| fval = c1 | ||
| else | ||
| fval = min(abs(dt/(restore_timescale*secday)),c1) |
There was a problem hiding this comment.
similar questions about this calculation
|
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. If I'm following, some boundary condition options can extend the interior solution to the domain haloes. |
|
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. |
|
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.
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. |
|
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 |
Deprecate restore_bgc Cleanup indentation in ice_forcing_bgc.F90 Update documentation
|
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. |
3 participants
PR checklist
Update halo and interior restoring
isievers, apcraig
results are bit-for-bit except boxrestore tests have changed answers.
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
The fields (set_boundary_flds, restore_flds) currently supported are
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
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
Several other minor modifications were made