1 deg CICE grid very slightly misaligned with MOM

[aekiss]

The 1deg CICE grid values are very slightly different from MOM in both latitude and longitude for both T and U points - see this notebook.

Perhaps there was loss of precision due to conversion to single precision and back to double?

Possibly related: #271

[aekiss]

would be good to have grid sanity checks to catch this sort of thing #129

[anton-seaice]

This should be fixed by writing a tool to generate the CICE grid from the MOM-supergrid, which would also give us known provenance of the CICE grid files. (The tool could go in om3-scripts)

OM3 has some sanity checks (with tolerance set by eps_imesh) but there could be more.

[anton-seaice]

@aekiss - I don't have permissions but I think this should go on the CMIP7 project board.

[anton-seaice]

The CICE grid.nc file needs these variables defined (per ice_grid) :

! (1) ULAT  (radians)    \\
! (2) ULON  (radians)    \\
! (3) HTN   (cm)         \\
! (4) HTE   (cm)         \\
! (7) ANGLE (radians)

I think we can derive all these from the mom supergrid.

Angle is defined as 'ANGLE = angle between local x direction and true east' which is consistent with MOM definition except in radians from -pi to pi. Its not clear which corner of the grid cell this is measured, although it is presumably at the U points (for consistency, and because CICE calculates an ANGLET quantity).

EDIT: More fields are needed than this for OM2, but they can all be derived from the mom supergrid, the right place to look is here in ice_grid

[anton-seaice]

Perhaps there was loss of precision due to conversion to single precision and back to double?

Similarly, it could have been a degrees to radian conversion using a constant with insufficient precision

[aekiss]

Plots of the grid absolute and relative errors are shown in the first part of https://github.com/aekiss/notebooks-MartinDix/blob/d6d147a/grid025_check.ipynb.

The relative error is very nearly constant (7.746366e-9) for all grids, suggesting a differing value of pi in the deg-rad conversion.

For that matter, could it just be because our numpy check uses a slightly different deg-rad conversion factor from cice?

[anton-seaice]

For that matter, it could just be because our numpy check uses a slightly different deg-rad conversion factor from cice.

This looks probably ok:

https://github.com/CICE-Consortium/Icepack/blob/f6ff8f7c4d4cb6feabe3651b13204cf43fc948e3/columnphysics/icepack_parameters.F90#L70

pi = 3.14159265358979323846_dbl_kind

>>> np.pi
3.141592653589793

[aidanheerdegen]

So there is a problem, or not?

[anton-seaice]

So there is a problem, or not?

We should fix this, but the discrepancy is <5e-6 so at most the differences in the results would be rounding error ( I think)

We should definitely fix the 0.25 degree grid.

[anton-seaice]

I deep dived on tarea too, to see what happens and to try and understand how we should be making the cice c-grid. I calculated it from the mom supergrid like this:

tarea = mom_ds.area.coarsen(nx=2, ny=2).sum()

(for sanity, I checked this is the same as mom_ds.area[::2,::2]+mom_ds.area[1::2,::2]+mom_ds.area[1::2,1::2]+mom_ds.area[::2,1::2])

I compared the calculated tarea from the cice file, and the differences south of 65N seem as expected (i.e. in the 14th significant figure).

However, in the tripole region they do become significant (in m^2):

In the cells around the tripole, the difference is >1%. Have I calculated the tarea wrong from the mom grid? Or is the cice file wrong? (The sum of the anomalies is 0)

p.s. FYI @ezhilsabareesh8

p.p.s: I guess it could different areas could be correct but calculated in different ways ... different assumptions about how round the earth is etc, but surely having them internally consistent would be good.

[ofa001]

I have mentioned to @MartinDix that we did find issues a 1 deg in strong current in Beaufort Sea was the stongest case, but corner points also were suspect.

[aekiss]

FYI the 1deg grid is briefly documented in Bi and Marsland 2010 and Bi et al 2013

[anton-seaice]

This issue was fixed here:

COSIMA/esmgrids@1b6837d

With some testing added too. We should be getting accuracy between mom & cice grids of ~1e-13.

This was added to the OM2 configs here through many PRs linked here: ACCESS-NRI/access-om2-configs#92

[dougiesquire]

I compared the calculated tarea from the cice file, and the differences south of 65N seem as expected (i.e. in the 14th significant figure).

However, in the tripole region they do become significant (in m^2)

I've been looking into inconsistencies between the ESM1.5 and OM2 grids and I think this difference in areas was actually introduced when the ocean grid was updated to the FMS mosaic format - see here for details.

[anton-seaice]

I've been looking into inconsistencies between the ESM1.5 and OM2 grids and I think this difference in areas was actually introduced when the ocean grid was updated to the FMS mosaic format - see here for details.

Is the OM2 ocean supergrid (in the FMS mosaic format?) incorrect then?

[dougiesquire]

Is the OM2 ocean supergrid (in the FMS mosaic format?) incorrect then?

That's what I'm wondering. I need look deeper into how those areas are calculated

[ofa001]

Hi @dougiesquire (@anton-seaice), its possible that the precision was increased with access-om2, from the much old access-esm1.5 (originally access-1-0/3 grid) , as you say, as I said when I checked the other day the grid points were nominally the same but the precision as different so yes I can imagine deriver quantities like area of grid boxes will also be impacted.

[aidanheerdegen]

but the precision as different so yes I can imagine deriver quantities like area of grid boxes will also be impacted.

Why only in the tripole? And at the North Pole of a rectilinear grid, so not actually a special location for the tripolar grid?

[dougiesquire]

Yeah, I don't think this particular issue is a precision thing. The difference is introduced in the conversion from the old FMS grid format to the mosaic format using FRE-NCtools transfer_to_mosaic_grid. It seems that the method that's used to calculate polygon cell areas in FRE-NCtools has an edge-case for polygons (cells) at/near a pole. There's ongoing work in FRE-NCtools to handle this edge case. My best guess at the moment is that our issue is related to that. (See the issue I linked above for more details).

I'll probably just open an issue on FRE-NCtools soon. Hopefully someone there can clarify.

[aidanheerdegen]

So that means the old grid has the "more correct" areas?

[dougiesquire]

So that means the old grid has the "more correct" areas?

They're both incorrect, but by slightly different amounts 🙃. If anything, the mosaic format grid (OM2) has the "more correct" areas. The following figure shows the areas along an arc of cells passing through and near the north pole. The bump at the pole shouldn't be there.

There's an option in the FRE-NCtools polygon area algorithm to calculate the area of polygons near the pole using a copy of the polygon that has been rotated away from the pole. Unfortunately this isn't exposed to transfer_to_mosaic_grid, but I added it and that gives the following areas through the pole:

These figures are taken from this notebook which contains more info.

[ofa001]

HI @dougiesquire @aidanheerdegen, these grid issues at the pole have been really interesting, one wonders what the implications have been for some of the odd issues we saw at the tripole seem in some of our runs which we never really.
liked. @aekiss looked at them a lot at 0.1 degree but we had them all the way back at 1 deg in ACCESS1-0/3 as well.

[aekiss]

There's this probably-unrelated issue which we've never figured out: #86

[ofa001]

Hi @aekiss that's was the link I want you to add in, its probably unrelated but its good to line them up at some point.