mom nuopc cap changes required for running with ice shelves

[claireyung]

This issue will document the changes required in the mom6 nuopc cap to allow ice shelves to run in MOM6. A draft codebase (definitely needs some fixes/clean up) is here in branch https://github.com/ACCESS-NRI/MOM6/tree/ice-shelf-dev and it compiles with ACCESS-OM3 e.g. here.

Key issues were

• Preventing ocean below ice shelf talking to atmosphere
• Properly initialising ice shelves

Important context: to avoid ocean talking to atmosphere, everything the coupler sees should just be open ocean i.e. sea ice mask (no sea ice below ice shelves, sea ice mask + ice sheet mask must exactly match full ocean mask). MOM should keep the fact that the ocean is larger and has ice shelf cavities to itself, and does meltwater fluxes by itself.

Changes made

1. Mesh mask mismatch error: Ocean nuopc mesh mask not matching MOM6 ocean mask needs to be allowed. We replace the error ERROR: ESMF mesh and MOM6 domain masks are inconsistent! with a warning instead of fatal error: call MOM_error(WARNING, err_msg)
2. So_omask export: Needed to change the So_omask export (omask variable) in mom_export subroutine of mom_cap_methods.F90. This is because MOM was using the internal mask, which differed from the mask given to the mediator to initialise e.g. the atmosphere-ocean fluxes, which is a static field defined in nuopc.runconfig. We had given nuopc the sea ice mask not the ocean mask to avoid the coupler sending fluxes into the cavities, but MOM was sending data back in these files resulting in NaNs in the mediator. The change is here

      omask(i,j) =  nint(ocean_grid%mask2dT(ig,jg)) - ceiling(ocean_state%ice_shelf_CSp%ISS%hmask(ig,jg))

and also required making ocean_state_type and ice_shelf_CS contents public (i.e. removing ; private for each of them) Please see note 1

3. Correctly initialising ice shelf in MOM nuopc cap: By comparing the FMS and nuopc cap I realised the ice shelf was being initialised differently, so added the following:

add correct modules:

use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces

initialise ice shelf fluxes and forces:

  if (OS%use_ice_shelf)  then
    call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%fluxes)
    call initialize_ice_shelf_forces(OS%ice_shelf_CSp, OS%grid, OS%US, OS%forces)
  endif

Please see note 2

Note 1:

• The subroutine State_SetExport also uses the ocean grid internal mask, applying a mask again as it exports. It doesn't do anything really, and I didn't need to change it.
• Potentially other places need this mask too. E.g. So_t currently has ocean in the cavities in the mediator output.... but it doesn't seem to be a problem because So_omask is correct. Maybe also in the input.

Note 2:

• MOM6 FMS cap just uses the fluxes and forces initialisations and not initialize_ice_shelf which I left above. If I did only the fluxes and forces I got a segfault.

Suggested improvements to code

• Currently I rely on making control structures ocean_state_type and ice_shelf_CS contents public to allow the nuopc cap to access them. Ideally this shouldn't happen, may need some feeder submodules that allow the required variables to be accessed
• I'm using ISS%hmask for the ice shelf - this could have values not 0 or 1 though (e.g. 2). Using fluxes%frac_shelf_h could also work (this was a previous attempt but the variable hadn't been properly initialized), but likewise can vary and have values between 0 and 1.
• I don't have an "if using ice shelf then apply this different omask" statement, this needs to be added to mom_export.

special thanks to @dougiesquire @anton-seaice for helpful discussions/hacking during the ACCESS workshop.

[anton-seaice]

1. Mesh mask mismatch error: Ocean nuopc mesh mask not matching MOM6 ocean mask needs to be allowed. We replace the error ERROR: ESMF mesh and MOM6 domain masks are inconsistent! with a warning instead of fatal error: call MOM_error(WARNING, err_msg)

Can we change this so that the mesh mask is compared to the (MOM6 mask minus ice shelf mask) ? (e.g. are the ice shelves initialised at this point in the cap ?)

• The subroutine State_SetExport also uses the ocean grid internal mask, applying a mask again as it exports. It doesn't do anything really, and I didn't need to change it.

I think we should change this for completeness, and for future compatibility for an ice shelf model

• I'm using ISS%hmask for the ice shelf - this could have values not 0 or 1 though (e.g. 2). Using fluxes%frac_shelf_h could also work (this was a previous attempt but the variable hadn't been properly initialized), but likewise can vary and have values between 0 and 1.

I am not 100% sure, but I think a fortran Where statement would help

WHERE (hmask>0) <<do some assignment>>

[anton-seaice]

As always - ocean team can do this work, as long as you have something which works !

[claireyung]

I feel like there is some abnormal behaviour regarding ice shelf restarts with the NUOPC cap compared to how it is in FMS-SIS2.

I'm trying to restart from the 10 day run but having a few issues, wondering if anyone has ideas....

• The ice shelf restart file Shelf.res.nc got saved in output000 rather than restart000 which was odd (nuopc thing? payu thing?). I manually copied it over, but something to fix later.
• It ran but seemed to initialise from a cold start and neglect the ocean restarts - CFL was zero at intialisation, was really big in day 11 (first day after restarting), melt rates look like they did at day 1 and not like day 10. Suspect this is because I am using input_filename = 'n' in &MOM_input_nml in input.nml, this was needed to initialise the ice shelf module at day 1 (and it works fine in MOM6-SIS2 with this parameter on all the time) otherwise it crashed telling me it could not find any restarts, but maybe it adversely results in the ocean not using restarts with the nuopc cap.
• If I get rid of input_filename = 'n' then I get

FATAL from PE     0: MOM_restart: Unable to find any restart files specified by    in directory ./.

which is the reason I added it in the first place. Key code is here I think:

  new_sim = .false.
  if ((dirs%input_filename(1:1) == 'n') .and. &
      (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true.

And below is also relevant, this is right before it crashes if you don't have input_filename = 'n': https://github.com/ACCESS-NRI/MOM6/blob/2025.07/src/ice_shelf/MOM_ice_shelf.F90#L1906

  elseif (.not.new_sim) then
    ! This line calls a subroutine that reads the initial conditions from a restart file.
    call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.")
    call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, CS%restart_CSp)

Is there something wrong with the nuopc cap that is preventing the restart file from being found or am I missing something that I'm supposed to tell the model?

[dougiesquire]

I think this is because we use the NUOPC cap with restart pointer files and this (and Payu) is not set up to handle the additional restart file. I encountered a similar thing with generic tracers and my hack at the time was to add an optional additional_restart_dir field to the NUOPC runconfig. If you can direct your ice shelf restarts to that location they should get handled by Payu correctly and be available to MOM for restarting the next run.

Note, this is a hack that I've been meaning to clean up for a long time...

[claireyung]

Thanks Dougie!
Not quite sure where that location of the additional_restart_dir should be - do you have an example config? Or do you mean change the code?

[anton-seaice]

See https://github.com/ACCESS-NRI/access-om3-configs/blob/4dcf8d6833237852fd1390486273d1a82dfdaee7/nuopc.runconfig#L121C6-L121C28

although i think it requires being able to set Shelf.res.nc to be saved in a directory named RESTART ?

[dougiesquire]

Yes, it would require code mods to write to additional_restart_dir, which can be set in the nuopc.runconfig file and defaults to ./RESTART. If that's not straightforward, perhaps it's time for us (me) to come up with a better solution for handling these additional restarts.

If you just wanted to hack something for testing, you could hardcode to write the restarts to ./RESTART.

[claireyung]

Thank you!! Here is some documentation of what works...

Changes that were needed to make it work with restarts (hacky!)

Commit 1: 1424586

• Hardcode ice shelf restart to save in './RESTART/' in nuopc cap. (And the / is important!). This allows payu to send the restart to the correct folder and be available for the next run.
• (ideally needs to use additional_restart_dir instead to be slightly less hard-coded)
• Also remove fluxes, forces) optional arguments from initialize_ice_shelf; these are initialised immediately afterwards but including them here incorrectly makes the shelf_sfc_mass_flux field be initialised when it should only be used in dynamic ice shelf cases (and results in a crash when it can't find the file if you're running from a restart.)
• Also hardcodes the name of the ice shelf restart Shelf.res.nc in initialize_ice_shelf (why was this necessary - it isn't in MOM6-SIS2)
• Important: for the very first run, I need to set input_filename = 'n' in the input.nml for the ice shelf to initialise. Then for all future restart runs this line should be commented out. This mid-run editing of input files is a bit annoying, but I also need to do it for the timestep too since for the first couple of days the timestep needs to be small, and then we can increase it to make the model run faster.

Commit 2: b2e31ba

• use frac_shelf_h instead of hmask in So_omask calculation. It seems that hmask does not get initialised from a restart (at least in static ice shelf mode) and therefore it has zeroes, making the So_omask include the cavities and an immediate NaN error on restart after this run (i.e. running from restart001 it crashes, but running from restart000 seemed fine since the omask was correct in restart000?) does So_omask only get calculated once at initialisation or can it vary? Seems odd that it ran fine from restart000 without ISS%hmask defined?
• with this, we should make the ice_shelf_CS contents private again
• what happens if frac_shelf_h is between 0 and 1?

[dougiesquire]

Thanks for recording this @claireyung. Just putting some of the things we just discussed down in writing.

• (ideally needs to use additional_restart_dir instead to be slightly less hard-coded)

This should be relatively straightforward, but I'm hopeful there's a better approach. We now have two restarts that don't play nicely with the CESM rpointer approach. This is good motivation to try come up with a less hacky solution.

• Also remove fluxes, forces) optional arguments from initialize_ice_shelf; these are initialised immediately afterwards but including them here incorrectly makes the shelf_sfc_mass_flux field be initialised when it should only be used in dynamic ice shelf cases (and results in a crash when it can't find the file if you're running from a restart.)

As you noticed, we should be able to remove the double invocation of initialize_ice_shelf in MOM_ocean_model_nuopc::ocean_model_init and MOM::initialize_MOM by passing ice_shelf_CSp to initialize_MOM in mom_ocean_model_nuopc.

• Also hardcodes the name of the ice shelf restart Shelf.res.nc in initialize_ice_shelf (why was this necessary - it isn't in MOM6-SIS2)

• Important: for the very first run, I need to set input_filename = 'n' in the input.nml for the ice shelf to initialise. Then for all future restart runs this line should be commented out. This mid-run editing of input files is a bit annoying, but I also need to do it for the timestep too since for the first couple of days the timestep needs to be small, and then we can increase it to make the model run faster.

Could you try not hardcoding Shelf.res.nc in initialize_ice_shelf and instead setting input_filename = 'r' for subsequent runs. This is still not a long-term solution, but it would be good if we could avoid this hardcoding while we work out a better general approach.

• use frac_shelf_h instead of hmask in So_omask calculation. It seems that hmask does not get initialised from a restart (at least in static ice shelf mode) and therefore it has zeroes, making the So_omask include the cavities and an immediate NaN error on restart after this run (i.e. running from restart001 it crashes, but running from restart000 seemed fine since the omask was correct in restart000?) does So_omask only get calculated once at initialisation or can it vary? Seems odd that it ran fine from restart000 without ISS%hmask defined?

I think we agreed that using frac_shelf_h is probably better anyway

• what happens if frac_shelf_h is between 0 and 1?

I think we will export another field for the fractional coverage. Fluxes will need to be weighted by the fractional coverage in the mediator. But we're a wee way from needing to worry about that yet I think.

[claireyung]

Thanks Dougie! I tried with no hardcoding of the Shelf.res.nc in ice shelf initialisation, and

input_filename = 'r'

but I got

FATAL from PE     0: MOM_restart: Unable to find any restart files specified by  r  in directory ./.

[aekiss]

Does something also need to be done to handle frazil differently in cavities?

See explanation here.
So-called "frazil" is actually a heat IOU: to prevent water falling below the local freezing point, heat is borrowed from the overlying cell. The overlying cell may in turn need to borrow heat from the cell above to remain liquid, so this process is iterated upwards. If any of this IOU remains after this upward iteration reaches the surface it is paid by the sea ice, where the extraction of heat may drive production of the frazil component of the sea ice (and reject brine into the ocean model).

So what happens inside a cavity, where there is no overlying sea ice model to extract the frazil heat from?

• Can we just assume the heat is extracted from the ice shelf? If so, is it sensible or latent heat? If latent (ie freezing frazil ice onto base of ice shelf), should there be an associated reduction in freshwater melt input to the ocean?
• Should there be associated brine rejection? At depth where the "frazil" first forms?