Compare CICE input grid and mesh files

CICE_mesh_mask_cmp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9b8c3db6-6ab7-428e-add4-3d950457b84b",
   "metadata": {},
   "source": [
    "# CICE is throwing an error that the masks are inconsistent between the provided and internally generated mesh files\n",
    "\n",
    "- In the mesh file, the mask is generated using `/g/data/ik11/inputs/access-om2/input_20201102/mom_1deg/ocean_mask.nc`\n",
    "- The ice mask, `kmt.nc`, is generated just by renaming `ocean_mask.nc`. In the current configuration `/g/data/ik11/inputs/access-om2/input_20201102/cice_1deg/kmt.nc` is used (which symlinks to `/g/data/ik11/inputs/access-om2/input_20201022/cice_1deg/kmt.nc`)\n",
    "\n",
    "## Useful links\n",
    "- https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/doc/source/user_guide/ug_implementation.rst#L144"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf8c091e-7e78-4597-ac95-d7b906260d4a",
   "metadata": {},
   "source": [
    "`ocean_mask.nc` and `kmt.nc` should match. Do they?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "683eb260-8bdd-4c9f-80f6-94f0a12e7e2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "34a9ed26-d91f-474c-8eea-edaf56bdc645",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocean_mask = xr.open_dataset(\n",
    "    \"/g/data/ik11/inputs/access-om2/input_20201102/mom_1deg/ocean_mask.nc\"\n",
    ")\n",
    "ice_kmt = xr.open_dataset(\n",
    "    \"/g/data/ik11/inputs/access-om2/input_20201022/cice_1deg/kmt.nc\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "80fd274a-6a88-4493-903a-ff8b103d14f6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(ocean_mask[\"mask\"] == ice_kmt[\"kmt\"]).all().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e84fda73-c30c-436c-b4c7-269227ec4936",
   "metadata": {},
   "source": [
    "So these mask match, but the internally generated mask in CICE doesn't match the mask in the mesh file. Why? Are the mesh and kmt masks consistent?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25764e7e-e695-47e8-9c63-9b8dff4a13db",
   "metadata": {},
   "source": [
    "Add lon and lat locations to kmt file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4845d58a-88d6-43d9-993e-383b7158a825",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b96ec05a-f710-4f6e-a178-c74ba814254e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ice_grid = xr.open_dataset(\n",
    "    \"/g/data/ik11/inputs/access-om3/0.x.0/grid.nc\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "469ace9e-de87-4a92-b770-13650036d26d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tcells are cell centres\n",
    "ice_kmt = ice_kmt.assign_coords({\n",
    "    \"lon\": ([\"ny\", \"nx\"], (np.rad2deg(ice_grid.tlon.values) + 360) % 360),\n",
    "    \"lat\": ([\"ny\", \"nx\"], np.rad2deg(ice_grid.tlat.values))\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2a8f6299-715b-4183-b6ac-34468b75e6c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = xr.open_dataset(\n",
    "    \"/g/data/ik11/inputs/access-om3/0.x.0/access-om2-1deg-ESMFmesh.nc\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "254801a8-2461-42c4-afb3-95b4aa54cd47",
   "metadata": {},
   "source": [
    "Check that mask in mesh and kmt files match"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f1abaf5b-c12b-4119-a153-583c50412fa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_mask_values(kmt, mesh):\n",
    "    \"\"\"Check mask values between kmt and mesh\"\"\"\n",
    "\n",
    "    eps = 1e-4\n",
    "    centerCoords = mesh.centerCoords.values\n",
    "\n",
    "    for elem in range(len(mesh.elementCount)):\n",
    "        lon, lat = centerCoords[elem]\n",
    "        \n",
    "        mesh_mask = mesh.elementMask.values[elem]\n",
    "        match_lat_lon = (abs(kmt.lon - lon) < eps) & (abs(kmt.lat - lat) < eps)\n",
    "        kmt_mask = int(kmt.kmt.where(match_lat_lon, drop=True).item())\n",
    "    \n",
    "        if mesh_mask != kmt_mask:\n",
    "            print(f\"mask mismatch at lon {lon}, lat {lat}\")\n",
    "            print(f\"  mesh: {mesh_mask}\")\n",
    "            print(f\"  kmt: {kmt_mask}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "df204a30-58fe-4625-bbe7-81967c28ccfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "check_mask_values(ice_kmt, mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "822e9721-c226-4d89-b359-99623f6f4ed3",
   "metadata": {},
   "source": [
    "How do things look for a configuration that actually runs?\n",
    "\n",
    "Need to read CICE bin files used in the CESM displaced pole configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "15c56e45-9ae5-4747-a93d-c047d18d6981",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_cice_grid_bin(file, nx, ny):\n",
    "    \"\"\"\n",
    "    Open a CICE grid binary file into an xarray Dataset\n",
    "\n",
    "    See:\n",
    "    - https://github.com/NCAR/CESM_MOM6_AQUA/blob/90690fd726123c717c55a4622de70d2243d1d4d7/Grid_Topo_Generation/gen_cesm_files.ipynb\n",
    "    - https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/cicecore/cicedyn/infrastructure/ice_read_write.F90#L540\n",
    "    - https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/cicecore/cicedyn/infrastructure/ice_grid.F90#L849\n",
    "    - https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/cicecore/cicedyn/infrastructure/ice_grid.F90#L387\n",
    "    \"\"\"\n",
    "\n",
    "    size = nx * ny\n",
    "    variables = [\"ulat\", \"ulon\", \"htn\", \"hte\", \"hus\", \"huw\", \"angle\"]\n",
    "\n",
    "    data_vars = {}\n",
    "    with open(file,'rb') as fobj:\n",
    "        for idx, variable in enumerate(variables):\n",
    "            data = np.fromfile(fobj, dtype=\"f8\", count=size).byteswap()\n",
    "            data = np.reshape(data, (nx, ny), order=\"F\").T\n",
    "            data_vars[variable] = ((\"ny\", \"nx\"), data)\n",
    "\n",
    "    return xr.Dataset(data_vars)\n",
    "\n",
    "def read_cice_topog_bin(file, nx, ny):\n",
    "    \"\"\"\n",
    "    Open a CICE topography/kmt binary file into an xarray Dataset\n",
    "    \n",
    "    See:\n",
    "    - https://github.com/NCAR/CESM_MOM6_AQUA/blob/90690fd726123c717c55a4622de70d2243d1d4d7/Grid_Topo_Generation/gen_cesm_files.ipynb\n",
    "    - https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/cicecore/cicedyn/infrastructure/ice_read_write.F90#L540\n",
    "    - https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/cicecore/cicedyn/infrastructure/ice_grid.F90#L849\n",
    "    - https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/cicecore/cicedyn/infrastructure/ice_grid.F90#L387\n",
    "    \"\"\"\n",
    "\n",
    "    size = nx * ny\n",
    "    with open(file,'rb') as fobj:\n",
    "        data = np.fromfile(fobj, dtype=\"i4\").byteswap()\n",
    "        data = np.reshape(data, (nx, ny), order=\"F\").T\n",
    "\n",
    "    return xr.Dataset({\"kmt\": ((\"ny\", \"nx\"), data)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "969dc776-aae0-4c96-9adc-8ea056cb1080",
   "metadata": {},
   "outputs": [],
   "source": [
    "nx_global = 320\n",
    "ny_global = 384\n",
    "\n",
    "ice_grid = read_cice_grid_bin(\n",
    "    \"/g/data/ik11/inputs/access-om3/0.x.0/cime/ocn/pop/gx1v6/grid/horiz_grid_20010402.ieeer8\",\n",
    "    nx_global,\n",
    "    ny_global\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2780042b-c5a6-4d13-85df-629c85e4ea40",
   "metadata": {},
   "source": [
    "Only ulat and ulon are available in this grid file, so we need to calculate the cell centres tlat and tlon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "636f66df-23af-4064-bdec-05740cabd0a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _halo_extrapolate(arr, ew_boundary_type, ns_boundary_type):\n",
    "    \"\"\"\n",
    "    Extrapolate to 'halo' cell as in CICE code ice_boundary.F90: subroutine ice_HaloExtrapolate.\n",
    "\n",
    "    nghost is currently hardcoded to 1\n",
    "\n",
    "    Stolen and adapted from \n",
    "    https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/configuration/tools/jra55_datasets/interp_jra55_ncdf_bilinear.py#L62\n",
    "    \"\"\"\n",
    "\n",
    "    nj, ni = arr.shape\n",
    "    \n",
    "    # E/W edges\n",
    "    if ew_boundary_type == 'cyclic':\n",
    "        arr[: ,0] = arr[:,-2] # -2, since -1 is ghost cell\n",
    "        arr[:,-1] = arr[:, 1] #  1, since  0 is ghost cell\n",
    "    else:\n",
    "        arr[:, 0] = 2.0 * arr[:, 1] - arr[:, 2]\n",
    "        arr[:,-1] = 2.0 * arr[:,-2] - arr[:,-3]\n",
    "\n",
    "    # South edge\n",
    "    if ns_boundary_type == 'cyclic':\n",
    "        arr[0,:] = arr[-2,:] # -2, since -1 is ghost cell\n",
    "    else:\n",
    "        arr[0,:] = 2.0 * arr[1,:] - arr[2,:]\n",
    "\n",
    "    # North edge treated a little different, depending on if boundary type is tripole\n",
    "    if ns_boundary_type == 'cyclic':\n",
    "        arr[-1,:] = arr[1,:] # 1, since 0 is ghost cell\n",
    "    elif (ns_boundary_type != 'cyclic'  and\n",
    "          ns_boundary_type != 'tripole' and\n",
    "          ns_boundary_type != 'tripoleT'):\n",
    "\n",
    "        arr[-1,:] = 2.0 * arr[-2,:] - arr[-3,:]\n",
    "    else:\n",
    "        pass\n",
    "\n",
    "    return arr\n",
    "    \n",
    "def get_tcells(ulat, ulon, ew_boundary_type, ns_boundary_type):\n",
    "    \"\"\"\n",
    "    Make TLAT/TLON from ULAT/ULON as in CICE code ice_grid.F90: subroutine Tlatlon\n",
    "    \n",
    "    - Tripole grid: \n",
    "        ew_boundary_type = \"cyclic\", \n",
    "        ns_boundary_type = \"tripole\" or \"tripoleT\"\n",
    "    - Displaced pole grid: \n",
    "        ew_boundary_type = \"cyclic\", \n",
    "        ns_boundary_type = \"open\"\n",
    "\n",
    "    Stolen and adapted from \n",
    "    https://github.com/ESCOMP/CICE/blob/4cb296c4003014fe57d6d00f86868a78a532fc95/configuration/tools/jra55_datasets/interp_jra55_ncdf_bilinear.py#L111\n",
    "    \"\"\"\n",
    "\n",
    "    ulatcos = np.cos(ulat)\n",
    "    ulatsin = np.sin(ulat)\n",
    "    uloncos = np.cos(ulon)\n",
    "    ulonsin = np.sin(ulon)\n",
    "\n",
    "    nj, ni = ulatcos.shape \n",
    "    nghost = 1\n",
    "    workshape = (nj+2*nghost, ni+2*nghost)\n",
    "\n",
    "    ulatcos1 = np.zeros(workshape, dtype='f8')\n",
    "    ulatsin1 = np.zeros(workshape, dtype='f8')\n",
    "    uloncos1 = np.zeros(workshape, dtype='f8')\n",
    "    ulonsin1 = np.zeros(workshape, dtype='f8')\n",
    "\n",
    "    # Fill middle of work arrays\n",
    "    ulatcos1[1:nj+1,1:ni+1] = ulatcos\n",
    "    ulatsin1[1:nj+1,1:ni+1] = ulatsin\n",
    "    uloncos1[1:nj+1,1:ni+1] = uloncos\n",
    "    ulonsin1[1:nj+1,1:ni+1] = ulonsin\n",
    "\n",
    "    # Fill halos\n",
    "    ulatcos1 = _halo_extrapolate(ulatcos1, ew_boundary_type, ns_boundary_type)\n",
    "    ulatsin1 = _halo_extrapolate(ulatsin1, ew_boundary_type, ns_boundary_type)\n",
    "    uloncos1 = _halo_extrapolate(uloncos1, ew_boundary_type, ns_boundary_type)\n",
    "    ulonsin1 = _halo_extrapolate(ulonsin1, ew_boundary_type, ns_boundary_type)\n",
    "\n",
    "    # Calculations as in ice_grid.F90: Tlatlon\n",
    "    x = uloncos1 * ulatcos1\n",
    "    y = ulonsin1 * ulatcos1\n",
    "    z = ulatsin1\n",
    "    \n",
    "    tx = (\n",
    "        x[0:nj, 0:ni] +   # x1\n",
    "        x[0:nj, 1:ni+1] + # x2\n",
    "        x[1:nj+1, 0:ni] + # x3\n",
    "        x[1:nj+1, 1:ni+1] # x4\n",
    "    ) / 4\n",
    "    ty = (\n",
    "        y[0:nj, 0:ni] +   # y1\n",
    "        y[0:nj, 1:ni+1] + # y2\n",
    "        y[1:nj+1, 0:ni] + # y3\n",
    "        y[1:nj+1, 1:ni+1] # y4\n",
    "    ) / 4\n",
    "    tz = (\n",
    "        z[0:nj, 0:ni] +   # z1\n",
    "        z[0:nj, 1:ni+1] + # z2\n",
    "        z[1:nj+1, 0:ni] + # z3\n",
    "        z[1:nj+1, 1:ni+1] # z4\n",
    "    ) / 4\n",
    "\n",
    "    da = np.sqrt(tx ** 2 + ty ** 2 + tz ** 2)\n",
    "    tz = tz / da\n",
    "    tlat = np.arcsin(tz)\n",
    "    tlon = np.arctan2(ty,tx)\n",
    "    \n",
    "    return tlat, tlon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b0cf65bc-c042-4c58-8234-3c53b62a05d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "tlat, tlon = get_tcells(\n",
    "    ice_grid.ulat,\n",
    "    ice_grid.ulon,\n",
    "    ew_boundary_type=\"cyclic\",\n",
    "    ns_boundary_type=\"open\"\n",
    ")\n",
    "ice_grid[\"tlat\"] = xr.DataArray(tlat, dims=(\"ny\", \"nx\"))\n",
    "ice_grid[\"tlon\"] = xr.DataArray(tlon, dims=(\"ny\", \"nx\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2f1ff295-6638-487a-9df9-23cf1baae930",
   "metadata": {},
   "outputs": [],
   "source": [
    "ice_kmt = read_cice_topog_bin(\n",
    "    \"/g/data/ik11/inputs/access-om3/0.x.0/cime/ocn/pop/gx1v6/grid/topography_20090204.ieeei4\",\n",
    "    nx_global,\n",
    "    ny_global\n",
    ")\n",
    "ice_kmt = ice_kmt.assign_coords({\n",
    "    \"lon\": ([\"ny\", \"nx\"], (np.rad2deg(ice_grid.tlon.values) + 360) % 360),\n",
    "    \"lat\": ([\"ny\", \"nx\"], np.rad2deg(ice_grid.tlat.values))\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "eacc201a-630d-4619-890b-b0f322c0bc06",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = xr.open_dataset(\n",
    "    \"/g/data/ik11/inputs/access-om3/0.x.0/cime/share/meshes/gx1v6_090205_ESMFmesh.nc\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5abbb9ad-d6a1-439f-bf8f-cc1a2d18e76b",
   "metadata": {},
   "source": [
    "Check that mask in mesh and kmt files match"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5d91e085-72ac-482b-9b60-cd637e4923a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "check_mask_values(ice_kmt > 0, mesh)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}