{ "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": "fb78f0e5", "metadata": {}, "source": [ "# Various Switches To Deal With Bugs and Implementation Choices" ] }, { "cell_type": "code", "execution_count": 2, "id": "2bd430d4", "metadata": {}, "outputs": [], "source": [ "# reproduce periodicity bug for U and nperio = 4\n", "nemoUbug = False\n", "# change metrics (and bathymetry) for straits\n", "straits = False\n", "# periodicity imposed on coordinates, metrics and areas\n", "coordperio = False\n", "# function used to compute mask_T from bathymetry\n", "maskbathynemo = False" ] }, { "cell_type": "markdown", "id": "26422e3f", "metadata": {}, "source": [ "# Input files, type of ORCA grid, and output file" ] }, { "cell_type": "code", "execution_count": 3, "id": "deb66d3d", "metadata": {}, "outputs": [], "source": [ "#model = 'eORCA1.2'\n", "#model = 'ORCA2.3'\n", "model = 'paleORCA'\n", "#model = 'ORCA2.4'\n", "#model = 'ORCA.05'http://localhost:8888/notebooks/Build_coordinates_mask_new_2409.ipynb#" ] }, { "cell_type": "code", "execution_count": 4, "id": "2475c8b7", "metadata": {}, "outputs": [], "source": [ "if model == 'eORCA1.2' :\n", " n_coord = 'eORCA1.2coordinates.nc' # 'eORCA_R1_coordinates_v1.0.nc' \n", " n_bathy = 'eORCA1.2bathy_meter.nc' # '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 = 'My_eORCA_R1_coordinates_mask_test.nc'\n", "\n", "if model == 'ORCA2.3' :\n", " nperio = 4\n", " n_coord = 'coordinates.nc' # 'ORCA2.3_coordinates.nc' # 'coordinates_ORCA2_LIM3_PISCES.nc' #\n", " n_bathy = 'bathy_meter_closea_topo.nc' # 'ORCA2.3_bathy_meter.nc' # 'bathy_meter_ORCA2_LIM3_PISCES.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'\n", " \n", "if model == 'paleORCA' :\n", " n_coord = 'coordinates_paleo.nc'\n", " n_bathy = 'bathy_meter_paleo.nc'\n", " nperio = 6\n", " n_out = 'My_coordinates_mask_paleo.nc'\n", "\n" ] }, { "cell_type": "markdown", "id": "8cd952a6", "metadata": {}, "source": [ "## Open input files\n", "Do not decode time which might be buggy in some files\n", "Suppress singleton dimensions if necessary" ] }, { "cell_type": "code", "execution_count": 5, "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": 6, "id": "9fb92919", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "failed to suppress time\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": 7, "id": "9dfd90e6", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Normal case\n" ] } ], "source": [ "Bathymetry = f_bathy['Bathymetry'].copy()\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": 8, "id": "2606d14b", "metadata": {}, "outputs": [], "source": [ "if straits :\n", " # Open straits for ORCA2 grid\n", " 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": 9, "id": "7e7a49d5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nperio specified : 6\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": 10, "id": "34cdcfe4", "metadata": {}, "outputs": [], "source": [ "Bathymetry = nemo.lbc (Bathymetry, nperio=nperio, cd_type='T')\n", "#Suppress ocean points at southernmost position for periodicity (3, 4, 5, 6)\n", "if nperio in (3 , 4 , 5 , 6) :\n", " Bathymetry[0,:] = 0.0 " ] }, { "cell_type": "markdown", "id": "0da91538", "metadata": {}, "source": [ "### Creates ```mask_T```" ] }, { "cell_type": "code", "execution_count": 11, "id": "49a654ef", "metadata": {}, "outputs": [], "source": [ "if maskbathynemo :\n", " # Use same formula as domzgr. Should be equivalent to Bathimetry > 0.0. \n", " mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='f4')\n", "else :\n", " mask_T = xr.where (Bathymetry > 0.0, 1, 0).astype (dtype='f4')" ] }, { "cell_type": "markdown", "id": "327d48ea", "metadata": {}, "source": [ "### Creates other masks" ] }, { "cell_type": "code", "execution_count": 12, "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=nemoUbug).astype (dtype='f4')\n", "mask_V = nemo.lbc (mask_V, nperio=nperio, cd_type='V').astype (dtype='f4')\n", "mask_F = nemo.lbc (mask_F, nperio=nperio, cd_type='F').astype (dtype='f4')\n", "mask_W = nemo.lbc (mask_W, nperio=nperio, cd_type='W').astype (dtype='f4')\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": 13, "id": "274e4321", "metadata": {}, "outputs": [], "source": [ "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].copy(), nperio=nperio, cd_type=cd_type).astype (dtype='f4')\n", " locals()[UtilName].name = UtilName\n", " locals()[UtilName].encoding['_FillValue'] = None" ] }, { "cell_type": "markdown", "id": "69ec3daa", "metadata": {}, "source": [ "### Add attributes" ] }, { "cell_type": "code", "execution_count": 14, "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": 15, "id": "a7ce1490", "metadata": {}, "outputs": [], "source": [ "if coordperio :\n", " nav_lon_grid_T = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='T')\n", " nav_lat_grid_T = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='T')\n", " nav_lon_grid_U = nemo.lbc (f_coord ['glamu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)\n", " nav_lat_grid_U = nemo.lbc (f_coord ['gphiu'].copy(), nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)\n", " nav_lon_grid_V = nemo.lbc (f_coord ['glamv'].copy(), nperio=nperio, cd_type='V')\n", " nav_lat_grid_V = nemo.lbc (f_coord ['gphiv'].copy(), nperio=nperio, cd_type='V')\n", " nav_lon_grid_W = nemo.lbc (f_coord ['glamt'].copy(), nperio=nperio, cd_type='W')\n", " nav_lat_grid_W = nemo.lbc (f_coord ['gphit'].copy(), nperio=nperio, cd_type='W')\n", " nav_lon_grid_F = nemo.lbc (f_coord ['glamf'].copy(), nperio=nperio, cd_type='F')\n", " nav_lat_grid_F = nemo.lbc (f_coord ['gphif'].copy(), nperio=nperio, cd_type='F')\n", "else :\n", " nav_lon_grid_T = f_coord ['glamt'].copy()\n", " nav_lat_grid_T = f_coord ['gphit'].copy()\n", " nav_lon_grid_U = f_coord ['glamu'].copy()\n", " nav_lat_grid_U = f_coord ['gphiu'].copy()\n", " nav_lon_grid_V = f_coord ['glamv'].copy()\n", " nav_lat_grid_V = f_coord ['gphiv'].copy()\n", " nav_lon_grid_W = f_coord ['glamt'].copy()\n", " nav_lat_grid_W = f_coord ['gphit'].copy()\n", " nav_lon_grid_F = f_coord ['glamf'].copy()\n", " nav_lat_grid_F = f_coord ['gphif'].copy()\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", "# remove _FillValue and missing_value\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": 16, "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": 17, "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": 18, "id": "8e4d997c", "metadata": {}, "outputs": [], "source": [ "# create new variables e1 e2 to keep f_coord the same\n", "for cd_type in ['t', 'u', 'v', 'f'] :\n", " for axis in ['1', '2'] :\n", " coordName = 'e' + axis + cd_type\n", " locals()[coordName]=f_coord[coordName].copy()\n", "# 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", " locals()[coordName]=xr.where(locals()[coordName] == 0.0, 1.0e2, locals()[coordName])" ] }, { "cell_type": "code", "execution_count": 19, "id": "175bae75", "metadata": {}, "outputs": [], "source": [ "# Correct areas for straits\n", "if straits :\n", " # ORCA R2 configuration\n", " if model == 'ORCA2.3' :\n", " # Gibraltar : e2u reduced to 20 km\n", " e2u[101,138:140] = 20.e3\n", " # Bab el Mandeb : e2u reduced to 30 km\n", " # e1v reduced to 18 km\n", " e1v[87,159] = 18.e3\n", " e2u[87,159] = 30.e3\n", " # Danish Straits: e2u reduced to 10 km\n", " e2u[115,144:146] = 10.e3\n", " # ORCA R1 configuration\n", " if model == 'eORCA1.2' :\n", " # Gibraltar : e2u reduced to 20 km\n", " e2u[240,281:283] = 20.e3\n", " # Bhosporus : e2u reduced to 10 km\n", " e2u[247,313:315] = 10.e3\n", " # Lombok : e1v reduced to 13 km\n", " e1v[163:165,43] = 13.e3\n", " # Sumba : e1v reduced to 8 km\n", " e1v[163:165,47] = 8.e3\n", " # Ombai : e1v reduced to 13 km\n", " e1v[163:165,52] = 13.e3\n", " # Timor Passage : e1v reduced to 20 km\n", " #e1v[163:165,55] = 20.e3\n", " # W Halmahera : e1v reduced to 30 km\n", " e1v[180:182,54] = 30.e3\n", " # E Halmahera : e1v reduced to 50 km\n", " e1v[180:182,57] = 50.e3\n", " # ORCA R05 configuration\n", " if model == 'ORCA.05' :\n", " # Reduced e2u at the Gibraltar Strait\n", " e2u[326,562:564] = 20.e3\n", " # Reduced e2u at the Bosphore Strait\n", " e2u[342,626:628] = 10.e3\n", " # Reduced e2u at the Sumba Strait\n", " e2u[231,92:94] = 40.e3\n", " # Reduced e2u at the Ombai Strait\n", " e2u[231,102] = 15.e3\n", " # Reduced e2u at the Palk Strait\n", " e2u[269,14] = 10.e3\n", " # Reduced e1v at the Lombok Strait\n", " e1v[231:233,86] = 10.e3\n", " # Reduced e1v at the Bab el Mandeb\n", " e1v[275,661] = 25.e3" ] }, { "cell_type": "code", "execution_count": 20, "id": "4f6cfbca", "metadata": {}, "outputs": [], "source": [ "if coordperio :\n", " area_grid_T = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='T')\n", " area_grid_U = nemo.lbc (e1u*e2u, nperio=nperio, cd_type='U',nemo_4U_bug=nemoUbug)\n", " area_grid_V = nemo.lbc (e1v*e2v, nperio=nperio, cd_type='V')\n", " area_grid_W = nemo.lbc (e1t*e2t, nperio=nperio, cd_type='W')\n", " area_grid_F = nemo.lbc (e1f*e2f, nperio=nperio, cd_type='F')\n", "else :\n", " area_grid_T = e1t*e2t\n", " area_grid_U = e1u*e2u\n", " area_grid_V = e1v*e2v\n", " area_grid_W = e1t*e2t\n", " area_grid_F = e1f*e2f\n", "\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'\n", "\n", "for cd_type in ['T', 'U', 'V', 'F', 'W'] :\n", " areaName = 'area_grid_' + cd_type\n", " locals()[areaName].encoding['_FillValue'] = None" ] }, { "cell_type": "markdown", "id": "3c6f7b08", "metadata": {}, "source": [ "## Construct bounds" ] }, { "cell_type": "code", "execution_count": 21, "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 = f_coord ['glamf'].copy() ; corner_lat = f_coord ['gphif'].copy()\n", " center_lon = f_coord ['glamt'].copy() ; center_lat = f_coord ['gphit'].copy()\n", " if cdgrd == 'U' : \n", " icnr = 0 ; jcnr = -1\n", " corner_lon = f_coord ['glamv'].copy() ; corner_lat = f_coord ['gphiv'].copy()\n", " center_lon = f_coord ['glamu'].copy() ; center_lat = f_coord ['gphiu'].copy()\n", " if cdgrd == 'V' : \n", " icnr = -1 ; jcnr = 0\n", " corner_lon = f_coord ['glamu'].copy() ; corner_lat = f_coord ['gphiu'].copy()\n", " center_lon = f_coord ['glamv'].copy() ; center_lat = f_coord ['gphiv'].copy()\n", " if cdgrd == 'F' : \n", " icnr = -1 ; jcnr = -1\n", " corner_lon = f_coord ['glamt'].copy() ; corner_lat = f_coord ['gphit'].copy()\n", " center_lon = f_coord ['glamf'].copy() ; center_lat = f_coord ['gphif'].copy()\n", " \n", " #y_grid, x_grid = corner_lon.shape ;\n", " nvertex = 4\n", " dims = ['y_grid_' + cdgrd, 'x_grid_' + cdgrd, 'nvertex_grid_' + cdgrd]\n", " \n", " bounds_lon = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)\n", " bounds_lat = xr.DataArray (np.zeros ((jpj, jpi, 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", " tmp = np.roll (corner_lon, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))\n", " bounds_lon[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]\n", " tmp = np.roll (corner_lat, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))\n", " bounds_lat[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]\n", " bounds_lon[:,:,nn] = nemo.lbc (bounds_lon[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)\n", " bounds_lat[:,:,nn] = nemo.lbc (bounds_lat[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)\n", " \n", " # Zero-size cells at closed boundaries if cell points provided,\n", " # otherwise they are closed cells with unrealistic bounds\n", " if not (nperio == 1 or nperio == 4 or nperio == 6) :\n", " for nn in range (nvertex) : \n", " bounds_lon[:,0,nn] = center_lon[:,0] # (West or jpni = 1), closed E-W\n", " bounds_lat[:,0,nn] = center_lat[:,0]\n", " bounds_lon[:,jpi,nn] = center_lon[:,jpi] # (East or jpni = 1), closed E-W\n", " bounds_lat[:,jpi,nn] = center_lat[:,jpi]\n", " if nperio != 2 :\n", " for nn in range (nvertex) :\n", " bounds_lon[0,:,nn] = center_lon[0,:] # (South or jpnj = 1), not symmetric\n", " bounds_lat[0,:,nn] = center_lat[0,:]\n", " if nperio < 3 :\n", " for nn in range (nvertex) :\n", " bounds_lon[jpj,:,nn] = center_lon[jpj,:] # (North or jpnj = 1), no north fold\n", " bounds_lat[jpj,:,nn] = center_lat[jpj,:]\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., nemo_4U_bug=nemoUbug)\n", " z_fld = np.repeat((z_fld == -1.0)[...,np.newaxis],4,axis=2)\n", " # circular shift of 2 indices in bounds third index\n", " bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)\n", " bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)\n", " bounds_lon[:,:,:] = np.where (z_fld, bounds_lon_tmp[:,:,:] , bounds_lon[:,:,:] )\n", " bounds_lat[:,:,:] = np.where (z_fld, bounds_lat_tmp[:,:,:] , bounds_lat[:,:,:] )\n", " \n", " # Invert cells at the symmetric equator\n", " if nperio == 2 :\n", " bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)\n", " bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)\n", " bounds_lon [0,:,:] = bounds_lon [0,:,:]\n", " bounds_lat [0,:,:] = bounds_lat [0,:,:]\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", " # remove _FillValue\n", " bounds_lon.encoding['_FillValue'] = None\n", " bounds_lat.encoding['_FillValue'] = None\n", "\n", " return bounds_lon, bounds_lat" ] }, { "cell_type": "code", "execution_count": 22, "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": 23, "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", "#replace nav_lon nav_lat with variables obtained from NEMO coordinates.nc file\n", "#by construction nav_lon nav_lat are those of the bathymetry \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$\"\n", "ds.attrs['SVN_HeadURL'] = \"$HeadURL$\"" ] }, { "cell_type": "code", "execution_count": 24, "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 }