Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

CICE grid generation from MOM supergrid #6

Closed
wants to merge 13 commits into from

Conversation

ezhilsabareesh8
Copy link
Contributor

@ezhilsabareesh8 ezhilsabareesh8 commented Feb 26, 2024

This pull request addresses issues COSIMA/access-om3#101 and COSIMA/access-om2#280 by adding a Python script to generate the CICE6 grid from the MOM6 supergrid.

The script facilitates the creation of a CICE grid, resolving a previously encountered longitude mismatch issue by correctly positioning the U points in the NE corner.


# Constants
ULAT = np.deg2rad(in_superGridFile['y'][2::2, 1::2])
ULON = np.deg2rad(in_superGridFile['x'][1::2, 2::2])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't this be

ULAT = np.deg2rad(in_superGridFile['y'][2::2, 2::2])
ULON = np.deg2rad(in_superGridFile['x'][2::2, 2::2])

?

it might be clearer to write:
np.deg2rad(in_superGridFile.y.sel(x=slice(start=2,step=2), y=slice(start=2, step=2)))

?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, but with isel, not sel

@anton-seaice
Copy link
Contributor

I'll mostly defer to Andrew review but we'll want some versioning information added to the output of both the source mom file and the script version used

@aekiss
Copy link
Contributor

aekiss commented Feb 27, 2024

Yes, versioning info would be good too, e.g. https://github.com/COSIMA/om3-scripts/blob/main/mesh_generation/generate_mesh.py#L397-L413

@aekiss
Copy link
Contributor

aekiss commented Feb 27, 2024

It would be a good idea to compare the new grid.nc with the old one to check that it differs only in the ways we want it to.

There might be something in here that could be useful:
https://github.com/aekiss/notebooks/blob/master/check_updates.ipynb

@aekiss
Copy link
Contributor

aekiss commented Feb 27, 2024

This is missing the following from /g/data/ik11/inputs/access-om2/input_20201102/cice_025deg/grid.nc - are they needed in CICE? I noticed there are even more things in /g/data/ik11/inputs/access-om2/input_20201102/cice_1deg/grid.nc but they apparently aren't needed, so maybe the following are also unnecessary too...?

	double angleT(ny, nx) ;
		angleT:units = "radians" ;
		angleT:title = "Rotation angle of T cells." ;
	double tarea(ny, nx) ;
		tarea:units = "m^2" ;
		tarea:title = "Area of T cells." ;
	double uarea(ny, nx) ;
		uarea:units = "m^2" ;
		uarea:title = "Area of U cells." ;

@anton-seaice
Copy link
Contributor

Those three quantities are computer internally in CICE:

https://github.com/CICE-Consortium/CICE/blob/64177e3e5658711f5f19030ddfcc6927f1dbda5d/cicecore/cicedyn/infrastructure/ice_grid.F90#L952

We might want to think a little bit about whether computing the areas in a projection aware way improves accuracy. I guess the first check would be whether the old uarea = HTN* HTE

@aekiss
Copy link
Contributor

aekiss commented Feb 27, 2024

I wonder whether it would be better to make CICE read the MOM supergrid directly?

@anton-seaice
Copy link
Contributor

It looks like code was added to cice5 to read tarea, tlat and tlon from the cice grid file rather than calculate them:

https://github.com/COSIMA/cice5/blob/edcfa6f9c76ed05b63196ce4b5355fa5a8f4fe3a/source/ice_grid.F90#L662

#ifdef AusCOM
! we also let cice read in the following fields to avoid possible mismatch
! between cice and mom4.
! ( 8) TLAT
! ( 9) TLON
! (10) TAREA
! (11) UAREA
! (12) ANGLET
! the following fields are also available in grid.nc but not read in 'cos they
! are not critical 
!  LONT_BONDS
!  LATT_BONDS
!  LONU_BONDS
!  LATU_BONDS
#endif

@anton-seaice
Copy link
Contributor

anton-seaice commented Feb 27, 2024

There is something similar in the cmeps driver for CICE:
https://github.com/CICE-Consortium/CICE/blob/64177e3e5658711f5f19030ddfcc6927f1dbda5d/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90#L731

 ! Initialize the cice mesh and the cice mask
       if (trim(grid_format) == 'meshnc') then
          ! In this case cap code determines the mask file
          call ice_mesh_setmask_from_maskfile(ice_maskfile, ice_mesh, rc=rc)
          if (ChkErr(rc,__LINE__,u_FILE_u)) return
          call ice_mesh_init_tlon_tlat_area_hm()
          if (ChkErr(rc,__LINE__,u_FILE_u)) return
       else
          ! In this case init_grid2 will initialize tlon, tlat, area and hm
          call init_grid2()
          call ice_mesh_check(gcomp,ice_mesh, rc=rc)
          if (ChkErr(rc,__LINE__,u_FILE_u)) return
       end if

Which looks like it will get the area from the mesh file.

@micaeljtoliveira
Copy link

@ezhilsabareesh8 Thanks for the changes! I think generate_cice_grid.py. was a better name, but the most important is that it's consistent.

@anton-seaice
Copy link
Contributor

@ezhilsabareesh8 - Can you test the new grid with OM2?

I think you will need to add:

! ( 8) TLAT
! ( 9) TLON
! (10) TAREA
! (11) UAREA
! (12) ANGLET

for it to work in OM2.

Should we have two versions of this grid? One for OM2 and one for OM3 (with less variables?)

I wonder whether it would be better to make CICE read the MOM supergrid directly?

Shall we split this into a different issue? Current OM3 implementation has CICE calculating area internally, which is problematic.

@ezhilsabareesh8
Copy link
Contributor Author

@anton-seaice I can test OM2 with the new grid.nc. It seems we don't need two versions of grid.nc, as the additional data will likely be ignored by CICE.

The newly generated grid.nc contains the following data:

- `ulat`: Latitude of U points (radians)
- `ulon`: Longitude of U points (radians)
- `tlat`: Latitude of T points (radians)
- `tlon`: Longitude of T points (radians)
- `htn`: Width of T cells on the N side (cm)
- `hte`: Width of T cells on the E side (cm)
- `angle`: Rotation angle of U cells (radians) 
'''

Additionally, we may need to include `TAREA` and `UREA` for OM2.

@anton-seaice
Copy link
Contributor

Oh yes - I was just thinking of two versions because its confusing if we have a grid.nc file loaded by OM3 that is then only using some of the variables.

@ezhilsabareesh8
Copy link
Contributor Author

@aekiss, I've incorporated the TAREA and UAREA variables in the latest commit (09be9fb). However, I've noticed a discrepancy between the output of the newly generated grid.nc and the old OM2 grid.nc for a 1degree grid. I'm uncertain whether the issue lies with the CICE grid or if the wrapping is not implemented correctly in the new grid. Could you please have a look at the following output?

(grid_ds['uarea'] - grid_old_ds['uarea']).plot()

Screenshot 2024-03-01 at 10 58 30 am

@anton-seaice
Copy link
Contributor

anton-seaice commented Mar 1, 2024

It looks like we need to something different at the fold join. (Calculated from mom supergrid on top, old one below)

Screenshot 2024-03-01 at 2 22 47 pm

@aekiss
Copy link
Contributor

aekiss commented Mar 4, 2024

re. the differences at the fold - not sure if this helps, but FYI we are using grid_type = "tripole" rather than tripoleT, so the northern poles and the fold between them all lie on U points rather than T points (see here). This is to match the alignment in MOM5 - see Elements of MOM sec 9.3 and fig 9.9 (copied below)
Screenshot 2024-03-04 at 2 48 43 pm

@aekiss
Copy link
Contributor

aekiss commented Mar 4, 2024

Sec 9.4.3 of Elements of MOM describes the special treatment of grid cell areas long the fold.

@aekiss
Copy link
Contributor

aekiss commented Mar 4, 2024

Re. the area mismatch, I guess the original 1deg areas in the supergrid were calculated using update_cell_area. Is the new code calculating areas the same way?

@navidcy
Copy link

navidcy commented Mar 4, 2024

Sorry to intrude in the discussion, somehow I get notifications for this.

The area of a quadrilateral on the sphere is uniquely defined given the lon-lat coords of its four vertices. It doesn’t need to be approximated. So, regardless how it was done in the past, it’s a good idea to drop any unnecessary approximations. The sum of all grid areas should be precisely 4πR^2 (up to machine precision).

@navidcy
Copy link

navidcy commented Mar 4, 2024

@navidcy
Copy link

navidcy commented Mar 4, 2024

There a way to avoid needing to ensure you give the vertices in counter-clockwise order. I can put more details when I am in front of a computer.

@anton-seaice
Copy link
Contributor

Sorry to intrude in the discussion, somehow I get notifications for this.

The area of a quadrilateral on the sphere is uniquely defined given the lon-lat coords of its four vertices. It doesn’t need to be approximated. So, regardless how it was done in the past, it’s a good idea to drop any unnecessary approximations. The sum of all grid areas should be precisely 4πR^2 (up to machine precision).

Thanks. We are trying to calculate the area of each cice grid cell by using the mom supergrid file. The areas in the mom supergrid seem to be calculated consistently with the way you are saying (https://github.com/ESMG/gridtools/blob/de0a18c1ce0807748aa70023300dfc415277bd4c/gridtools/spherical.py#L36).

This is a good thing for us to try though, if the areas calculated based on the lat/lons only are the same as the ones derived from the mom supergrid, that would be a good validation of the method.

@anton-seaice
Copy link
Contributor

Re. the area mismatch, I guess the original 1deg areas in the supergrid were calculated using update_cell_area. Is the new code calculating areas the same way?

All we are doing is summing the area of the four relevant cells from the supergrid to make the cice grid, I can't see anything inconsistent with the update_cell_area in that?

This is similar to what Navid is saving I think - we could use this method to calculate area based on lat/lon only, and not from the supergrid.

@anton-seaice
Copy link
Contributor

anton-seaice commented Mar 4, 2024

The 1deg grid has meridional grid refinement near the equator, and coarsening elsewhere

download-1

(figure from here)

This may play a role in these differences at mid-latitudes.

Does this mean the coordinates move away from regular 1 degree spacing?

p.s. Its also odd that the patterning below 65N is there in the uarea calc and not the tarea calc.

@anton-seaice
Copy link
Contributor

The patterning found above in t-area seems to only be there in the 1-degree grid, so it might be an artifact of the Meridional refinement and the grids getting out of sync at some point. Its not there in the 0.25 degree grid, so maybe we can ignore it.

The difference in u-area, is still strange. In our current grid at the join:

image

and in the calculating from the mom supergrid at the join:

image

I think the second one is correct. (Amusingly the files in '/g/data/ik11/grids/' also have one of each option)

This figure kind of demonstrates it (viewed straight down on the north pole):

image

u-area is the sum of the areas of a quarter of each grid cell which surrounds each u-point (i.e. NE corner point). The areas north of the join (The thick line) are smaller than those south of the join with a notable change. So the u-areas at the join are the average of the smaller and larger cell areas which exist either side of the join, and therefore create the 'step' shown in the second figure.

As an aside, I stumbled across this code which Nic wrote and does what we want:

https://github.com/COSIMA/esmgrids/blob/130c8d00774c624c196fc2ca15f5cad78c4cdbc1/esmgrids/mom_grid.py

@aekiss
Copy link
Contributor

aekiss commented Mar 6, 2024

Does this mean the coordinates move away from regular 1 degree spacing?

yes, the meridional (y) spacing is refined to 1/3° near the equator on the 1° grid.

Note that the y spacing is not regular at most latitudes at 0.25° or 0.1° either, because a Mercator grid is used between 65S and 65N to keep the cell aspect ratio equal to 1. See https://github.com/COSIMA/ACCESS-OM2-1-025-010deg-report/blob/master/figures/grid/grid.ipynb

p.s. Its also odd that the patterning below 65N is there in the uarea calc and not the tarea calc.

yes, indeed

@aekiss
Copy link
Contributor

aekiss commented Mar 6, 2024

I think I agree with your analysis @anton-seaice, so long as the 65°N junction between the Mercator and tripolar grids falls on U points - is that the case?

@anton-seaice
Copy link
Contributor

anton-seaice commented Mar 6, 2024

I think I agree with your analysis @anton-seaice, so long as the 65°N junction between the Mercator and tripolar grids falls on U points - is that the case?

Good point - I checked - all three resolutions have the join at u-points

@ezhilsabareesh8
Copy link
Contributor Author

The angle ANGLET in the current implementation is calculated directly from the super grid as Refer here

ANGLE = np.deg2rad(in_superGridFile['angle_dx'][2::2, 2::2]) ; ANGLET= np.deg2rad(in_superGridFile['angle_dx'][1::2, 1::2])

which is consistent with the OM2 grid.nc file. However, it appears that this differs from the weighted average approach employed in CICE6 here.

Should we consider adapting the angleT calculation method to utilize the weighted average approach used in CICE6?

TLON = np.deg2rad(in_superGridFile['x'][1::2, 1::2])

HTN = in_superGridFile['dx'] * 100.0 # convert to cm
HTN = HTN[1::2, ::2] + HTN[1::2, 1::2]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Screenshot 2024-03-07 at 4 15 14 pm

In the mom_grid, dx has units ('nyp','nx'). These correspond to the edges in y and the middle in x. So, the y coordinate (i.e. rows) should be 2::2

…t edge respectively. Therefore for cice we want the index to start from 2.
# add an extra column at the end. Also u-cells cross the
# tri-polar fold so add an extra row at the top.
area_ext = np.append(AREA[:], AREA[:, 0:1], axis=1)
area_ext = np.append(area_ext[:], area_ext[-1:, :], axis=0)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this line only works because the top row is symmetrical? It doesn't actual seem to fold it... ?

@anton-seaice
Copy link
Contributor

I did a plot with every var and every resolution (left to right: 1deg, 0.25deg, 0.1deg). I masked difference which were less than 2e-6

  • Need to investigate: Angle, 0.1 ULAT/ULON
  • UAREA doesn't match well at any resolution, which is problematic.

I think we should merge this PR soon. We still have to add CF metadata and test with OM2, but these could be new issues to stop this PR totally blowing out? We also need to decide if we are sticking with np syntax (e.g. [2::2,1::2] or moving to written xarray style (e.g. ds.coarsen and ds.isel(x=slice(...)))

image

@anton-seaice
Copy link
Contributor

I turned on compression:

Before

-rw-r--r-- 1 as2285 tm70 742M Mar  7 16:27 grid_01deg.nc
-rw-r--r-- 1 as2285 tm70 119M Mar  7 16:27 grid_025deg.nc
-rw-r--r-- 1 as2285 tm70 8.3M Mar  7 16:27 grid_1deg.nc

After:

-rw-r--r-- 1 as2285 tm70  123M Mar  8 10:47 grid_01deg.nc
-rw-r--r-- 1 as2285 tm70   15M Mar  8 10:46 grid_025deg.nc
-rw-r--r-- 1 as2285 tm70 1019K Mar  8 10:46 grid_1deg.nc

I left chunking at the default sizes, which gives 0.1 deg chunks of 900, 1200 (y,x) but the other resolutions are full size

else:
nc.inputfile = f"{in_superGridPath}"
nc.timeGenerated = f"{datetime.now()}"
nc.created_by = f"{os.environ.get('USER')}"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nc.inputfile = f"{in_superGridPath}"
nc.timeGenerated = f"{datetime.now()}"
nc.created_by = f"{os.environ.get('USER')}"

These three lines can go outside the if/else as they are the same in both cases.

Can we add a hash for the input file? Similar to the payu manifests md5 hash?

@anton-seaice
Copy link
Contributor

anton-seaice commented Mar 8, 2024

I compared every resolution & every var between the version from this script, and the version in '/g/data/ik11/inputs/access-om2/input_20200530/'

  • ulat/ulon/angle for 0.25 degree resolution has been updated to be referencing the correct corner now
  • angleT in 0.01 degree looks wrong, as Ezhil noted
  • as discussed, uarea has changed at 1deg/0/25 deg resolution

image

(SRC: https://github.com/anton-seaice/sandbox/blob/743147926a9f7afabefc90276e45e8e02bcfa3e5/cice_grid_valiation.ipynb)

@anton-seaice
Copy link
Contributor

I suggest - to remove duplication in code and make the result more extendable we close this PR and re-write this using https://github.com/COSIMA/esmgrids.

The script can be based on https://github.com/COSIMA/access-om2/blob/master/tools/make_cice_grid.py but we need add version control, fix AngleT in esmgrids and add CF metadata.

Using esmgrids makes this easier to extend in the future (e.g. to make the ocean_grid file for analysis).

@anton-seaice
Copy link
Contributor

I suggest - to remove duplication in code and make the result more extendable we close this PR and re-write this using https://github.com/COSIMA/esmgrids.

The script can be based on https://github.com/COSIMA/access-om2/blob/master/tools/make_cice_grid.py but we need add version control, fix AngleT in esmgrids and add CF metadata.

Using esmgrids makes this easier to extend in the future (e.g. to make the ocean_grid file for analysis).

@aekiss thoughts ?

@anton-seaice anton-seaice mentioned this pull request Mar 20, 2024
@anton-seaice
Copy link
Contributor

Thanks for all your work on this Ezhil. I've taken our lessons and re-worked using esmgrids in a new PR :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants