Some updates to the ROMS sea-ice code

[talebu]

Hi everyone

After working with the ROMS sea ice model for some months now, I have made some changes to the official ROMS code in order for the sea ice model to work correctly. All my changes are available at https://github.com/taleMET/roms-seaice/tree/develop. A summary of the changes can be found below:

1. Allocation and initialization of ice model boundary states for all ice state variables are moved to ROMS/Modules/mod_boundary.F. Read of boundary data for ice state variables is added to ROMS/Nonlinear/get_data.F and ROMS/Nonlinear/set_data.F, and the boundary state variables are added to the varinfo file.
2. The ice strength Pstar is now read from the ice.in file. Without this, the ice was completely slushy and would move almost only dependent on the ocean currents.
3. For the closed lateral boundary condition for internal ice temperature (T_i) the temperature is set to be the same as in the point next to it, instead of using the computed enthalpy and the ice thickness. This leads to a more stable boundary, where without this, some blobs of warmer ice would come in from the boundaries.
4. In the ice bulk flux scheme, the transfer coefficients for sensible heat (Ch) and latent heat (Ce) are set to constants (1.75E-3_r8). This is consistent with the code developed by Kate earlier. In the previous implementation (before my changes), these constants are set to Ch=Wstar(i)*Tstar(i)/(-Wspeed*delT(i)+0.0098_r8*blk_ZT(ng)) and Ce=Wstar(i)*Qstar(i)/(-Wspeed*delQ(i)+eps). However, this would lead to very large numbers as the denominator could possibly be 0 (or a very small number). I also think these equations might be wrong, and that the intent was something else when this was implemented.

Let me know if there are any questions or comments about this.

Tale

[dylanrs617]

Thanks for posting this, I will test the new code. I was using the older version of your branch that had the BRY fixes in it. One recent bug I've discovered is that the ROMS ice model is not bit for bit and will seg fault or impose grid scale noise depending on the tiling. I've attached some of the setup information (header, ocean.in, ice.in), it is an idealized model of the Alaskan Beaufort shelf break that undergoes baroclinic instability. My branch can be found at https://github.com/dylanschlichting/rom ... seaice-dvd

The noise does not appear if I choose any of the following tiling: 2x2, 4x4,16x5, but 16x20 produces the plot below of the vertical velocity within the first hour and 16x16 produces a seg fault. I am positive this is a parallel bug, I have several other grids that run stably for over 4 months. But for a 400x1000 mesh, we want to throw more than 80 cores at the problem. See attached for a comparison of 16x5 versus 16x20 tiling. My target tiling is 16x32 because it is easy to get up to 512 cores on our HPC, NERSC. I have not tried to replicate this with the lake ice test case, but this seems like an important problem that needs to be looked into.

w_plot_16x5_tiling.png (44.29 KiB) Viewed 1374 times

w_noise_16x20_tiling.png (42.09 KiB) Viewed 1374 times

Last, see attached for a snapshot of a 500 m resolution grid on a smaller domain, 200x300 km. Periodic E-W boundaries, wall at S and Rad at N for the 3d variables. The model starts to radiate these waves away from the jet, which is initialized at ~Y = 110 km. If you're wondering what those waves are, I haven't figured that out yet but the radiation BRY conditions seem to be working, I don't see any reflection. Could be forcing related, could be inertia gravity waves that are partially resolved at 500 m resolution. This was to provide more evidence that the model works, but it will not run stably if the "incorrect" tiling is used.

[talebu]

Hi

Thanks for testing my code!

I experienced some problems earlier when using the ice model. Then the model would blow up without any particular reason. This happened both to the lake ice test and to my model. The fix for this was to compile with FFLAGS := -mcmodel=large. I am not sure if this would fix your problem, but it is worth a try.

Also, do you get the same problem if you compile without ice activated?

This might also be caused by the instabilities from wrong implementation of the bulk flux scheme. Will you try to use my latest code, and see if this fixes the problem?

[dylanrs617]

Hi Tale,

I solved the bit for bit problem. This was hard to track down, but I turned the ice model off and still had the issue. Needed to provide the forcing for bulk fluxes in a .nc file instead of all analytical. I couldn't track down exactly what analytical settings were causing the problem, but it runs now regardless of tiling. I can run on 512 cores too with your updated code and compiling option as you suggested. I will let it run overnight and post an update tomorrow.

Thanks,
Dylan

[dylanrs617]

It ran stably overnight. The waves are gone and the ice model supresses submesoscale activity until the eddies are hardly visibly in the vorticity field. I can't attach an animation, so here are a few frames. I'm happy with this, it looks reasonable. I'm not sure whether it was my overhaul of bulk fluxes or Pstar being read in from the ice.in file, but this looks great. The setup details can be found in this repo, which will be updated frequently over the next few months: https://github.com/dylanschlichting/beaufort_jet.

Thanks for your work on the code!
- Dylan

[talebu]

This looks very cool! Thanks for testing my code. I will let you know if I find any other bugs, and if I manage to fix them.

[arango]

Thank you for the updates to the sea ice model. I am currently busy with other projects and cannot review these updates. I will take a look when I get a chance. I rewrote that model in a very compact form so it can be easily expanded. I am glad that others have found this simple model helpful.