From cf5ffe36d94b0c0616211c792e641804963a151e Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Fri, 20 Feb 2026 12:17:21 +0100 Subject: [PATCH] make the oasis writer work for amip runs --- ocp_tool/oasis_writer.py | 124 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) diff --git a/ocp_tool/oasis_writer.py b/ocp_tool/oasis_writer.py index 9b00194..507e803 100644 --- a/ocp_tool/oasis_writer.py +++ b/ocp_tool/oasis_writer.py @@ -72,6 +72,9 @@ def write_oasis_grid_files( # Copy runoff files into OASIS files (skip for AMIP mode - no ocean coupling) if config.ocean.grid_name != 'AMIP': _append_runoff_to_oasis_files(config, file_types) + else: + # Add AMIP forcing-reader grid for SST/SIC forcing + _append_amip_grid_to_oasis_files(config, file_types) def _write_single_oasis_file( @@ -205,6 +208,127 @@ def _append_runoff_to_oasis_files(config: OCPConfig, file_types: list) -> None: print(f"\n {'='*50} \n") +def _append_amip_grid_to_oasis_files(config: OCPConfig, file_types: list) -> None: + """Append AMIP forcing-reader grid (360x180 regular) to OASIS files for SST/SIC forcing.""" + print(' Adding AMIP forcing-reader grid (360x180) to OASIS files') + + # Create regular lat-lon grid: 180 lats x 360 lons + nlats = 180 + nlons = 360 + + # Cell centers: lon = 0.5 to 359.5, lat = -89.5 to 89.5 + lons_1d = np.arange(0.5, 360.0, 1.0) # 0.5, 1.5, 2.5, ..., 359.5 + lats_1d = np.arange(-89.5, 90.0, 1.0) # -89.5, -88.5, ..., 89.5 + lons, lats = np.meshgrid(lons_1d, lats_1d) + + # Cell corners (4 corners per cell) + lons_corners = np.arange(0.0, 360.0, 1.0) # 0, 1, 2, ..., 359 + lats_corners = np.arange(-90.0, 91.0, 1.0) # -90, -89, ..., 90 + + # Build corner arrays: shape (nlats, nlons, 4) + corners_lon = np.zeros((nlats, nlons, 4)) + corners_lat = np.zeros((nlats, nlons, 4)) + + for i in range(nlats): + for j in range(nlons): + # Corners in counter-clockwise order + # Use 360° for the right edge of the last cell, not 0° + lon_right = lons_corners[j+1] if j < nlons-1 else 360.0 + corners_lon[i, j, 0] = lons_corners[j] # bottom-left + corners_lon[i, j, 1] = lon_right # bottom-right + corners_lon[i, j, 2] = lon_right # top-right + corners_lon[i, j, 3] = lons_corners[j] # top-left + + corners_lat[i, j, 0] = lats_corners[i] # bottom-left + corners_lat[i, j, 1] = lats_corners[i] # bottom-right + corners_lat[i, j, 2] = lats_corners[i+1] # top-right + corners_lat[i, j, 3] = lats_corners[i+1] # top-left + + # Calculate cell areas (in m^2) + earth_radius = 6371000.0 # meters + cell_area = np.zeros((nlats, nlons)) + for i in range(nlats): + lat_min = lats_corners[i] * np.pi / 180.0 + lat_max = lats_corners[i+1] * np.pi / 180.0 + dlon_rad = 1.0 * np.pi / 180.0 + # Area = R^2 * dlon * (sin(lat_max) - sin(lat_min)) + cell_area[i, :] = earth_radius**2 * dlon_rad * (np.sin(lat_max) - np.sin(lat_min)) + + # Write to each file type + for file_type in file_types: + oasis_file = config.output_paths.oasis / f'{file_type}.nc' + nc = Dataset(str(oasis_file), 'a') + + grid_name = 'AMIP' + xname = f'x_{grid_name}' + yname = f'y_{grid_name}' + + # Create dimensions + if xname not in nc.dimensions: + nc.createDimension(xname, nlons) + if yname not in nc.dimensions: + nc.createDimension(yname, nlats) + + if file_type == 'grids': + # Write lon/lat centers + lonname = f'{grid_name}.lon' + latname = f'{grid_name}.lat' + + if lonname not in nc.variables: + var_lon = nc.createVariable(lonname, 'float64', (yname, xname)) + var_lon.units = 'degrees_east' + var_lon.standard_name = 'Longitude' + var_lon[:] = lons + + if latname not in nc.variables: + var_lat = nc.createVariable(latname, 'float64', (yname, xname)) + var_lat.units = 'degrees_north' + var_lat.standard_name = 'Latitude' + var_lat[:] = lats + + # Write corners + crnname = f'crn_{grid_name}' + cloname = f'{grid_name}.clo' + claname = f'{grid_name}.cla' + + if crnname not in nc.dimensions: + nc.createDimension(crnname, 4) + + if cloname not in nc.variables: + var_clo = nc.createVariable(cloname, 'float64', (crnname, yname, xname)) + var_clo.units = 'degrees_east' + var_clo.standard_name = 'Longitude corners' + var_clo[:] = np.transpose(corners_lon, (2, 0, 1)) # (nlats, nlons, 4) -> (4, nlats, nlons) + + if claname not in nc.variables: + var_cla = nc.createVariable(claname, 'float64', (crnname, yname, xname)) + var_cla.units = 'degrees_north' + var_cla.standard_name = 'Latitude corners' + var_cla[:] = np.transpose(corners_lat, (2, 0, 1)) # (nlats, nlons, 4) -> (4, nlats, nlons) + + elif file_type == 'areas': + areaname = f'{grid_name}.srf' + if areaname not in nc.variables: + var_area = nc.createVariable(areaname, 'float64', (yname, xname)) + var_area.units = 'm2' + var_area.standard_name = 'Grid cell area' + var_area[:] = cell_area + + elif file_type == 'masks': + maskname = f'{grid_name}.msk' + if maskname not in nc.variables: + var_mask = nc.createVariable(maskname, 'int32', (yname, xname)) + var_mask.units = '1' + var_mask.standard_name = 'Grid mask (0=land, 1=ocean)' + # All land (0) for AMIP forcing reader + var_mask[:] = np.zeros((nlats, nlons), dtype=np.int32) + + nc.close() + + print(' AMIP forcing-reader grid added successfully') + print(f"\n {'='*50} \n") + + def interpolate_vegin_data( config: OCPConfig, grid: GaussianGrid