1_Differential_heating_Meridional_transport.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chauffage différentiel et transport méridien d'énergie\n",
    "\n",
    "**Auteur : FERRY Frédéric (DESR/ENM/C3M) - septembre 2022**\n",
    "\n",
    "D'après Stull : https://geo.libretexts.org/Bookshelves/Meteorology_and_Climate_Science/Book%3A_Practical_Meteorology_(Stull)/11%3A_General_Circulation/11.02%3A_Section_3-\n",
    "\n",
    "Le fichier de données au format netcdf (moyennes mensuelles CERES EBAF 4.1) doit être téléchargé et placé dans le répertoire ceres : https://asdc.larc.nasa.gov/project/CERES/CERES_EBAF_Edition4.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Dossiers des données, figures et animations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_data='./ceres/'\n",
    "dir_figs='./figs/'\n",
    "dir_anim='./anim/'\n",
    "if not os.path.exists(dir_figs):\n",
    "    os.makedirs(dir_figs)\n",
    "if not os.path.exists(dir_anim):\n",
    "    os.makedirs(dir_anim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "file=\"CERES_EBAF_Edition4.1_200003-202111.nc\"\n",
    "data=xr.open_dataset(dir_data+file)\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Selection de la période 2001-2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "year1=input(\"année de départ : \")\n",
    "year2=input(\"année de fin : \")"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Selection des données sur la période choisie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "data=data.sel(time=slice(year1,year2))\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bilan radiatif au sommet de l'atmosphère"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Tableaux des mois et des saisons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "months=['January','February','March','April', 'May', 'June',\n",
    "        'July', 'August', 'September', 'October', 'November', 'December']\n",
    "seasons = ['DJF','JJA','MAM','SON']"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Fonction de conversion des longitudes 0 - 360 en -180 - 180."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def lonflip(da):\n",
    "    lon_name = 'lon'\n",
    "    da['_longitude_adjusted'] = xr.where(\n",
    "        da[lon_name] > 180,\n",
    "        da[ lon_name] - 360,\n",
    "        da[lon_name])\n",
    "    da = (\n",
    "        da\n",
    "        .swap_dims({lon_name: '_longitude_adjusted'})\n",
    "        .sel(**{'_longitude_adjusted': sorted(da._longitude_adjusted)})\n",
    "        .drop(lon_name))\n",
    "    da = da.rename({'_longitude_adjusted': lon_name})\n",
    "    return da"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul des climatologies des données au sommet de l'atmosphère."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Rayonnement solaire incident\n",
    "solar_mean = lonflip(data['solar_mon'].mean('time')) # moyenne annuelle\n",
    "solar_mean_z = solar_mean.mean(axis=1)  # moyenne annuelle et zonale\n",
    "solar_mean_month = lonflip(data['solar_mon'].groupby('time.month').mean('time')) # moyenne mensuelle\n",
    "solar_mean_month_z = solar_mean_month.mean(axis=2)  # moyenne mensuelle et zonale\n",
    "solar_mean_season = lonflip(data['solar_mon'].groupby('time.season').mean('time')) # moyenne saisonnière\n",
    "solar_mean_sz = solar_mean_season.mean(axis=2)  # moyenne saisonnière et zonale\n",
    "\n",
    "# Rayonnement visible (solaire réfléchi)\n",
    "toa_sw_mean = lonflip(data['toa_sw_all_mon'].mean('time'))  # moyenne annuelle\n",
    "toa_sw_mean_z = toa_sw_mean.mean(axis=1) # moyenne annuelle et zonale\n",
    "toa_sw_mean_month = lonflip(data['toa_sw_all_mon'].groupby('time.month').mean('time')) # moyenne mensuelle\n",
    "toa_sw_mean_month_z = toa_sw_mean_month.mean(axis=2)  # moyenne mensuelle et zonale\n",
    "toa_sw_mean_season = lonflip(data['toa_sw_all_mon'].groupby('time.season').mean('time')) # moyenne saisonnière\n",
    "toa_sw_mean_sz = toa_sw_mean_season.mean(axis=2)  # moyenne saisonnière et zonale\n",
    "\n",
    "# Rayonnement IR émis\n",
    "toa_lw_mean = lonflip(data['toa_lw_all_mon'].mean('time'))  # moyenne annuelle\n",
    "toa_lw_mean_z = toa_lw_mean.mean(axis=1) # moyenne annuelle et zonale\n",
    "toa_lw_mean_month = lonflip(data['toa_lw_all_mon'].groupby('time.month').mean('time')) #  moyenne mensuelle\n",
    "toa_lw_mean_month_z = toa_lw_mean_month.mean(axis=2) # moyenne mensuelle et zonale\n",
    "toa_lw_mean_season = lonflip(data['toa_lw_all_mon'].groupby('time.season').mean('time'))  # moyenne saisonnière\n",
    "toa_lw_mean_sz = toa_lw_mean_season.mean(axis=2) # moyenne saisonnière et zonale\n",
    "\n",
    "# Royonnement net TOA\n",
    "toa_net_mean = lonflip(data['toa_net_all_mon'].mean('time'))# moyenne annuelle\n",
    "toa_net_mean_z = toa_net_mean.mean(axis=1) # moyenne annuelle et zonale\n",
    "toa_net_mean_month = lonflip(data['toa_net_all_mon'].groupby('time.month').mean('time')) # moyenne mensuelle\n",
    "toa_net_mean_month_z = toa_net_mean_month.mean(axis=2) # moyenne mensuelle et zonale\n",
    "toa_net_mean_season = lonflip(data['toa_net_all_mon'].groupby('time.season').mean('time')) # moyenne saisonnière\n",
    "toa_net_mean_sz = toa_net_mean_season.mean(axis=2) # moyenne saisonnière et zonale\n",
    "\n",
    "lat=solar_mean.lat.values\n",
    "lon=solar_mean.lon.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bilan radiatif au sommet de l'atmosphère en moyenne zonale"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Bilan radiatif en moyenne annuelle et zonale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,8))\n",
    "fig.suptitle('Top Of Atmosphere Radiative Budget - zonal and annual mean', fontsize=16)\n",
    "plt.title('CERES '+year1+'-'+year2, loc='center')\n",
    "\n",
    "plt.xlabel('Latitude')\n",
    "plt.ylabel('Radiation ($W/m^2$)')\n",
    "plt.xticks(np.arange(-90, 120, 30), ('90S', '60S', '30S', 'Eq', '30N', '60N', '90N'))\n",
    "plt.plot(lat, solar_mean_z, label='Incoming Solar Radiation', color='red', linewidth=2)\n",
    "plt.plot(lat, toa_sw_mean_z, label='Reflected Solar Radiation', color='purple', linewidth=2)\n",
    "plt.plot(lat, solar_mean_z-toa_sw_mean_z, label='Absorbed SW Radiation', color='orange', linewidth=2)\n",
    "plt.plot(lat, toa_lw_mean_z, label='Emitted LW Radiation', color='blue', linewidth=2)\n",
    "plt.plot(lat, toa_net_mean_z, label='Net Radiation', color='black', linewidth=2)\n",
    "plt.axvline(0, color='grey', linestyle=\"-\", linewidth=0.5)\n",
    "plt.axhline(0, color='grey', linestyle=\"-\", linewidth=0.5)\n",
    "plt.fill_between(lat, toa_net_mean_z, where=toa_net_mean_z > 0, color='lightcoral')\n",
    "plt.fill_between(lat, toa_net_mean_z, where=toa_net_mean_z < 0, color='lightcyan')\n",
    "plt.legend(loc='lower center', ncol=3)\n",
    "\n",
    "figname=dir_figs+'TOA_rad_budget_zonal_annual'\n",
    "plt.savefig(figname+'.png',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximation analytique du bilan radiatif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "<p>On approxime les variations méridiennes de la température par la fonction $T_m(lat) = a + b [cos^3(lat)(1+\\frac{3}{2}*sin^2(lat)]$</p>\n",
    "    \n",
    "<p>$a=a_1-\\gamma*z$</p>\n",
    "<p>$b=b_1*(1-\\frac{z}{z_T})$</p>\n",
    "<br>\n",
    "\n",
    "<p>$z$ : altitude</p>\n",
    "<p>$a_1 = –12 °C$</p>\n",
    "<p>$\\gamma = 3,14 °C.km^{–1}$</p>\n",
    "<p>$z_T = 11 km$ : altitude de la tropopause</p>\n",
    "<p>$b_1 = 40°C$ : différence de température équateur-pôle près de la surface</p>\n",
    "\n",
    "<br>\n",
    "<b>Coder cette fonction</b>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Profil méridien idéalisé de température."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def T(z):\n",
    "    return"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Tracé du profil méridien de température pour les altitudes 0 km, 1.5 km, 5.5 km, 10 km et 15 km."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,8))\n",
    "\n",
    "plt.xlabel('Latitude')\n",
    "plt.ylabel('Temperature (°C)')\n",
    "#plt.xticks(np.arange(-90, 120, 30), ('90S', '60S', '30S', 'Eq', '30N', '60N', '90N'))\n",
    "plt.xticks(np.arange(-90, 120, 30))\n",
    "plt.title('Idealized variation of annual-average temperature')\n",
    "plt.axvline(0, color='black', linestyle=\"--\")\n",
    "\n",
    "for alt in [0,1.5,5.5,10,15]:\n",
    "    plt.plot(lat, T(alt), label=str(alt)+'km')\n",
    "plt.legend()\n",
    "figname='./figs/ideal_temp'\n",
    "plt.savefig(figname+'.png',bbox_inches='tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "<p>On approxime les variations méridiennes du rayonnement incident par :</p>\n",
    "<p>$E_{sol}=E_0+E_1*cos(2*lat)$</p>\n",
    "<p>$E_0 = 298 W/m^2, E1 = 123 W/m^2$</p>\n",
    "<br>\n",
    "\n",
    "<p>On considère le rayonnement réfléchi constant avec la latitude :</p>\n",
    "<p>$E_2 = 110 W/m^2$</p>\n",
    "<p>Le rayonnement absorbé est donc :</p>\n",
    "<p>$E_{in} = E_{sol} - E_2$</p>\n",
    "<br>\n",
    "\n",
    "<p> On suppose que l'émission infrarouge est caractéristique d'une température idéalisée de moyenne troposphère $T_m$ définie précédemment (convertie en Kelvin) à l'altitude $z_m = 5,5 km$ et on approxime cette émission par la loi de Stefan-Boltzmann :</p>\n",
    "\n",
    "<p>$E_{out}=\\epsilon \\sigma T^4$</p>\n",
    "\n",
    "<p>$\\epsilon=0,9$</p>\n",
    "<p>$\\sigma=5,67*10^{-8}Wm^{-2}K^{-4}$</p>\n",
    "<br>\n",
    "    \n",
    "<b>Tracer les courbes théoriques de $E_{sol}$, $E_{out}$, $E_{in}$ et le bilan net $E_{in}-E_{out}$. Superposer ces courbes aux courbes des observations CERES. L'approximation analytique est-elle acceptable ?</b>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Bilan en moyenne annuelle et zonale, approximation algébrique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chauffage différentiel et transport méridien d'énergie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "<p>La circonférence d'un parallèle (cercle de latitude constante) est plus petite près des pôles qu'à l'équateur. Pour prendre en compte cet effet, on multiplie le bilan radiatif net algébrique par la circonférence d'un parallèle pour convertir le rayonnement (unité $W/m^2$) en $E_ϕ$ (unité $W/m$). $E_ϕ$ peut être interprété comme la puissance transférée vers (où depuis) une section de 1m qui entoure la Terre le long d'un parallèle</p>\n",
    "\n",
    "<br>\n",
    "    \n",
    "<b> Coder $E_ϕ$ et tracer sa distribution méridienne.</b>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Bilan net en moyenne annuelle et zonale avec prise en compte des distances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Transport méridien de chaleur."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "<p>La quantité $E_ϕ$ peut être définie comme un chauffage différentiel. Le transport méridien à chaque pôle est nul car le pôle est un point où convergent les méridiens. En sommant $E_ϕ$ pour toutes les ceintures de latitudes du pôle nord à n'importe quelle latitude on peut calculer le transport méridien total nécessaire pour que la circulation globale compense tous les déséquilibres méridiens au nord de cette latitude </p>\n",
    "\n",
    "<br>\n",
    "    \n",
    "<b>Utiliser la fonction scipy.integrate.cumtrapz pour implémenter le calcul du transport $Tr⁡(lat)=-\\sum_{lat=90°}^{-90°} E_ϕ \\Delta y$ et tracer sa distribution méridienne.</b>\n",
    "    \n",
    "<br>\n",
    "\n",
    "<b>La largeur de la ceinture de latitude $\\Delta y$ est relée à la variation de latitude $\\Delta lat$ par : $\\Delta y(km) = (111 km/°)*\\Delta lat (°)$</b>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "<p>Comparer avec le transport méridien réel ci-dessous :\n",
    "    K. E. Trenberth and J. M. Caron, 2001: “J. Climate”, 14, 3433-3443.</p>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"tr.png\" width=\"600\">"
   ]
  }
 ],
 "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}