source: TOOLS/MOSAIX/Build_coordinates_mask.ipynb @ 5925

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "10fcba33",
   "metadata": {},
   "source": [
    "# Constructs ```coordinates_mask``` file\n",
    "The ```coordinates_mask``` file is needed to run MOSAIX, to build interpolation weights for the coupled model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1cbb871",
   "metadata": {},
   "source": [
    "Olivier Marti - olivier.marti@lsce.ipsl - 10/06/2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "731d770c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy  as np\n",
    "import nemo\n",
    "import datetime, os, platform"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26422e3f",
   "metadata": {},
   "source": [
    "# Input files, type of ORCA grid, and output file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "deb66d3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#model   = 'eORCA1.2'\n",
    "model   = 'ORCA2.3'\n",
    "#model   = 'ORCA2.4'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2475c8b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "if model == 'eORCA1.2' :\n",
    "    n_coord = 'eORCA_R1_coordinates_v1.0.nc' \n",
    "    n_bathy = 'eORCA_R1_bathy_meter_v2.2.nc'\n",
    "    #n_coord = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_coordinates_v1.0.nc'\n",
    "    #n_bathy = 'https://vesg.ipsl.upmc.fr/thredds/dodsC/work/p86mart/GRAF/DATA/eORCA_R1_bathy_meter_v2.2.nc'\n",
    "    nperio  = 6\n",
    "    n_out   = 'eORCA_R1_coordinates_mask_test.nc'\n",
    "\n",
    "if model == 'ORCA2.3' :\n",
    "    nperio  = 4\n",
    "    n_coord = 'coordinates_ORCA2_LIM3_PISCES.nc' #'ORCA2.3_coordinates.nc' #'coordinates.nc' # \n",
    "    n_bathy = 'bathy_meter_ORCA2_LIM3_PISCES.nc' #'ORCA2.3_bathy_meter.nc' #'bathy_meter_closea_topo.nc' # \n",
    "    n_out   = 'My_test.nc'\n",
    "    \n",
    "if model == 'ORCA2.4' :\n",
    "    nperio  = 4\n",
    "    n_coord = 'ORCA2.4_coordinates.nc'\n",
    "    n_bathy = 'ORCA2.4_bathy_meter.nc'\n",
    "    n_out   = 'ORCA2.4_coordinates_mask_test.nc'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cd952a6",
   "metadata": {},
   "source": [
    "## Open input files\n",
    "Do not decode time which might be buggy in soem files\n",
    "Suppress singleton dimensions if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0b7c7279",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "f_coord = xr.open_dataset (n_coord, decode_times=False).squeeze()\n",
    "f_bathy = xr.open_dataset (n_bathy, decode_times=False).squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9fb92919",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "success\n"
     ]
    }
   ],
   "source": [
    "# Suppress time if necessary\n",
    "try    :\n",
    "    del f_coord['time']\n",
    "    print ('success')\n",
    "except :\n",
    "    pass\n",
    "    print ('failed to suppress time')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d61fb6b2",
   "metadata": {},
   "source": [
    "## Creates masks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "299c0a66",
   "metadata": {},
   "source": [
    "### Read bathymetry "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9dfd90e6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exception\n"
     ]
    }
   ],
   "source": [
    "Bathymetry = f_bathy['Bathymetry']\n",
    "\n",
    "try :\n",
    "    Bathymetry = Bathymetry.rename ({'nav_lat':'nav_lat_grid_T', 'nav_lon':'nav_lon_grid_T'})\n",
    "    print ('Normal case')\n",
    "except : \n",
    "    print ('Exception')\n",
    "    nav_lon_grid_T = f_bathy ['nav_lon']\n",
    "    nav_lat_grid_T = f_bathy ['nav_lat']\n",
    "    Bathymetry = xr.DataArray (Bathymetry, coords = { \"nav_lat_grid_T\": ([\"y\", \"x\"], nav_lat_grid_T),\n",
    "                                                      \"nav_lon_grid_T\": ([\"y\", \"x\"], nav_lon_grid_T) } )\n",
    "    \n",
    "Bathymetry = Bathymetry.rename ({'y':'y_grid_T', 'x':'x_grid_T'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2606d14b",
   "metadata": {},
   "outputs": [],
   "source": [
    "if model == 'ORCA2.3' :\n",
    "    # orca_r2: Gibraltar strait open\n",
    "    Bathymetry[101,139] = 284.\n",
    "    # orca_r2: Bab el Mandeb strait open\n",
    "    Bathymetry[87,159] = 137."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71c9c821",
   "metadata": {},
   "source": [
    "### Determine which periodicity is needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7e7a49d5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nperio specified :  4\n"
     ]
    }
   ],
   "source": [
    "jpj, jpi = Bathymetry.shape\n",
    "try : \n",
    "    if nperio != None :\n",
    "        print ('nperio specified : ', nperio)\n",
    "except :\n",
    "    nperio = None\n",
    "    if jpi ==  182 : nperio = 4 # ORCA2. \\!/ We choose legacy orca2\n",
    "    if jpi ==  362 : nperio = 6 # ORCA1.\n",
    "    if jpi == 1442 : nperio = 6 # ORCA025.\n",
    "    #\n",
    "    if nperio == None :\n",
    "        sys.exit ('nperio not found, and cannot by guessed' )\n",
    "    else :\n",
    "        print ('nperio set as {:d} (deduced from data size {:d} {:d})'.format (nperio, jpj, jpi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "34cdcfe4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Bathymetry = nemo.lbc (Bathymetry, nperio=nperio, cd_type='T')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0da91538",
   "metadata": {},
   "source": [
    "### Creates ```mask_T```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "49a654ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='int8')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "327d48ea",
   "metadata": {},
   "source": [
    "### Creates other masks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f0b656a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_U = mask_T * mask_T.shift (x_grid_T=-1)\n",
    "mask_V = mask_T * mask_T.shift (y_grid_T=-1)\n",
    "mask_F = mask_T * mask_T.shift (y_grid_T=-1) * mask_T.shift (x_grid_T=-1) * mask_T.shift (y_grid_T=-1, x_grid_T=-1)\n",
    "mask_W = mask_T\n",
    "\n",
    "mask_U = nemo.lbc (mask_U, nperio=nperio, cd_type='U', nemo_4U_bug=True).astype (dtype='int8')\n",
    "mask_V = nemo.lbc (mask_V, nperio=nperio, cd_type='V').astype (dtype='int8')\n",
    "mask_F = nemo.lbc (mask_F, nperio=nperio, cd_type='F').astype (dtype='int8')\n",
    "mask_W = nemo.lbc (mask_W, nperio=nperio, cd_type='W').astype (dtype='int8')\n",
    "\n",
    "mask_U = mask_U.rename ( {'y_grid_T'      : 'y_grid_U'      , 'x_grid_T'      : 'x_grid_U', \n",
    "                          'nav_lat_grid_T': 'nav_lat_grid_U', 'nav_lon_grid_T': 'nav_lon_grid_U'} )\n",
    "mask_V = mask_V.rename ( {'y_grid_T'      : 'y_grid_V'      , 'x_grid_T'      : 'x_grid_V',\n",
    "                          'nav_lat_grid_T': 'nav_lat_grid_V', 'nav_lon_grid_T': 'nav_lon_grid_V'} )\n",
    "mask_W = mask_W.rename ( {'y_grid_T'      : 'y_grid_W'      , 'x_grid_T'      : 'x_grid_W',\n",
    "                          'nav_lat_grid_T': 'nav_lat_grid_W', 'nav_lon_grid_T': 'nav_lon_grid_W'} )\n",
    "mask_F = mask_F.rename ( {'y_grid_T'      : 'y_grid_F'      , 'x_grid_T'      : 'x_grid_F',\n",
    "                          'nav_lat_grid_T': 'nav_lat_grid_F', 'nav_lon_grid_T': 'nav_lon_grid_F'} )\n",
    "\n",
    "mask_T.name = 'mask_T'\n",
    "mask_U.name = 'mask_U'\n",
    "mask_V.name = 'mask_V'\n",
    "mask_F.name = 'mask_F'\n",
    "mask_W.name = 'mask_W'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb36da6",
   "metadata": {},
   "source": [
    "### Masks with duplicate points removed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9ce5946",
   "metadata": {},
   "outputs": [],
   "source": [
    "#maskutil_T = nemo.lbc_mask (mask_T, nperio=nperio, cd_type='T')\n",
    "#maskutil_U = nemo.lbc_mask (mask_U, nperio=nperio, cd_type='U')\n",
    "#maskutil_V = nemo.lbc_mask (mask_V, nperio=nperio, cd_type='V')\n",
    "#maskutil_F = nemo.lbc_mask (mask_F, nperio=nperio, cd_type='F')\n",
    "#maskutil_W = nemo.lbc_mask (mask_W, nperio=nperio, cd_type='W')\n",
    "\n",
    "#maskutil_T.name = 'maskutil_T'\n",
    "#maskutil_U.name = 'maskutil_U'\n",
    "#maskutil_V.name = 'maskutil_V'\n",
    "#maskutil_F.name = 'maskutil_F'\n",
    "#maskutil_W.name = 'maskutil_W'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "274e4321",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plus compact que ci-dessus, mais un peu incompréhensible .... !!!\n",
    "for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n",
    "    MaskName = 'mask_'     + cd_type\n",
    "    UtilName = 'maskutil_' + cd_type\n",
    "    locals()[UtilName] = nemo.lbc_mask (locals()[MaskName], nperio=nperio, cd_type=cd_type)\n",
    "    locals()[UtilName].name = UtilName"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69ec3daa",
   "metadata": {},
   "source": [
    "### Add attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f46e1c93",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_T.attrs    ['cell_measures'] = 'area: area_grid_T'\n",
    "mask_U.attrs    ['cell_measures'] = 'area: area_grid_U'\n",
    "mask_V.attrs    ['cell_measures'] = 'area: area_grid_V'\n",
    "mask_W.attrs    ['cell_measures'] = 'area: area_grid_W'\n",
    "mask_F.attrs    ['cell_measures'] = 'area: area_grid_F'\n",
    "\n",
    "maskutil_T.attrs['cell_measures'] = 'area: area_grid_T'\n",
    "maskutil_U.attrs['cell_measures'] = 'area: area_grid_U'\n",
    "maskutil_V.attrs['cell_measures'] = 'area: area_grid_V'\n",
    "maskutil_W.attrs['cell_measures'] = 'area: area_grid_W'\n",
    "maskutil_F.attrs['cell_measures'] = 'area: area_grid_F'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "580e2776",
   "metadata": {},
   "source": [
    "## Creates grid coordinates and surfaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a7ce1490",
   "metadata": {},
   "outputs": [],
   "source": [
    "nav_lon_grid_T = nemo.lbc (f_coord ['glamt'], nperio=nperio, cd_type='T')\n",
    "nav_lat_grid_T = nemo.lbc (f_coord ['gphit'], nperio=nperio, cd_type='T')\n",
    "nav_lon_grid_U = nemo.lbc (f_coord ['glamu'], nperio=nperio, cd_type='U',nemo_4U_bug=True)\n",
    "#nav_lon_grid_U = f_coord ['glamu']\n",
    "nav_lat_grid_U = nemo.lbc (f_coord ['gphiu'], nperio=nperio, cd_type='U',nemo_4U_bug=True)\n",
    "#nav_lat_grid_U = f_coord ['gphiu']\n",
    "nav_lon_grid_V = nemo.lbc (f_coord ['glamv'], nperio=nperio, cd_type='V')\n",
    "nav_lat_grid_V = nemo.lbc (f_coord ['gphiv'], nperio=nperio, cd_type='V')\n",
    "nav_lon_grid_W = nemo.lbc (f_coord ['glamt'], nperio=nperio, cd_type='W')\n",
    "nav_lat_grid_W = nemo.lbc (f_coord ['gphit'], nperio=nperio, cd_type='W')\n",
    "nav_lon_grid_F = nemo.lbc (f_coord ['glamf'], nperio=nperio, cd_type='F')\n",
    "nav_lat_grid_F = nemo.lbc (f_coord ['gphif'], nperio=nperio, cd_type='F')\n",
    "\n",
    "nav_lon_grid_T = nav_lon_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )\n",
    "nav_lat_grid_T = nav_lat_grid_T.rename( {'y':'y_grid_T', 'x':'x_grid_T'} )\n",
    "nav_lon_grid_U = nav_lon_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )\n",
    "nav_lat_grid_U = nav_lat_grid_U.rename( {'y':'y_grid_U', 'x':'x_grid_U'} )\n",
    "nav_lon_grid_V = nav_lon_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )\n",
    "nav_lat_grid_V = nav_lat_grid_V.rename( {'y':'y_grid_V', 'x':'x_grid_V'} )\n",
    "nav_lon_grid_W = nav_lon_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )\n",
    "nav_lat_grid_W = nav_lat_grid_W.rename( {'y':'y_grid_W', 'x':'x_grid_W'} )\n",
    "nav_lon_grid_F = nav_lon_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )\n",
    "nav_lat_grid_F = nav_lat_grid_F.rename( {'y':'y_grid_F', 'x':'x_grid_F'} )\n",
    "\n",
    "for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n",
    "    for dir_type in ['lat', 'lon'] :\n",
    "        coord_name = 'nav_' + dir_type + '_grid_' + cd_type\n",
    "        locals()[coord_name].encoding['_FillValue'] = None\n",
    "        locals()[coord_name].encoding['missing_value'] = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "516b4c5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "nav_lon_grid_T.name = 'nav_lon_grid_T'\n",
    "nav_lat_grid_T.name = 'nav_lat_grid_T'\n",
    "nav_lon_grid_U.name = 'nav_lon_grid_U'\n",
    "nav_lat_grid_U.name = 'nav_lat_grid_U'\n",
    "nav_lon_grid_V.name = 'nav_lon_grid_V'\n",
    "nav_lat_grid_V.name = 'nav_lat_grid_V'\n",
    "nav_lon_grid_F.name = 'nav_lon_grid_F'\n",
    "nav_lat_grid_F.name = 'nav_lat_grid_F'\n",
    "nav_lon_grid_W.name = 'nav_lon_grid_W'\n",
    "nav_lat_grid_W.name = 'nav_lat_grid_W'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cffad9b4",
   "metadata": {},
   "source": [
    "### Add attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "de0e6514",
   "metadata": {},
   "outputs": [],
   "source": [
    "nav_lon_grid_T.attrs['standard_name'] = 'longitude'\n",
    "nav_lon_grid_T.attrs['long_name']     = 'Longitude'\n",
    "nav_lon_grid_T.attrs['units']         = 'degrees_east'\n",
    "nav_lat_grid_T.attrs['standard_name'] = 'latitude'\n",
    "nav_lat_grid_T.attrs['long_name']     = 'Latitude'\n",
    "nav_lat_grid_T.attrs['units']         = 'degrees_north'\n",
    "\n",
    "nav_lon_grid_U.attrs['standard_name'] = 'longitude'\n",
    "nav_lon_grid_U.attrs['long_name']     = 'Longitude'\n",
    "nav_lon_grid_U.attrs['units']         = 'degrees_east'\n",
    "nav_lat_grid_U.attrs['standard_name'] = 'latitude'\n",
    "nav_lat_grid_U.attrs['long_name']     = 'Latitude'\n",
    "nav_lat_grid_U.attrs['units']         = 'degrees_north'\n",
    "\n",
    "nav_lon_grid_V.attrs['standard_name'] = 'longitude'\n",
    "nav_lon_grid_V.attrs['long_name']     = 'Longitude'\n",
    "nav_lon_grid_V.attrs['units']         = 'degrees_east'\n",
    "nav_lat_grid_V.attrs['standard_name'] = 'latitude'\n",
    "nav_lat_grid_V.attrs['long_name']     = 'Latitude'\n",
    "nav_lat_grid_V.attrs['units']         = 'degrees_north'\n",
    "\n",
    "nav_lon_grid_W.attrs['standard_name'] = 'longitude'\n",
    "nav_lon_grid_W.attrs['long_name']     = 'Longitude'\n",
    "nav_lon_grid_W.attrs['units']         = 'degrees_east'\n",
    "nav_lat_grid_W.attrs['standard_name'] = 'latitude'\n",
    "nav_lat_grid_W.attrs['long_name']     = 'Latitude'\n",
    "nav_lat_grid_W.attrs['units']         = 'degrees_north'\n",
    "\n",
    "nav_lon_grid_F.attrs['standard_name'] = 'longitude'\n",
    "nav_lon_grid_F.attrs['long_name']     = 'Longitude'\n",
    "nav_lon_grid_F.attrs['units']         = 'degrees_east'\n",
    "nav_lat_grid_F.attrs['standard_name'] = 'latitude'\n",
    "nav_lat_grid_F.attrs['long_name']     = 'Latitude'\n",
    "nav_lat_grid_F.attrs['units']         = 'degrees_north'\n",
    "\n",
    "nav_lon_grid_T.attrs['bounds'] = 'bounds_lon_grid_T'\n",
    "nav_lat_grid_T.attrs['bounds'] = 'bounds_lat_grid_T'\n",
    "nav_lon_grid_U.attrs['bounds'] = 'bounds_lon_grid_U'\n",
    "nav_lat_grid_U.attrs['bounds'] = 'bounds_lat_grid_U'\n",
    "nav_lon_grid_V.attrs['bounds'] = 'bounds_lon_grid_V'\n",
    "nav_lat_grid_V.attrs['bounds'] = 'bounds_lat_grid_V'\n",
    "nav_lon_grid_W.attrs['bounds'] = 'bounds_lon_grid_W'\n",
    "nav_lat_grid_W.attrs['bounds'] = 'bounds_lat_grid_W'\n",
    "nav_lon_grid_F.attrs['bounds'] = 'bounds_lon_grid_F'\n",
    "nav_lat_grid_F.attrs['bounds'] = 'bounds_lat_grid_F'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "8e4d997c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove zero values from areas\n",
    "# need to be define for the extended grid south of -80S\n",
    "# some point are undefined but you need to have e1 and e2 .NE. 0\n",
    "for cd_type in ['t', 'u', 'v', 'f'] :\n",
    "    for axis in ['1', '2'] :\n",
    "        coordName = 'e' + axis + cd_type\n",
    "        f_coord[coordName]=xr.where(f_coord[coordName] == 0.0, 1.0e2, f_coord[coordName])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "175bae75",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correct areas for straits\n",
    "if model == 'ORCA2.3' :\n",
    "    # orca_r2: Gibraltar    : e2u reduced to 20 km\n",
    "    f_coord['e2u'][101,138:140] = 20.e3\n",
    "    # orca_r2: Bab el Mandeb: e2u reduced to 30 km\n",
    "    #                         e1v reduced to 18 km\n",
    "    f_coord['e1v'][87,159] = 18.e3\n",
    "    f_coord['e2u'][87,159] = 30.e3\n",
    "    # orca_r2: Danish Straits : e2u reduced to 10 km\n",
    "    f_coord['e2u'][115,144:146] = 10.e3\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "4f6cfbca",
   "metadata": {},
   "outputs": [],
   "source": [
    "area_grid_T = nemo.lbc (f_coord ['e1t']*f_coord ['e2t'], nperio=nperio, cd_type='T')\n",
    "area_grid_U = nemo.lbc (f_coord ['e1u']*f_coord ['e2u'], nperio=nperio, cd_type='U',nemo_4U_bug=True)\n",
    "#area_grid_U = f_coord ['e1u']*f_coord ['e2u']\n",
    "area_grid_V = nemo.lbc (f_coord ['e1v']*f_coord ['e2v'], nperio=nperio, cd_type='V')\n",
    "area_grid_W = nemo.lbc (f_coord ['e1t']*f_coord ['e2t'], nperio=nperio, cd_type='W')\n",
    "area_grid_F = nemo.lbc (f_coord ['e1f']*f_coord ['e2f'], nperio=nperio, cd_type='F')\n",
    "\n",
    "area_grid_T.name = 'area_grid_T'\n",
    "area_grid_U.name = 'area_grid_U'\n",
    "area_grid_V.name = 'area_grid_V'\n",
    "area_grid_F.name = 'area_grid_F'\n",
    "area_grid_W.name = 'area_grid_W'\n",
    "\n",
    "area_grid_T = area_grid_T.rename ({'y':'y_grid_T', 'x':'x_grid_T'})\n",
    "area_grid_U = area_grid_U.rename ({'y':'y_grid_U', 'x':'x_grid_U'})\n",
    "area_grid_V = area_grid_V.rename ({'y':'y_grid_V', 'x':'x_grid_V'})\n",
    "area_grid_W = area_grid_W.rename ({'y':'y_grid_W', 'x':'x_grid_W'})\n",
    "area_grid_F = area_grid_F.rename ({'y':'y_grid_F', 'x':'x_grid_F'})\n",
    "\n",
    "area_grid_T.attrs['standard_name'] = 'cell_area'\n",
    "area_grid_T.attrs['units']         = 'm2'\n",
    "area_grid_T.attrs['coordinates']   = 'nav_lat_grid_T nav_lon_grid_T'\n",
    "area_grid_U.attrs['standard_name'] = 'cell_area'\n",
    "area_grid_U.attrs['units']         = 'm2'\n",
    "area_grid_U.attrs['coordinates']   = 'nav_lat_grid_U nav_lon_grid_U'\n",
    "area_grid_V.attrs['standard_name'] = 'cell_area'\n",
    "area_grid_V.attrs['units']         = 'm2'\n",
    "area_grid_V.attrs['coordinates']   = 'nav_lat_grid_V nav_lon_grid_V'\n",
    "area_grid_W.attrs['standard_name'] = 'cell_area'\n",
    "area_grid_W.attrs['units']         = 'm2'\n",
    "area_grid_W.attrs['coordinates']   = 'nav_lat_grid_W nav_lon_grid_W'\n",
    "area_grid_F.attrs['standard_name'] = 'cell_area'\n",
    "area_grid_F.attrs['units']         = 'm2'\n",
    "area_grid_F.attrs['coordinates']   = 'nav_lat_grid_F nav_lon_grid_F'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c6f7b08",
   "metadata": {},
   "source": [
    "## Construct bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "69c404d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_bounds (cdgrd) :\n",
    "    '''\n",
    "    Constructs lon/lat bounds\n",
    "    Bounds are numerated counter clockwise, from bottom left\n",
    "    See NEMO file OPA_SRC/IOM/iom.F90, ROUTINE set_grid_bounds, for more details\n",
    "    '''\n",
    "    # Define offset of coordinate representing bottom-left corner\n",
    "    if cdgrd in ['T', 'W'] : \n",
    "        icnr = -1 ; jcnr = -1\n",
    "        corner_lon = nav_lon_grid_F ; corner_lat = nav_lat_grid_F\n",
    "    if cdgrd == 'U'        : \n",
    "        icnr =  0 ; jcnr = -1\n",
    "        corner_lon = nav_lon_grid_V ; corner_lat = nav_lat_grid_V\n",
    "    if cdgrd == 'V'        : \n",
    "        icnr = -1 ; jcnr =  0\n",
    "        corner_lon = nav_lon_grid_U ; corner_lat = nav_lat_grid_U\n",
    "    if cdgrd == 'F'        : \n",
    "        icnr = -1 ; jcnr = -1\n",
    "        corner_lon = nav_lon_grid_T ; corner_lat = nav_lat_grid_T\n",
    "\n",
    "    y_grid, x_grid = corner_lon.shape ; nvertex = 4\n",
    "    dims = ['y_grid_' + cdgrd, 'x_grid_' + cdgrd, 'nvertex_grid_' + cdgrd]\n",
    "    \n",
    "    bounds_lon = xr.DataArray (np.zeros ((y_grid, x_grid, nvertex)), dims=dims)\n",
    "    bounds_lat = xr.DataArray (np.zeros ((y_grid, x_grid, nvertex)), dims=dims)\n",
    "    \n",
    "    idx = [(jcnr,icnr), (jcnr,icnr+1), (jcnr+1,icnr+1), (jcnr+1,icnr)]\n",
    "    \n",
    "    # Compute cell vertices that can be defined, \n",
    "    # and complete with periodicity\n",
    "    for nn in range (nvertex) :\n",
    "        bounds_lon[:,:,nn] = np.roll  (corner_lon, shift=(idx[nn]), axis=(-2,-1))\n",
    "        bounds_lat[:,:,nn] = np.roll  (corner_lat, shift=(idx[nn]), axis=(-2,-1))\n",
    "        bounds_lon[:,:,nn] = nemo.lbc (bounds_lon[:,:,nn], nperio=nperio, cd_type=cdgrd)\n",
    "        bounds_lat[:,:,nn] = nemo.lbc (bounds_lat[:,:,nn], nperio=nperio, cd_type=cdgrd)\n",
    "    \n",
    "    # Rotate cells at the north fold\n",
    "    if nperio >= 3 :\n",
    "        # Working array for location of northfold\n",
    "        z_fld = nemo.lbc (np.ones ((jpj, jpi)), nperio=nperio, cd_type=cdgrd, psgn=-1. )\n",
    "        z_fld = (z_fld == -1.0)\n",
    "        \n",
    "        for nn in range (nvertex) :\n",
    "            nr = np.mod (nn+2, 4)\n",
    "            bounds_lon[:,:,nn] = np.where (z_fld, bounds_lon[:,:,nr], bounds_lon[:,:,nn])\n",
    "            bounds_lat[:,:,nn] = np.where (z_fld, bounds_lat[:,:,nr], bounds_lat[:,:,nn])\n",
    "    \n",
    "    # Invert cells at the symmetric equator\n",
    "    if nperio == 2 :\n",
    "        for nn in range (nvertex) :\n",
    "            nr = np.mod (nn+2, 4)\n",
    "            bounds_lon [0, :, nn] = bounds_lon [0, :, nr]\n",
    "            bounds_lat [0, :, nn] = bounds_lat [0, :, nr]\n",
    "    \n",
    "    bounds_lon.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd\n",
    "    bounds_lat.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd\n",
    "    bounds_lon.attrs['units']       = 'degrees_east'\n",
    "    bounds_lat.attrs['units']       = 'degrees_north'\n",
    "    bounds_lon.name = 'bounds_lon_grid_' + cdgrd\n",
    "    bounds_lat.name = 'bounds_lat_grid_' + cdgrd\n",
    "\n",
    "    return bounds_lon, bounds_lat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "87a8f592",
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds_lon_grid_T, bounds_lat_grid_T = set_bounds ('T')\n",
    "bounds_lon_grid_U, bounds_lat_grid_U = set_bounds ('U')\n",
    "bounds_lon_grid_V, bounds_lat_grid_V = set_bounds ('V')\n",
    "bounds_lon_grid_W, bounds_lat_grid_W = set_bounds ('W')\n",
    "bounds_lon_grid_F, bounds_lat_grid_F = set_bounds ('F')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dbf22b3",
   "metadata": {},
   "source": [
    "## Save Data in a NetCDF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "ca87fa1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.Dataset ({\n",
    "    'mask_T'      : mask_T     ,\n",
    "    'mask_U'      : mask_U     ,\n",
    "    'mask_V'      : mask_V     ,\n",
    "    'mask_W'      : mask_W     ,\n",
    "    'mask_F'      : mask_F     ,\n",
    "    'area_grid_T' : area_grid_T,\n",
    "    'area_grid_U' : area_grid_U,\n",
    "    'area_grid_V' : area_grid_V,\n",
    "    'area_grid_W' : area_grid_W,\n",
    "    'area_grid_F' : area_grid_F,\n",
    "    'maskutil_T'  : maskutil_T ,\n",
    "    'maskutil_U'  : maskutil_U ,\n",
    "    'maskutil_V'  : maskutil_V ,\n",
    "    'maskutil_W'  : maskutil_W ,\n",
    "    'maskutil_F'  : maskutil_F ,\n",
    "    'bounds_lon_grid_T': bounds_lon_grid_T,\n",
    "    'bounds_lat_grid_T': bounds_lat_grid_T,\n",
    "    'bounds_lon_grid_U': bounds_lon_grid_U,\n",
    "    'bounds_lat_grid_U': bounds_lat_grid_U,\n",
    "    'bounds_lon_grid_V': bounds_lon_grid_V,\n",
    "    'bounds_lat_grid_V': bounds_lat_grid_V,\n",
    "    'bounds_lon_grid_W': bounds_lon_grid_W,\n",
    "    'bounds_lat_grid_W': bounds_lat_grid_W,\n",
    "    'bounds_lon_grid_F': bounds_lon_grid_F,\n",
    "    'bounds_lat_grid_F': bounds_lat_grid_F,\n",
    "})\n",
    "for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n",
    "    for dir_type in ['lat', 'lon'] :\n",
    "        coord_name = 'nav_' + dir_type + '_grid_' + cd_type\n",
    "        ds.coords[coord_name]=locals()[coord_name]\n",
    "\n",
    "ds.attrs['name']         = 'coordinates_mask'\n",
    "ds.attrs['description']  = 'coordinates and mask for MOSAIX'\n",
    "ds.attrs['title']        = 'coordinates_mask'\n",
    "ds.attrs['source']       = 'IPSL Earth system model'\n",
    "ds.attrs['group']        = 'ICMC IPSL Climate Modelling Center'\n",
    "ds.attrs['Institution']  = 'IPSL https.//www.ipsl.fr'\n",
    "ds.attrs['Model']        = model\n",
    "ds.attrs['timeStamp']    = '{:%Y-%b-%d %H:%M:%S}'.format (datetime.datetime.now ())\n",
    "ds.attrs['history']      = 'Build from ' + n_coord + ' and ' + n_bathy\n",
    "ds.attrs['directory']    = os.getcwd     ()\n",
    "ds.attrs['user']         = os.getlogin   ()\n",
    "ds.attrs['HOSTNAME']     = platform.node ()\n",
    "ds.attrs['Python']       = 'Python version: ' +  platform.python_version ()\n",
    "ds.attrs['xarray']       = 'xarray version: ' +  xr.__version__\n",
    "ds.attrs['OS']           = platform.system  ()\n",
    "ds.attrs['release']      = platform.release ()\n",
    "ds.attrs['hardware']     = platform.machine ()\n",
    "\n",
    "ds.attrs['SVN_Author']   = \"$Author:  $\"\n",
    "ds.attrs['SVN_Date']     = \"$Date:  $\"\n",
    "ds.attrs['SVN_Revision'] = \"$Revision:  $\"\n",
    "ds.attrs['SVN_Id']       = \"$Id: Build_coordinates_mask.py $\"\n",
    "ds.attrs['SVN_HeadURL']  = \"$HeadURL: http://forge.ipsl.jussieu.fr/igcmg/svn/TOOLS/MOSAIX/Build_coordinates_mask.py $\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bfcc709",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Why is this necessary ?\n",
    "ds['nav_lon_grid_T'].attrs['bounds'] = nav_lon_grid_T.attrs['bounds']\n",
    "ds['nav_lat_grid_T'].attrs['bounds'] = nav_lat_grid_T.attrs['bounds']\n",
    "ds['nav_lon_grid_U'].attrs['bounds'] = nav_lon_grid_U.attrs['bounds']\n",
    "ds['nav_lat_grid_U'].attrs['bounds'] = nav_lat_grid_U.attrs['bounds']\n",
    "ds['nav_lon_grid_V'].attrs['bounds'] = nav_lon_grid_V.attrs['bounds']\n",
    "ds['nav_lat_grid_V'].attrs['bounds'] = nav_lat_grid_V.attrs['bounds']\n",
    "ds['nav_lon_grid_W'].attrs['bounds'] = nav_lon_grid_W.attrs['bounds']\n",
    "ds['nav_lat_grid_W'].attrs['bounds'] = nav_lat_grid_W.attrs['bounds']\n",
    "ds['nav_lon_grid_F'].attrs['bounds'] = nav_lon_grid_F.attrs['bounds']\n",
    "ds['nav_lat_grid_F'].attrs['bounds'] = nav_lat_grid_F.attrs['bounds']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "41391c22",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.to_netcdf (n_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75143732",
   "metadata": {},
   "source": [
    "# That's all folk's !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88457886",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}