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
96 changes: 96 additions & 0 deletions grid_generation/Gen_CICE_grid.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
"""
Script: generate_cice_grid.py
ezhilsabareesh8 marked this conversation as resolved.
Show resolved Hide resolved

Description:
This script generates a CICE grid from the MOM super grid provided in the input NetCDF file.

Contact:
Name: Ezhilsabareesh Kannadasan

Usage:
python generate_cice_grid.py <input_superGridPath> <output_file>
- input_superGridPath: Path to the MOM super grid NetCDF file.
- output_file: Path to the output CICE grid NetCDF file.
"""
import numpy as np
import xarray as xr
from netCDF4 import Dataset
import sys

def generate_cice_grid(in_superGridPath, output_file):
# Read input files
in_superGridFile = xr.open_dataset(in_superGridPath)

# 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

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

# MPI
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

HTE = in_superGridFile['dy'] * 100.0 # convert to cm
HTE = HTE[::2, 1::2] + HTE[1::2, 1::2]

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

# Close input files
in_superGridFile.close()

# Create a new NetCDF file
nc = Dataset(output_file, 'w', format='NETCDF4')

# Define dimensions
ny, nx = ULAT.shape
nc.createDimension('ny', ny)
nc.createDimension('nx', nx)

# Define variables
ulat = nc.createVariable('ulat', 'f8', ('ny', 'nx'))
ulon = nc.createVariable('ulon', 'f8', ('ny', 'nx'))
tlat = nc.createVariable('tlat', 'f8', ('ny', 'nx'))
tlon = nc.createVariable('tlon', 'f8', ('ny', 'nx'))
htn = nc.createVariable('htn', 'f8', ('ny', 'nx'))
hte = nc.createVariable('hte', 'f8', ('ny', 'nx'))
angle = nc.createVariable('angle', 'f8', ('ny', 'nx'))

# Add attributes
ulat.units = "radians"
ulat.title = "Latitude of U points"
ulon.units = "radians"
ulon.title = "Longitude of U points"
tlat.units = "radians"
tlat.title = "Latitude of T points"
tlon.units = "radians"
tlon.title = "Longitude of T points"
htn.units = "cm"
htn.title = "Width of T cells on North side."
ezhilsabareesh8 marked this conversation as resolved.
Show resolved Hide resolved
hte.units = "cm"
hte.title = "Width of T cells on East side."
ezhilsabareesh8 marked this conversation as resolved.
Show resolved Hide resolved
angle.units = "radians"
angle.title = "Rotation angle of U cells."

# Write data to variables
ulat[:] = ULAT
ulon[:] = ULON
tlat[:] = TLAT
tlon[:] = TLON
htn[:] = HTN
hte[:] = HTE
angle[:] = ANGLE

# Close the file
nc.close()

print("NetCDF file created successfully.")

if __name__ == "__main__":
if len(sys.argv) != 3:
print("Usage: python script.py <input_superGridPath> <output_file>")
ezhilsabareesh8 marked this conversation as resolved.
Show resolved Hide resolved
sys.exit(1)

input_superGridPath = sys.argv[1]
output_file = sys.argv[2]

generate_cice_grid(input_superGridPath, output_file)