C grid gx3 tests: Kimmritz et al. 2016 C1 configuration

[JFLemieux73]

No description provided.

[JFLemieux73]

I have bad news...I just ran gx3 for one year with B grid EVP, B grid VP and C grid EVP.

I used this commit: bffe4e8

ndte=1200
visc_method = 'avg_zeta'
elasticDamp = 0.12d0
deltaminEVP = 2e-9
deltaminVP = 2e-9
capping = 0.

Here is the thickness diff (Bvp-Bevp):

They are very close (range: -0.1 to 0.1). Here is now Cevp-Bevp. Notice that the range is now -1 to 1 m.

That's not cool!

[JFLemieux73]

I don't think it is related to elastic waves. Is it possible the issue is related to advection or ridging (maybe the deformation subroutine has a bug?).

I am going to use boxwallblockp5 for debugging.

[JFLemieux73]

Here are a few things I would like to test/investigate:
-thermo on or off
-test new deformation_T2 subroutine
-test kstrength = 0 or 1
-test remap and upwind
-issue with metric terms?

[dabail10]

I think that adding the viscosities to the history might not be worth it. The B-grid values are just local to the subroutine stress. I think with the gbox12 or even gbox80, you should just print these. I can try a test where I print them for a B and C grid. I did this awhile back and I remember that the issue was the values of these where you had a a jump from ice to no ice. This was improved when using the avg_strength method. Keep in mind for the B-grid, you are never averaging the viscous coefficients and you only use the strength at the T-points.

[JFLemieux73]

Thanks Dave. For the moment I don't need to look at viscosities. After my tests with boxwallp5 I don,t think there is an issue with viscosites. We will see...

[JFLemieux73]

I ran only one month of gx3 and I see the same problem. I look at monthly means below (January 2005). The max diff hi (B vp) - hi (B evp) is only around 1 cm. Here is the difference between hi (C evp) and hi (B evp):

[eclare108213]

I suggest trying a run with thermo and advection and ridging all off, to compare divu. If that's okay, then turn ridging on and compare again.

[JFLemieux73]

There are large differences. In the first test I just change deform_option to 2 (new method in my deformT branch). Good idea Elizabeth. The new deform method does not solve the problem. The difference between option 1 and 2 is around 1 cm but the diff between C evp (option 2) and B evp is similar to the fig above.

[JFLemieux73]

Hello,
Questions for @eclare108213 @apcraig @dabail10 :
-has anyone looked at the new grid lengths (e.g. dxE) for the gx3 test case. Do they make sense?
-Am I right that the B grid starts from a complete restart but that the C grid does not? I don't think it would explain the big differences in thickness but it could help me for debugging.

[JFLemieux73]

Follow up question. Let's say I want to have the same initial conditions. What would be the best way to set the velocities and stresses to zero at the beginning? In dyn _prep2 but only done for the first time level?

[eclare108213]

Good morning @JFLemieux73. Can you post a figure that shows dxE etc that look odd to you? I haven't looked at those yet.

You can set the B- and C -rid test cases with gx3 to have the same initial conditions using ice_in -- you shouldn't need to go into the code itself to do that. There is probably also a set_nml option to do that -- I'll have to look to see what they are. gx3 does usually start from a B-grid restart file, but you can use no-ice or default (slab) conditions instead for both grids.

And those latter options start from rest.

[JFLemieux73]

Thanks Elizabeth. I have not looked at dxE or other grid lengths. As our tests are good on uniform grids I am wondering if the weird results on the gx3 grid can be due to these terms.

[dabail10]

Both the B and C grids start from the same gx3 restart. The Cgrid velocities and stresses are initialized to zero though. This would account from some initial differences in the first couple of days.

Dave

[JFLemieux73]

Thanks Dave. I have another question. As Elizabeth suggested, I tried to run without advection, ridging and thermo. I set
ktherm=kridge=ktransport = -1

The problem is that uvelE=vvelN=0 during all the simulation....basically nothing happens. I guess tauair is not calculated. Any idea what I could do?

[apcraig]

I've been following progress. A few thoughts. Do all the problems in the gx3 grid start next to the land boundary? It looks like the box cases are not showing the same issues. The box cases (except the kmtislands) have very limited land geometry. For instance, the land is all constant in i or j or are a corner. The gx3 will have the opposite of a corner (a "peninsula") which is not in the box cases.

We are masking the computation of uvelN, vvelE, uvel, vvel after we grid_average_X2Y. Is this working properly? We should check in the gx3 case whether the B and C uvel and vvel fields have values of zero at exactly the same locations. I believe they should.

We impose the BC in the C grid solver via an implementation that imposes a mirror image inside the boundary. Is this not working correctly in non-uniform grids? We should check. Also, do we need to leverage that information somewhere else where we're not. I would say we should leverage it in the uvel, vvel calculation, but I think the masking is doing it via a different method.

Finally, the grid parameters can be written out to the history file with f_dxe, f_dye, f_dxn, f_dyn, f_nmask, f_emask, etc. We should probably add that to the default namelist file.

[dabail10]

That makes sense about the wind stresses and the fluxes. The subroutine icepack_atm_boundary which is called within icepack_step_therm1. This is called from step_therm1 in CICE. I believe step_therm1 is not called when ktherm = -1. So, we need to bring the call to icepack_atm_boundary up or move the if condition on ktherm into icepack.

[eclare108213]

I don't recommend refactoring the logic in the code for this test.
It's easier to just set a non-physical wind stress as a default (e.g. constant, u=v diagonal), since this isn't intended to be a realistic run anyhow - you're just testing the velocity and stress calculations with complicated boundaries. And that's what the kmtislands land mask was intended to do...

[JFLemieux73]

I added some more diagnostics for debugging. This is the Arctic volume for the first month (ktherm=2, ktransport=kridge=1):

And the Arctic RMS ice speed:

So it looks like it is a bit faster with the C-grid. I guess this leads to more deformations and therefore more growth.

[JFLemieux73]

I have replaced the initial condition to ice_ic = 'default'. The same things happen: thicker and slightly faster for the C grid. Here is the volume for one month:

[JFLemieux73]

Here is the slab experiment repeated but I used kstrength=0 (Hibler) with P*=0.0 (no rheology). Here is the volume for both B and C grids:

[eclare108213]

I remember reading a paper a while back, maybe by Martin Losch??? which characterized differences in (E)VP discretization in terms of an apparent ice strength. Do you remember that? Maybe you were an author, @JFLemieux73! I don't remember whether the comparison was B vs C for EVP, or if it was VP vs EVP.... If we rule out other potential causes like varying grid lengths, then perhaps that is the reason. But we need to rule out other things first.

[JFLemieux73]

Do you mean:
On solving the momentum equations of dynamic sea ice models with implicit solvers and the elastic–viscous–plastic technique?

This is Losch and Danilov 2012.

[eclare108213]

Sounds right!

[dabail10]

The gx3 forcing is now JRA! My tests were five year runs using JRA forcing from 1995 to 2000. I could certainly do gx1 runs for better confidence.

[eclare108213]

The gx3 forcing is now JRA!

Right - I'm still trying to keep up with all of you -

[dabail10]

So, advection = 'upwind' versus 'remap' for the Bgrid passes the test in the NH, but fails in the SH. :)

[JFLemieux73]

Thanks Dave. Interesting. Does it mean our QC test should only be used to test non BFB changes to the code but not important modifications to the physics or numerics (e.g. advection scheme, grid)?

One thing to note also is that as opposed to remap, C grid with upwind has less volume than B grid with upwind.

[dabail10]

I do believe the QC test is for science changes as well. It is interesting that the "two-stage" test passes, but not the quadratic skill test. This was originally designed by Andrew and Matt Turner. I will run some gx1 tests just to make sure we are not being caught by statistical significance with the gx3.

[apcraig]

Before we draw a lot of conclusions about QC in these tests, we should repeat them with gx1. QC was developed with gx1 resolution in mind, and I think that's important. It could be the results are the same, but we should check/confirm.

[JFLemieux73]

One more thing to consider...we are seeing changes in volume of ~1e12 m3. If you look at Koldunov et al. 2019 (Fast EVP solutions in a ...) Figure 4 b and d you can see that the impact of ndte is significant. The delta volume are at least twice what we are seeing here...although at gx3 resolution the impact is probably a lot less. But what I am trying to say is that subtleties in deformations (like changing the grid...) can have a strong impact on the volume.

[dabail10]

The gx1 C grid versus B grid passes the two-state test, but fails the quadratic skill test in both hemisphere. The C grid solution is generally thicker close to land, but then a thinner away from the coast. I think the cice.t-test.py script has these backward again, but I have provided the ncdiff between the two runs at 2010-01-01. The scale is -6 to 6m. I have stayed with the default values for ndte and so on for these tests. I agree that the solutions would change with ndte=1200.

[JFLemieux73]

This thicker ice is puzzling...Here are a few more tests.
ndte=1200,
kstrength=0,
visc_method = 'avg_zeta'
elasticDamp = 0.12d0
deltaminEVP = 2e-9
capping = 0.

I think we can see the change in volume pretty quickly. Here are results for the first 60 days.

[JFLemieux73]

In the next experiments, everything is the same except that P*=0.0

[JFLemieux73]

Basically we see the same thing with rheology turned off. In a sense this is good news; it looks like the problem is not related to the new rheology discretization (including the metric terms). I also tested with coriolis='zero'. Qualitatively speaking we obtain the same results (not shown).

[JFLemieux73]

So what explains the thicker ice with the C grid?
-the different staggering
-the interpolation for the advection
-wind stress (see other results above)

Other ideas?

[JFLemieux73]

I also ran long experiments. I wanted to run 5 years but for some reasons it stopped after 4 years (on 20081231). I got the error message:
(abort_ice) error = (ice_read_nc_xy) ERROR: not enough records airtmp

@dabail10 @apcraig...any idea what's going on?

[JFLemieux73]

As I don't have 5 years here is the thickness field for the B grid after 3 years (January 2008):

And for the C grid:

[JFLemieux73]

They are not so different. They are a lot more similar than in Dave's results. @dabail10 can you check is something is wrong in your runs or in your graphs? I find that the ice is very thin in your B grid results shown above.

Or maybe the fact that I get more similar results is due to the use of other settings (i.e., avg_zeta, ndte, ...).

[dabail10]

This is a problem with the forcing and leap years. If you set leap_years to .false. there is enough for the gx3. Do you have a diff here? Keep in mind that I did not change deltamin or ndte in my runs. I wanted to see "out-of-the-box" behaviour.

Dave

[JFLemieux73]

Thanks Dave. I will restart the runs. Here is C minus B with a range of -0.5 to 0.5 m:

[apcraig]

I think we should close this. Anyone disagree?

[eclare108213]

One question that came up in this issue is whether how the variable grid spacing is handled could be a source of the differences. There is now a box test for that, which I recommend running as a sanity check. Otherwise, if @JFLemieux73 and @dabail10 are satisfied that the gx3 tests are sufficiently close and/or understood, then we can close this issue.

[JFLemieux73]

I will try to find some time this week to run the new box test to see if the variable grid spacing is handled properly in the C grid code. I am still puzzled by the thicker ice with the C grid in gx3, gx1 runs. Another way to verify if this is related to the variable grid would be to set up a pan-Arctic polar stereographic test (dx=constant) with the 'same' forcing as gx3,gx1. It would be interesting to see it we would still get thicker ice with the C grid. Would it be difficult to set up such a test?

[eclare108213]

I think it would be difficult to make the forcing sufficiently similar. What you could do, though, is use idealized forcing that is the same for the dx=constant and variable grids. But that's essentially what the box test is! You might could make the test more "challenging" for the numerics by using more complicated forcing, or more complicated land masks. The box tests have lots of options now.

[JFLemieux73]

Ok I will test the variable grid spacing. Could someone help me with the setup options?

[eclare108213]

There's an overview of the various namelist options here
https://cice-consortium-cice.readthedocs.io/en/main/user_guide/ug_implementation.html#rectangular-grids
with hopefully enough detail to set up a variable-resolution case. Is it enough?

[JFLemieux73]

Ok thanks for the info.

[JFLemieux73]

See issue #786 for results about the variable grid spacing.

[eclare108213]

work on this issue is continuing in #791

[apcraig]

Closing, active discussions of C-grid testing in #792. Feel free to reopen is needed.