{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MesoNHforC2OMODO\n",
    "### Notebook for Analyzing the MesoNHforC2OMODO Simulations of Deep Convection    \n",
    "\n",
    "*Jean-Pierre Chaboureau - University of Toulouse*   \n",
    "*CloudsforAfrica -  July 2026*\n",
    "\n",
    "This practical will guide participants through analysis of deep-convection simulations from the Meso-NH database for C2OMODO (available at https://doi.org/10.25326/789). The database contains multiple files, each around 50 Gb, from simulations of convective events at resolutions of 1-3 km or 200 m. The outputs include model variables such as 3D\n",
    "winds, cloud and precipitation data, as well as synthetic satellite outputs (e.g., radar reflectivity and brightness temperature in the infrared and microwave spectra). The goal is to analyze model output and pseudo-observation variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import xarray as xr\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import matplotlib.ticker as mticker\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "from matplotlib import colors\n",
    "from matplotlib.colors import BoundaryNorm\n",
    "import numpy as np\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file=\"MesoNH-ice3_DARWINF22_1km.nc\"\n",
    "file=\"MesoNH-ice3_CADDIWAF7_1km.nc\"\n",
    "file=\"MesoNH-ice3_CADDIWAF7_1km_8timesteps.nc\"\n",
    "data = xr.open_dataset(file)\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time=data.coords[\"time\"].values\n",
    "longitude=data.coords[\"longitude\"].values\n",
    "latitude=data.coords[\"latitude\"].values\n",
    "altitude=data.data_vars[\"ALT\"][:,:,:]\n",
    "level=data.coords[\"level\"].values\n",
    "lon1d=longitude[0,:]\n",
    "lat1d=latitude[:,0]\n",
    "ntime=len(time)\n",
    "\n",
    "def plot_axes(ax):\n",
    "    ax.coastlines()\n",
    "    gl = ax.gridlines(\n",
    "    draw_labels=True, linewidth=0.8,\n",
    "    color=\"gray\", alpha=0.5, linestyle=\"--\")\n",
    "# Set grid interval (every 1° here)\n",
    "    gl.xlocator = mticker.FixedLocator(range(-180, 181, 1))\n",
    "    gl.ylocator = mticker.FixedLocator(range(-90, 91, 1))\n",
    "# Format labels\n",
    "    gl.xformatter = LONGITUDE_FORMATTER\n",
    "    gl.yformatter = LATITUDE_FORMATTER\n",
    "# Style labels\n",
    "    gl.xlabel_style = {\"size\": 10}\n",
    "    gl.ylabel_style = {\"size\": 10}\n",
    "# Only label left and bottom axes\n",
    "    gl.top_labels = False\n",
    "    gl.right_labels = False\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animation of a 2D field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "\n",
    "yvar = 'IWP' # 'WICE' # 'aos_183TBT' #  'INPRR' # 'MSG2_IR108' # \n",
    "\n",
    "varname = data.data_vars[yvar].long_name + \\\n",
    "       ' (' + data.data_vars[yvar].units + ')'\n",
    "fig, ax = plt.subplots(figsize=(10, 10),\n",
    "    subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
    "cf = ax.pcolormesh(lon1d, lat1d, data.data_vars[yvar][0,:,:], #vmin=0, vmax=15,\n",
    "        shading=\"auto\", cmap='viridis', transform=ccrs.PlateCarree())\n",
    "cb = fig.colorbar(cf, label=varname,  orientation='vertical', shrink=0.5, pad=0.08)\n",
    "#\n",
    "plot_axes(ax)\n",
    "def update(frame):\n",
    "    cf.set_array(data.data_vars[yvar][frame, :, :].values[:-1, :-1].ravel())\n",
    "    return [cf]\n",
    "plt.close(fig)\n",
    "ani = FuncAnimation(fig, update, frames=ntime, interval=100, blit=False)\n",
    "HTML(ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time series of minimum/maximum values of a 2D field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yvar= 'IWP' # 'WICE' # 'INPRR' \n",
    "vardata = np.max(np.max(data.data_vars[yvar][:,:,:],axis=2),axis=1)\n",
    "#yvar= 'MSG2_IR108' # 'aos_183TBT' # \n",
    "#vardata = np.min(np.min(data.data_vars[yvar][:,:,:],axis=2),axis=1)\n",
    "varname = data.data_vars[yvar].long_name + \\\n",
    "       ' (' + data.data_vars[yvar].units + ')'\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "im = ax.plot(time,vardata, marker='x',linestyle='dotted')\n",
    "for it in range(ntime): ax.text(time[it],vardata[it],it,fontsize=12)\n",
    "ax.set_ylabel(varname)\n",
    "# plt.savefig('figTimeSeries.png', dpi=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select a time index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "it=2\n",
    "py_date =np.datetime64(time[it], 's').astype(datetime)\n",
    "print(f\"{py_date.strftime('%d-%m-%Y %H:%M:%S UTC')}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Horizontal cross-section of a 2D field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lvar2D=True\n",
    "if Lvar2D:\n",
    "    # for a 2D field\n",
    "    yvar = 'IWP' # 'aos_183TBT' # 'MSG2_IR108' # 'INPRR' #    \n",
    "    vardata  = data.data_vars[yvar][it,:,:]\n",
    "else:\n",
    "    # for a 3D field\n",
    "    ik=60 ; print('Altitude =',np.mean(altitude[ik,:,:].values),' m')\n",
    "    yvar = 'MRI' # cloud ice\n",
    "    vardata = data.data_vars[yvar][it,ik,:,:]\n",
    "#\n",
    "varname = data.data_vars[yvar].long_name + \\\n",
    "       ' (' + data.data_vars[yvar].units + ')'\n",
    "#\n",
    "fig, ax = plt.subplots(figsize=(10, 10),\n",
    "    subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
    "cmap='jet' # 'seismic' # 'viridis' 'Blues'\n",
    "cf = ax.pcolormesh(lon1d,lat1d,vardata,cmap=cmap, #vmax=50,\n",
    "                   transform=ccrs.PlateCarree())\n",
    "cb = fig.colorbar(cf, label=varname, orientation='vertical', shrink=0.5, pad=0.08)\n",
    "#\n",
    "ind_vardata = np.unravel_index(np.nanargmax(vardata, axis=None), vardata.shape)\n",
    "ax.text(lon1d[ind_vardata[1]],lat1d[ind_vardata[0]],'X')\n",
    "print('index of the location of max '+varname,ind_vardata)\n",
    "ax.plot([lon1d[0],lon1d[-1]], [lat1d[ind_vardata[0]],lat1d[ind_vardata[0]]],'-',color='black')\n",
    "#\n",
    "Ladd_contour=False\n",
    "if Ladd_contour:\n",
    "    ik=36 ; print('Altitude =',np.mean(altitude[ik,:,:].values),' m')\n",
    "    wwind= data.data_vars[\"WT\"][it,ik,:,:] # vertical velocity\n",
    "    ax.contour(lon1d,lat1d,wwind,colors='red',levels=[5])\n",
    "#ind_vardata = np.unravel_index(np.nanargmax(wwind, axis=None), vardata.shape)\n",
    "#print('index of the location of max '+yvar2D,ind_vardata)\n",
    "#\n",
    "Ladd_wind=False\n",
    "if Ladd_wind:\n",
    "    ik=36 ; print('Altitude =',np.mean(altitude[ik,:,:].values),' m')\n",
    "    uwind= data.data_vars[\"UT\"][it,ik,:,:]\n",
    "    vwind= data.data_vars[\"VT\"][it,ik,:,:]\n",
    "    ifreq=16\n",
    "    ax.quiver(lon1d[::ifreq],lat1d[::ifreq],\n",
    "          uwind[::ifreq,::ifreq],vwind[::ifreq,::ifreq],color='yellow')\n",
    "#\n",
    "py_date =np.datetime64(time[it], 's').astype(datetime)\n",
    "props = dict(boxstyle='round',facecolor='white', alpha=0.8)\n",
    "ax.text(0.05,0.92,f\"{py_date.strftime('%d-%m-%Y %H:%M:%S UTC')}\",\n",
    "        transform=ax.transAxes,fontsize=12,bbox=props)\n",
    "plot_axes(ax)\n",
    "#\n",
    "# zoom in [lon_min,lon_max,lat_min,lat_max]\n",
    "maplimits=[-20.5,-19.5,13.5,14.5] # [128,130,-14,-12.5]\n",
    "Lzoom=False\n",
    "if Lzoom: ax.set_extent(maplimits, crs=ccrs.PlateCarree())\n",
    "#\n",
    "#\n",
    "# plt.savefig('figHXS2Dfield.png', dpi=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Histogram of a 2D field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "vardata  = data.data_vars[yvar][it,:,:].values.flatten()\n",
    "varname = data.data_vars[yvar].long_name + \\\n",
    "       ' (' + data.data_vars[yvar].units + ')'\n",
    "ax.hist(vardata,bins=20,histtype='step') #,range=(0,20))\n",
    "ax.set_xlabel(varname,fontsize=16)\n",
    "ax.set_ylabel(\"# Samples\",fontsize=16)\n",
    "ax.set_yscale('log')\n",
    "# plt.savefig('figHistogram.png', dpi=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vertical section across a convective core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ij=ind_vardata[0]\n",
    "cloud= data.data_vars[\"MRC\"][it,:,ij,:] + data.data_vars[\"MRI\"][it,:,ij,:]\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "levels=[0.001,0.01,0.1,1,5,10]\n",
    "cmap=mpl.cm.get_cmap(name='hot_r')\n",
    "norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n",
    "#cf = ax.contourf(lon1d,level,cloud,cmap=cmap,norm=norm,levels=levels,extend='both')\n",
    "cf = ax.pcolormesh(lon1d,level,cloud,cmap=cmap,norm=norm)\n",
    "cb = fig.colorbar(cf, label='Cloud mixing ratio g kg-1', orientation='vertical', shrink=0.5, pad=0.02)\n",
    "#\n",
    "Ladd_temperature=False\n",
    "if Ladd_temperature:\n",
    "    temperature= data.data_vars[\"TEMP\"][it,:,ij,:]\n",
    "    ax.contour(lon1d,level,temperature,colors='red') #,levels=[1,10])\n",
    "# \n",
    "Ladd_wind=False\n",
    "if Ladd_wind:\n",
    "    uwind= data.data_vars[\"UT\"][it,:,ij,:]\n",
    "    wwind= 1.*data.data_vars[\"WT\"][it,:,ij,:]\n",
    "    ifreq=8\n",
    "    ax.quiver(lon1d[::ifreq],level[::ifreq],\n",
    "              uwind[::ifreq,::ifreq],wwind[::ifreq,::ifreq],color='blue')\n",
    "#\n",
    "ax.set_xlim(-20.5,-19.5)\n",
    "ax.set_xlabel('longitude (degree)',fontsize=16)\n",
    "ax.set_ylim(0,20000)\n",
    "ax.set_ylabel('altitude (m)',fontsize=16)\n",
    "#\n",
    "py_date =np.datetime64(time[it], 's').astype(datetime)\n",
    "props = dict(boxstyle='round',facecolor='white', alpha=0.8)\n",
    "ax.text(0.05,0.92,f\"{py_date.strftime('%d-%m-%Y %H:%M:%S UTC')}\",\n",
    "        transform=ax.transAxes,fontsize=12,bbox=props)\n",
    "# plt.savefig('figVXScloud.png', dpi=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D-histogram between IWP and BT at 183 GHz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "var1=data.data_vars[\"aos_183TBT\"][0,:,:].values.flatten()\n",
    "var2=data.data_vars[\"IWP\"][0,:,:].values.flatten()\n",
    "ax.hist2d(var1,var2, bins=40, norm=colors.LogNorm())\n",
    "ax.set_xlabel(\"aos_183TBT (K)\",fontsize=16)\n",
    "ax.set_ylabel(\"IWP (kg/m2)\",fontsize=16)\n",
    "# plt.savefig('figIWPvs183TBT.png', dpi=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select a second time index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "itp=it+1\n",
    "print('Difference in time:',time[itp]-time[it])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Horizontal cross-section of the 1-min change of a 2D field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yvar = 'IWP' # 'aos_183TBT' 'aos_325TBT' 'INPRR' 'MSG2_IR108' \n",
    "vardata = data.data_vars[yvar][itp,:,:]- data.data_vars[yvar][it,:,:]\n",
    "varname = 'Dif ' + data.data_vars[yvar].long_name + \\\n",
    "       ' (' + data.data_vars[yvar].units + ')'\n",
    "fig, ax = plt.subplots(figsize=(10, 10),\n",
    "    subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
    "cf = ax.pcolormesh(lon1d,lat1d,vardata,cmap='seismic',\n",
    "                   vmin=-7,vmax=7,transform=ccrs.PlateCarree())\n",
    "cb = fig.colorbar(cf, label=varname,  orientation='vertical', shrink=0.5, pad=0.08)\n",
    "wmax= np.nanmax(data.data_vars[\"WT\"][itp,:,:,:],axis=0)\n",
    "ax.plot([lon1d[0],lon1d[-1]], [lat1d[ind_vardata[0]],lat1d[ind_vardata[0]]],'-',color='black')\n",
    "#\n",
    "Ladd_contour=False\n",
    "if Ladd_contour:\n",
    "    yvarbis = 'aos_183TBT' \n",
    "    vardatabis = data.data_vars[yvarbis][itp,:,:]- data.data_vars[yvarbis][it,:,:]\n",
    "    ax.contour(lon1d,lat1d,vardatabis,colors='green',levels=[-5,5])\n",
    "#\n",
    "Lzoom=False\n",
    "if Lzoom: ax.set_extent(maplimits, crs=ccrs.PlateCarree())\n",
    "#\n",
    "py_date =np.datetime64(time[it], 's').astype(datetime)\n",
    "pyp_date =np.datetime64(time[itp], 's').astype(datetime)\n",
    "props = dict(boxstyle='round',facecolor='white', alpha=0.8)\n",
    "ax.text(0.05,0.89,\"Variable Difference\"\n",
    "        +\"\\nat \"+f\"{pyp_date.strftime('%H:%M:%S UTC')}\"+\" and\"\n",
    "        +\"\\nat \"+f\"{py_date.strftime('%H:%M:%S UTC')}\",\n",
    "        transform=ax.transAxes,fontsize=12,bbox=props)\n",
    "plot_axes(ax)\n",
    "# plt.savefig('figDelta183BT.png', dpi=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D-histogram between change in IWP and change in BT at 183 GHz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "var1=data.data_vars[\"aos_183TBT\"][itp,:,:].values.flatten()-data.data_vars[\"aos_183TBT\"][it,:,:].values.flatten()\n",
    "var2=data.data_vars[\"IWP\"][itp,:,:].values.flatten()-data.data_vars[\"IWP\"][it,:,:].values.flatten()\n",
    "ax.hist2d(var1,var2, bins=40, norm=colors.LogNorm())\n",
    "ax.set_xlabel(\"delta aos_183TBT (K)\",fontsize=16)\n",
    "ax.set_ylabel(\"delta IWP (kg/m2)\",fontsize=16)\n",
    "# plt.savefig('figdeltaIWPvsdelta183TBT.png', dpi=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vertical section across a convective core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ij=ind_vardata[0]\n",
    "varname='Cloud mixing ratio g kg-1'\n",
    "cloud = (data.data_vars[\"MRC\"][itp,:,ij,:] + data.data_vars[\"MRI\"][itp,:,ij,:] ) \\\n",
    "      - (data.data_vars[\"MRC\"][it,:,ij,:] + data.data_vars[\"MRI\"][it,:,ij,:])\n",
    "levels=[-0.1,-0.05,-0.01,0.01,0.05,0.1]\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "cmap='seismic'\n",
    "#cf = ax.contourf(lon1d,level,cloud,cmap=cmap,levels=levels,extend='both')\n",
    "cf = ax.pcolormesh(lon1d,level,cloud,cmap=cmap)\n",
    "cb = fig.colorbar(cf, label=varname, orientation='vertical', shrink=0.5, pad=0.02)\n",
    "#\n",
    "Ladd_temperature=False\n",
    "if Ladd_temperature:\n",
    "    temperature= data.data_vars[\"TEMP\"][it,:,ij,:]\n",
    "    ax.contour(longitude[ij,:],level,temperature,colors='red') #,levels=[1,10])\n",
    "    ax.set_xlabel('longitude (degree)',fontsize=16) ; ax.set_ylabel('altitude (m)',fontsize=16)\n",
    "#\n",
    "Ladd_wind=True\n",
    "if Ladd_wind:\n",
    "    uwind= data.data_vars[\"UT\"][it,:,ij,:]\n",
    "    wwind= 1.*data.data_vars[\"WT\"][it,:,ij,:]\n",
    "    ifreq=4\n",
    "    ax.quiver(longitude[ij,::ifreq],level[::ifreq],\n",
    "              uwind[::ifreq,::ifreq],wwind[::ifreq,::ifreq],color='black')\n",
    "#\n",
    "ax.set_xlim(-20.5,-19.5)\n",
    "ax.set_xlabel('longitude (degree)',fontsize=16)\n",
    "ax.set_ylim(0,20000)\n",
    "ax.set_ylabel('altitude (m)',fontsize=16)\n",
    "#\n",
    "py_date =np.datetime64(time[it], 's').astype(datetime)\n",
    "props = dict(boxstyle='round',facecolor='white', alpha=0.8)\n",
    "ax.text(0.05,0.92,f\"{py_date.strftime('%d-%m-%Y %H:%M:%S UTC')}\",\n",
    "        transform=ax.transAxes,fontsize=12,bbox=props)\n",
    "# plt.savefig('figVXScloud.png', dpi=100)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Format de la Cellule Texte Brut",
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}