MOM supergrid and other grid types

[dabail10]

There is some interest in the community for allowing CICE to read the MOM6 supergrid file directly. In particular, check out this thread on the community forum:

We would use a combination of grid_type and grid_format here. I would suggest adding a category of grid_format = 'momnc'. In most cases these are grid_type = 'tripole'. Currently we have to write kmt/grid files from the MOM supergrid file. This process involves creating SCRIP grid files.

NUOPC also brings the the possibility of mesh files. This is currently only used for prescribed ice.

[anton-seaice]

Hi Dave - who was requesting the mom supergrid? I cant see a thread on the forum?

I think this is a reasonable option for ACCESS too. Would CESM use it if was available? ( @DeniseWorthen same question for UFS really?)

Anyway - ill bring it up at the next consortium meeting if we have time.

[DeniseWorthen]

I'm not sure how much use this capability would be for UFS. The supergrid only provides the lat/lon locations, so you still need to be able to supply the land mask and calculate the angle. One issue w/ the land mask is that since MOM6 can change the land mask at run-time (using the topoedits), you could end up w/ an inconsistent land mask, depending on how you generate the mask for CICE.

In UFS, we have a utility which creates the required files for CICE using the super-grid and ancillary files. That same utility creates a lot of other required files for our coupled system. For example, it generates the mapped ocean mask used to generate ATM ICs, a set of ESMF weights used to post-process the tripole grid to rectilinear or to downscale a tripole ocean/ice restart to lower tripole resolution, a set of files used to generate the WW3 grid etc. We need all those files anyway, so I'm not sure how much use we would make of a direct-supergrid functionality in CICE.

[dabail10]

This came from a post on the CESM-CICE forum. Paul Hall was developing something to read the MOM supergrid. Is this true that the MOM supergrid does not have the angles? I find this hard to believe.

https://bb.cgd.ucar.edu/cesm/threads/custom-grids-with-cice.7707/#post-46369

[DeniseWorthen]

You can only retrieve anglet from the supergrid.

[dabail10]

Got it. I guess one has to interpolate to the U and V points. We normally supply ANGLE at the B-grid U points and then interpolate to the T-grid points. However, the supergrid is twice as many points in X and Y as the model grid. So a visual here:

Model grid:

------V-----Q
|           |
|     T     U   
|           |
-------------

Supergrid:

T------T------T
|      |      |
T----- T------T
|      |      |
T------T------T

So, ANGLET should be available at the U and V points?

[dabail10]

Here are the contents of a supergrid file for our 2/3 degree tripole grid (540x480). Note that the supergrid is 1080x960.

netcdf ocean_hgrid_221123 {
dimensions:
	nyp = 961 ;
	nxp = 1081 ;
	ny = 960 ;
	nx = 1080 ;
	string = 255 ;
variables:
	char tile(string) ;
	double y(nyp, nxp) ;
		y:units = "degrees" ;
	double x(nyp, nxp) ;
		x:units = "degrees" ;
	double dy(ny, nxp) ;
		dy:units = "meters" ;
	double dx(nyp, nx) ;
		dx:units = "meters" ;
	double area(ny, nx) ;
		area:units = "m2" ;
	double angle_dx(nyp, nxp) ;
		angle_dx:units = "meters" ;

// global attributes:
		:file_name = "ocean_hgrid.nc" ;
		:Description = "MOM6 2/3 degree tripolar grid (ORCA type)" ;
		:Author = "Fred Castruccio (fredc@ucar.edu)" ;
		:Created = "23/11/2022 14:41:36 -0700" ;
		:type = "MOM6 supergrid file" ;

[DeniseWorthen]

Can you reproduce the internally calculated angles (sin_rot,cos_rot) using the angle_dx? I have never found them to agree.

[dabail10]

I guess that doesn't surprise me. I haven't tried that yet. The code in CICE uses the vector interpolation from the center of the Earth and should be pretty accurate. I don't know how the supergrid angles are computed. We generally try to allow for some errors here. I think it is a threshold of 0.01 or something like that. How big are the errors you are talking about? I think the most important thing is that the grid variables are consistent between MOM6 and CICE6 internal variables and the mesh files for both.

[DeniseWorthen]

Our grid generation code computes the rotation angles from latCt and lonCt in the same way that the actual MOM code does at initialization, so they are identical, except for 1 deg, where there are 10-6 differences along the right hand side seam. The angles I get from using the angle_dx array are quite a bit different. I'll try to remember to generate a figure for you. It's never been clear to me what the angle_dx array in the super grid represents.

[dabail10]

Got it. I will ask for further clarification.

[anton-seaice]

I'm not sure how much use this capability would be for UFS. The supergrid only provides the lat/lon locations, so you still need to be able to supply the land mask and calculate the angle.

Yeah the landmask does make it hard. I think we would still need a mask file for cice (or use the mesh mask file used in the coupler).

One issue w/ the land mask is that since MOM6 can change the land mask at run-time (using the topoedits), you could end up w/ an inconsistent land mask, depending on how you generate the mask for CICE.

Does the MOM cap include a comparison to the ESMF mesh mask file? There is this line (i'm not a 100% sure what it does) in MOM to compare to the nuopc/cmeps mask and similar in CICE. So between the two checks would it catch this?

In UFS, we have a utility which creates the required files for CICE using the super-grid and ancillary files.

Is this on github somewhere? It would be good to compare :)

We generally try to allow for some errors here. I think it is a threshold of 0.01 or something like that.

This looks like its only checked for lat/lon. Not angle or area.

How big are the errors you are talking about? I think the most important thing is that the grid variables are consistent between MOM6 and CICE6 internal variables and the mesh files for both.

I started down this quest because the areas seems quite different between this CICE calculation (which assumes rectangular grid cells) & the supergrid areas (presumably caluclated using the surface of a sphere).

[DeniseWorthen]

The utility is within UFS-UTILS , where it is called cpld_gridgen

[dabail10]

It is true that there is a separate bathymetry file that would be used for the mask or the landfast ice parameterization. We will have to add code for that as well.

[dabail10]

For the supergrid I have learned the following:

tlat(i,j) = supergrid_lat(2i,2j)
tlon(i,j) = supergrid_lat(2i,2j)
angle = supergrid_angle_dx(2i,2j) <- this might be the angle wrt true north (checking now)

For the B-grid:

ulat(i,j) = supergrid_lat(2*i+1,2*j+1)
ulon(i,j) = supergrid_lon(2*i+1,2*j+1)

For the C-grid:

ulat(i,j) = supergrid_lat(2*i+1,2*j)
ulon(i,j) = supergrid_lon(2*i+1,2*j)
vlat(i,j) = supergrid_lat(2*i,2*j+1)
vlon(i,j) = supergrid_lon(2*i,2*j+1)

[DeniseWorthen]

This is what I get comparing asin(sin_rot) from MOM6 output (black) and angleT using the supergrid (red). So the magnitudes don't match, and there is a funny little bump at the poles on either side using the super-grid. I'd be happy if you could explain this to me, because it has always puzzled me!

[dabail10]

This is really helpful. I will pass it on to the MOM6 people here. I wonder if this is an issue with your grid generator. I know we have our own tool.

[dabail10]

I think this is the same thing. However, here the arcsin(sin_rot) gives me the bump at 90N. I am plotting the top row of the grid. Were you doing a line down the middle? Mine don't quite line up because sin_rot has coordinates of lat/lon.

[DeniseWorthen]

In this case, I am plotting what MOM produces itself at run time after reading the super-grid vs what just reading the angle array directory (taking into account the factor of two). MOM itself does not utilize the angle from the super-grid, instead initializing it in

https://github.com/mom-ocean/MOM6/blob/129e1bda02d454fb280819d1d87ae16347fd044c/src/initialization/MOM_shared_initialization.F90#L560

[dabail10]

Can you attach your script that makes the plot?

[DeniseWorthen]

I was plotting the top row of the 1/4deg grid. You need to get rid of the lat/lon axes in the mom6 output. Try this NCO command

ncks -C -v sin_rot mom6-output-file.nc sin_rot.nc

[DeniseWorthen]

I have a file called 'tripole.mx025.nc' where I've written the anglet retrieved from the supergrid.

The variable I read from the supergrid is ang(0:nx,0:ny), where nx,ny are the super grid dimensions and ni,nj are the reduced grid dimensions:

  do j = 1,nj
     do i = 1,ni
        i2 = 2*i ; j2 = 2*j
        angle(i,j) =   ang(i2,j2)*deg2rad
        !deg
        anglet(i,j) =  -ang(i2-1,j2-1)*deg2rad

....

Then using that file, in ferret:

yes? set data tripole.mx025.nc
yes? set data sin_rot.nc
yes? plot/j=1080 asin(sin_rot[d=2]),anglet[d=1]

[dabail10]

OK. I think this is the same plot for our 2/3 degree tripole grid.

[dabail10]

Note that the orange curve is from the supergrid (angle_dx).

[DeniseWorthen]

hm. puzzling. Where did you get the 2/3deg configuration? I didn't think MOM6 provided anything other than 1deg,1/2deg and 1/4deg.

[dabail10]

This was generated at NCAR by the ORCA grid generator.

[DeniseWorthen]

Could you put the supergrid, mask and bathymetry file somewhere on derecho for me?

[dabail10]

Everything is in my home directory:

/glade/u/home/dbailey

angle_mom.py
ocean_hgrid_221123.nc
sin_rot.nc

Dave

[DeniseWorthen]

OK, thanks. Let me take a look.

[dabail10]

Ok. I figured out the issue with my python script. I had the indices reversed. Here is the latest plot. The orange line is the supergrid and the blue line is what is computed within MOM. There is some difference near the pole, but this is intentional. Fred Castruccio here at NCAR does this computation for the supergrid using his own code. Arguably we should be reading in this angle. The other option is to add the angle computation code from MOM into CICE. However, I would prefer the supergrid angle (orange line) myself.

[DeniseWorthen]

But in both cases, you are reading a supergrid file that has been created by the ORCA grid generator?

[dabail10]

Yes.

[eclare108213]

Not sure if this is helpful, but I dug up the Smith et al 1995 paper describing POP's grids. This is the basis for the curvilinear coordinate grid calculations (originally) in CICE. NB this paper has theorems and proofs... It also talks about B grids and C grids. Their "CD grid" means something different: Continuously-differentiable Dipole grid -- but the great-circle calculations are what's in CICE (I think).
LA-UR-95-1146.pdf

[anton-seaice]

I have a semi-functioning draft of a change for this. @dabail10 @DeniseWorthen - if you can provide a mom_supergrid file you use, I can test them. Especially if you have any that are not tripolar grids ?

[dabail10]

Thanks for opening this. I have attached a MOM supergrid file for our 2/3 degree global grid. Ironically this was generated using the ORCA grid generator. :) The main issue is that there is a field called angle_dx here, which does not agree with the internally generated angles in MOM. So, technically we should ignore this field in the supergrid file and use the same code that MOM has for the angles.

ocean_hgrid_221123.nc.gz

[DeniseWorthen]

I have copies of our 4 tripole MOM6 cases on Derecho in /glade/derecho/scratch/worthen. These are all "as is" from GFDL.

1/4deg mom6.025_hgrid.nc
1/2deg mom6.050_hgrid.nc
1deg mom6.100_hgrid.nc
5deg mom6.500_hgrid.nc

[anton-seaice]

There is another interesting quirk here.

MOM6 uses the mom_mosaic areas for the t-cell/a-grid and b-grid areas but calculates them internally for the c-grid u & v areas. See NOAA-GFDL/MOM6#740. The two areas will be quite different near the poles.

It probably doesn't matter much if we read them from the mom_mosaic file or calculate them, as all the coupling is done on the a-grid points anyway. I think my preference is the first option. Do either of you have a preference?

[dabail10]

I would like to be consistent with the MOM computation. This is also true of the angles.

[anton-seaice]

Thanks - we should be able to get equality to single point precision for both areas and angles.

MOM6 has a method to change channel widths at run-time . We won't be able to match that if CHANNEL_CONFIG is set to anything other than none in MOM_input