{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "93470e7b",
   "metadata": {},
   "source": [
    "# Volcanic Forcing on AMOC Transitions - Figure 3\n",
    "\n",
    "This notebook generates a comprehensive multi-panel figure analyzing the impact of volcanic forcing on Atlantic Meridional Overturning Circulation (AMOC) transitions in a complex climate model. The analysis examines how volcanic events affect AMOC stability at different climate states.\n",
    "\n",
    "## Figure Structure\n",
    "\n",
    "**LEFT SIDE (Interstadial Age Test):**\n",
    "- **Panel A**: Overview showing AMOC ensemble with three zoom regions at different interstadial ages\n",
    "- **Panels B, C, D**: Detailed zoom views showing 50-member ensembles at three different interstadial ages:\n",
    "  - Panel B: Younger interstadial state\n",
    "  - Panel C: Mid-age interstadial state  \n",
    "  - Panel D: Older interstadial state (closer to stadial transition)\n",
    "\n",
    "**RIGHT SIDE (Stadial Age Test):**\n",
    "- **Panel E**: Overview showing AMOC ensemble with three zoom regions at different stadial ages\n",
    "- **Panels F, G, H**: Detailed zoom views showing 50-member ensembles at three different stadial ages:\n",
    "  - Panel F: Earlier in stadial period\n",
    "  - Panel G: Mid-stadial period\n",
    "  - Panel H: Later in stadial period (closer to interstadial transition)\n",
    "\n",
    "## Key Features in Each Zoom Panel\n",
    "\n",
    "Each zoom panel (B-D and F-H) includes:\n",
    "- **50 ensemble members** comparing:\n",
    "  - Blue lines: Noise-only ensemble (no volcanic forcing)\n",
    "  - Red lines: Volcanic ensemble (with volcanic forcing at specified time)\n",
    "- **Histogram inset** showing final AMOC state distribution\n",
    "- **Transition probabilities** calculated as percentage of ensemble members transitioning below basin boundary threshold\n",
    "- **Basin boundary** marked as green shaded region (critical threshold for AMOC state transitions)\n",
    "- **Volcanic event timing** indicated by red vertical line and annotation\n",
    "\n",
    "## Scientific Context\n",
    "\n",
    "The analysis tests the hypothesis that volcanic forcing can trigger AMOC transitions, and that the effectiveness depends on:\n",
    "1. **Climate state age**: Proximity to natural transition threshold\n",
    "2. **Background state**: Whether system is in interstadial (warm) or stadial (cold) state\n",
    "\n",
    "Results show volcanic forcing is most effective at triggering transitions when the climate system is already close to a natural tipping point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "3bef220a-eb8e-4182-9e52-b844de098fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required libraries for data processing, visualization, and climate data analysis\n",
    "# - cftime: Handle non-standard calendars used in climate models\n",
    "# - xarray: Work with multi-dimensional labeled arrays (NetCDF data)\n",
    "# - matplotlib/seaborn: Create publication-quality figures\n",
    "# - pandas/numpy: Data manipulation and numerical operations\n",
    "\n",
    "import datetime\n",
    "import datetime as dt\n",
    "import glob\n",
    "\n",
    "import cftime\n",
    "import matplotlib.dates as mdates\n",
    "import matplotlib.gridspec as gridspec\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import seaborn as sns\n",
    "from typing import Union, List\n",
    "\n",
    "import xarray as xr\n",
    "from xclim import ensembles\n",
    "\n",
    "from matplotlib.patches import Rectangle#from scipy.integrate import trapz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "6a7d2152-3a50-4b41-9eac-aa508edd2e80",
   "metadata": {},
   "outputs": [],
   "source": [
    "def setaxesfsize(axl, fontsize):\n",
    "    \"\"\"\n",
    "    Set font size for tick labels on x and y axes of a given matplotlib axis object.\n",
    "\n",
    "    Parameters:\n",
    "    axl (matplotlib.axes.Axes): The axis object to modify.\n",
    "    fontsize (int): The desired font size for the tick labels.\n",
    "\n",
    "    Returns:\n",
    "    None\n",
    "    \"\"\"\n",
    "    ticklabelbot = axl.get_xticklabels()\n",
    "    ticklabelleft = axl.get_yticklabels()\n",
    "\n",
    "    for labelx in ticklabelbot:\n",
    "        labelx.set_fontsize(fontsize)\n",
    "\n",
    "    for labely in ticklabelleft:\n",
    "        labely.set_fontsize(fontsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "acfbe862-e58d-4f43-9b84-0dc399499d6f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.coding.times.CFDatetimeCoder at 0x779b45aa7560>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Configure time decoder for climate model data\n",
    "# Creates a CFDatetimeCoder instance to properly decode time coordinates from NetCDF files\n",
    "# using cftime (handles non-standard calendars like 'noleap' used in climate models)\n",
    "\n",
    "time_coder = xr.coders.CFDatetimeCoder(use_cftime=True)\n",
    "time_coder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "6e6b6c01-003c-4e25-a45b-313326fde1f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_dataset(filename: str) -> Union[xr.DataArray, None]:\n",
    "    \"\"\"\n",
    "    Load a dataset from a NetCDF file.\n",
    "\n",
    "    Args:\n",
    "        filename (str): The name of the NetCDF file to load.\n",
    "\n",
    "    Returns:\n",
    "        xr.DataArray: The loaded and processed dataset.\n",
    "        None: If there was an error loading the file.\n",
    "\n",
    "    Raises:\n",
    "        FileNotFoundError: If the specified file does not exist.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        dataset = xr.load_dataset(filename, decode_times=time_coder)\n",
    "        return dataset.to_dataarray().rename('MOC').drop_vars('variable').squeeze()\n",
    "    except FileNotFoundError:\n",
    "        print(f\"Error: File '{filename}' not found.\")\n",
    "        return None\n",
    "    except Exception as e:\n",
    "        print(f\"Error loading dataset: {str(e)}\")\n",
    "        return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "92dfa2e0-8d1b-47a7-985c-15717eb38e3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_mean_last10yr(data: xr.DataArray, nruns: int, t_from: int = -72, t_to: int = -60, endmon: int = -1) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    Calculate the mean of the last 10 years for each ensemble member.\n",
    "\n",
    "    Args:\n",
    "        data (xr.DataArray): The input data array.\n",
    "        nruns (int): The number of ensemble runs.\n",
    "        t_from (int, optional): Start index for time slice. Defaults to -72.\n",
    "        t_to (int, optional): End index for time slice. Defaults to -60.\n",
    "        endmon (int, optional): End month for data selection. Defaults to -1.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: Array of mean values for the last 10 years of each ensemble member.\n",
    "\n",
    "    Raises:\n",
    "        ValueError: If nruns is less than or equal to 0.\n",
    "    \"\"\"\n",
    "    if nruns <= 0:\n",
    "        raise ValueError(\"Number of runs (nruns) must be greater than 0.\")\n",
    "\n",
    "    try:\n",
    "        return np.array([\n",
    "            data[i][0:endmon].isel(time=slice(t_from, t_to)).mean(dim=\"time\").values\n",
    "            for i in range(nruns)\n",
    "        ])\n",
    "    except Exception as e:\n",
    "        print(f\"Error calculating mean: {str(e)}\")\n",
    "        return np.array([])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "1393e09c-b75f-4b4b-8cc2-709b1a8e67f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load ensemble data from NetCDF files\n",
    "# This cell loads AMOC maximum strength data (10-year running mean) for multiple experiments:\n",
    "# - ctrl230i: Control run (no volcanic forcing)\n",
    "# - ens_eq_no*: Noise-only ensembles at different model years (3701, 6701, 7701, 8501, 8701, 8901)\n",
    "# - ens_eq_*: Volcanic forcing ensembles at corresponding years\n",
    "# Each ensemble contains 50 members comparing volcanic vs. noise-only scenarios\n",
    "\n",
    "root_dir = '/home/guido/Dropbox/ipynb/D-O_analysis/volcanic/'\n",
    "\n",
    "# Load each variable separately from its corresponding NetCDF file\n",
    "ctrl230i_moc_max_rm10a = load_dataset(root_dir + 'ctrl230i_moc_max_rm10a.nc')\n",
    "\n",
    "ens_eq_no3701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_no3701_moc_max_rm10a.nc')\n",
    "ensmean_eq_no3701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_no3701_moc_max_rm10a.nc')\n",
    "ens_eq_no6701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_no6701_moc_max_rm10a.nc')\n",
    "ensmean_eq_no6701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_no6701_moc_max_rm10a.nc')\n",
    "ens_eq_no7701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_no7701_moc_max_rm10a.nc')\n",
    "ensmean_eq_no7701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_no7701_moc_max_rm10a.nc')\n",
    "\n",
    "ens_eq_3701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_3701_moc_max_rm10a.nc')\n",
    "ensmean_eq_3701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_3701_moc_max_rm10a.nc')\n",
    "ens_eq_6701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_6701_moc_max_rm10a.nc')\n",
    "ensmean_eq_6701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_6701_moc_max_rm10a.nc')\n",
    "ens_eq_7701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_7701_moc_max_rm10a.nc')\n",
    "ensmean_eq_7701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_7701_moc_max_rm10a.nc')\n",
    "\n",
    "ens_eq_no8501_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_no8501_moc_max_rm10a.nc')\n",
    "ensmean_eq_no8501_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_no8501_moc_max_rm10a.nc')\n",
    "ens_eq_no8701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_no8701_moc_max_rm10a.nc')\n",
    "ensmean_eq_no8701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_no8701_moc_max_rm10a.nc')\n",
    "ens_eq_no8901_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_no8901_moc_max_rm10a.nc')\n",
    "ensmean_eq_no8901_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_no8901_moc_max_rm10a.nc')\n",
    "\n",
    "ensmean_eq_8901_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_8901_moc_max_rm10a.nc')\n",
    "\n",
    "ens_eq_8501_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_8501_moc_max_rm10a.nc')\n",
    "ens_eq_8901_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_8901_moc_max_rm10a.nc')\n",
    "\n",
    "ensmean_eq_8501_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_8501_moc_max_rm10a.nc')\n",
    "ensmean_eq_8701_moc_max_rm10a = load_dataset(root_dir + 'ensmean_eq_8701_moc_max_rm10a.nc')\n",
    "ens_eq_8701_moc_max_rm10a = load_dataset(root_dir + 'ens_eq_8701_moc_max_rm10a.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ba217332-068a-4fd2-a0a5-45e68b04f9a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define number of ensemble members for each experiment\n",
    "# All experiments use 50-member ensembles to ensure robust statistics\n",
    "# for transition probability calculations\n",
    "\n",
    "nruns_3701 = 50\n",
    "nruns_6701 = 50\n",
    "nruns_7701 = 50\n",
    "nruns_8901 = 50\n",
    "nruns_8501 = 50\n",
    "nruns_8701 = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "d057a7d5-8c2a-4b03-ad3f-ea679474732c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate final AMOC state for each ensemble member\n",
    "# Computes the mean AMOC strength over the last 10 years of each simulation\n",
    "# This is used to determine if a transition occurred (below basin boundary threshold)\n",
    "# and to create histogram distributions showing final states\n",
    "\n",
    "ens_eq_no3701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_no3701_moc_max_rm10a, nruns_3701)\n",
    "ens_eq_3701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_3701_moc_max_rm10a, nruns_3701)\n",
    "\n",
    "ens_eq_no6701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_no6701_moc_max_rm10a, nruns_6701)\n",
    "ens_eq_6701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_6701_moc_max_rm10a, nruns_6701)\n",
    "\n",
    "ens_eq_no7701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_no7701_moc_max_rm10a, nruns_7701)\n",
    "ens_eq_7701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_7701_moc_max_rm10a, nruns_7701)\n",
    "\n",
    "ens_eq_no8501_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_no8501_moc_max_rm10a, nruns_8501)\n",
    "ens_eq_8501_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_8501_moc_max_rm10a, nruns_8501)\n",
    "\n",
    "ens_eq_no8701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_no8701_moc_max_rm10a, nruns_8701)\n",
    "\n",
    "ens_eq_8701_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_8701_moc_max_rm10a, nruns_8701)\n",
    "\n",
    "ens_eq_8901_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_8901_moc_max_rm10a, nruns_8901)\n",
    "ens_eq_no8901_moc_max_mean_last10yr = calculate_mean_last10yr(ens_eq_no8901_moc_max_rm10a, nruns_8901)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "b7ed5a56-549f-4e51-bf38-1b7bd8a7970e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 1477180 characters]",
      "text/plain": [
       "<Figure size 1600x1600 with 14 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create main figure with multi-panel layout\n",
    "# This cell sets up the complete Figure 3 structure using matplotlib gridspec:\n",
    "# - 140×200 grid allows precise control of subplot positioning\n",
    "# - Left side (panels A-D): Interstadial age test at years 3701, 6701, 7701\n",
    "# - Right side (panels E-H): Stadial age test at years 8501, 8701, 8901\n",
    "# - Each side has overview plot + 3 zoom panels with histogram insets\n",
    "\n",
    "#### Equatorial\n",
    "fig = plt.figure(figsize=(16, 16))\n",
    "\n",
    "# set the style\n",
    "sns.set_style(\"ticks\")\n",
    "sns.set_context(\"paper\")\n",
    "\n",
    "zoomcolor = (\n",
    "    0.8,\n",
    "    0.8,\n",
    "    0.8,\n",
    "    0.3,\n",
    ")  # lightgray with 0.3 transparency, we are going to zoom out this area\n",
    "\n",
    "# Create a 1x6 grid (this can be adjusted based on your needs)\n",
    "gs = gridspec.GridSpec(140, 200)\n",
    "\n",
    "# Three horizontal axes spanning the entire width\n",
    "zoomcolor = (0.8, 0.8, 0.8, 0.3)  # lightgray with 0.3 transparency, we are going to zoom out this area\n",
    "tlmid = 97\n",
    "trmid = 103\n",
    "ax0 = fig.add_subplot(gs[0:20, :tlmid])\n",
    "ax1 = fig.add_subplot(gs[26:60, 0:70], facecolor=zoomcolor)\n",
    "axins1 = fig.add_subplot(gs[26:60, 70:tlmid], facecolor=zoomcolor)\n",
    "ax2 = fig.add_subplot(gs[66:100, 0:70], facecolor=zoomcolor)\n",
    "axins2 = fig.add_subplot(gs[66:100, 70:tlmid], facecolor=zoomcolor)\n",
    "ax3 = fig.add_subplot(gs[106:140, 0:70], facecolor=zoomcolor)\n",
    "axins3 = fig.add_subplot(gs[106:140, 70:tlmid], facecolor=zoomcolor)\n",
    "\n",
    "# Duplicate axes on the right\n",
    "ax0_dup = fig.add_subplot(gs[0:20, trmid:])\n",
    "ax1_dup = fig.add_subplot(gs[26:60, trmid:170], facecolor=zoomcolor)\n",
    "axins1_dup = fig.add_subplot(gs[26:60, 170:197], facecolor=zoomcolor)\n",
    "ax2_dup = fig.add_subplot(gs[66:100, trmid:170], facecolor=zoomcolor)\n",
    "axins2_dup = fig.add_subplot(gs[66:100, 170:197], facecolor=zoomcolor)\n",
    "ax3_dup = fig.add_subplot(gs[106:140, trmid:170], facecolor=zoomcolor)\n",
    "axins3_dup = fig.add_subplot(gs[106:140, 170:197], facecolor=zoomcolor)\n",
    "\n",
    "ax = [ax0, ax1, axins1, ax2, axins2, ax3, axins3, ax0_dup, ax1_dup, axins1_dup, ax2_dup, axins2_dup, ax3_dup, axins3_dup]\n",
    "\n",
    "##############################\n",
    "\n",
    "\n",
    "##### Top Left  ##################################################################################\n",
    "################## Main Plot\n",
    "lw = 1.5\n",
    "lalp = 0.7\n",
    "ax0ymax = 32\n",
    "ax0ymin = 0\n",
    "axlbls = 14  # axes label size\n",
    "axts = 12  # axes tick label size\n",
    "# v2601 = cftime.num2date(2601*365, 'days since 0000-01-01 00:00:00', calendar='noleap')\n",
    "# ax0.vlines(v2601,ax0ymin,ax0ymax, color='red', alpha=0.7, linewidths=2.0, zorder=1)\n",
    "#vgrid0 = [yr * 365 for yr in range(2701, 10000, 1000)]\n",
    "#[\n",
    "#    ax0.vlines(\n",
    "#        vgrid0[i], ax0ymin, ax0ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1\n",
    "#    )\n",
    "#    for i in [1,4,5]\n",
    "#]\n",
    "vgrid0 = [\n",
    "    cftime.num2date(yr * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\")\n",
    "    for yr in [3701, 6701, 7901]\n",
    "]\n",
    "\n",
    "[\n",
    "    ax0.vlines(\n",
    "        mdates.date2num(vgrid0[i]), ax0ymin, ax0ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1\n",
    "    )\n",
    "    for i in [0,1,2]\n",
    "]\n",
    "\n",
    "# put in boxes for warm to cold and cold to warm transition probability study\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "\n",
    "################ these are the insets that we will zoom in on\n",
    "xzoom1_start = 3601\n",
    "xzoom1_end = 4001\n",
    "xzoom2_start = 6601\n",
    "xzoom2_end = 7001\n",
    "xzoom3_start = 7601\n",
    "xzoom3_end = 8001\n",
    "##############################\n",
    "\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "boxdates1_start = mdates.date2num(cftime.num2date(xzoom1_start * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "boxdates1_end = mdates.date2num(cftime.num2date(xzoom1_end * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "\n",
    "# add the grey rectangles of interest for zoom plots\n",
    "rect_width = boxdates1_end - boxdates1_start\n",
    "rectxmin = 7\n",
    "rectxmax = 29\n",
    "rect_height = rectxmax - rectxmin\n",
    "rect = Rectangle(\n",
    "    (boxdates1_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    color=zoomcolor,\n",
    ")\n",
    "border = Rectangle(\n",
    "    (boxdates1_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    facecolor=\"none\",\n",
    "    edgecolor=\"black\",\n",
    "    linewidth=0.7,\n",
    ")\n",
    "ax0.add_patch(rect)\n",
    "ax0.add_patch(border)\n",
    "\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "boxdates2_start = mdates.date2num(cftime.num2date(xzoom2_start * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "boxdates2_end = mdates.date2num(cftime.num2date(xzoom2_end * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "rect_width = boxdates2_end - boxdates2_start\n",
    "rectxmin = 7\n",
    "rectxmax = 29\n",
    "rect_height = rectxmax - rectxmin\n",
    "rect = Rectangle(\n",
    "    (boxdates2_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    color=zoomcolor,\n",
    ")\n",
    "border = Rectangle(\n",
    "    (boxdates2_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    facecolor=\"none\",\n",
    "    edgecolor=\"black\",\n",
    "    linewidth=0.7,\n",
    ")\n",
    "ax0.add_patch(rect)\n",
    "ax0.add_patch(border)\n",
    "\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "boxdates3_start = mdates.date2num(cftime.num2date(xzoom3_start * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "boxdates3_end = mdates.date2num(cftime.num2date(xzoom3_end * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "rect_width = boxdates3_end - boxdates3_start\n",
    "rectxmin = 7\n",
    "rectxmax = 29\n",
    "rect_height = rectxmax - rectxmin\n",
    "rect = Rectangle(\n",
    "    (boxdates3_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    color=zoomcolor,\n",
    ")\n",
    "border = Rectangle(\n",
    "    (boxdates3_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    facecolor=\"none\",\n",
    "    edgecolor=\"black\",\n",
    "    linewidth=0.7,\n",
    ")\n",
    "ax0.add_patch(rect)\n",
    "ax0.add_patch(border)\n",
    "\n",
    "\n",
    "cntcol = \"black\"\n",
    "cntalp = 0.5\n",
    "# volcanic run colors\n",
    "vcols = [\"purple\", \"blue\", \"red\"]\n",
    "# eruption colors\n",
    "ecols = [\"magenta\", \"cyan\", \"orange\"]\n",
    "\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax0,\n",
    "    color=\"blue\",\n",
    "    linewidth=1.0,\n",
    "    alpha=0.7,\n",
    "    label=\"Control\",\n",
    ")\n",
    "ax0.axhline(\n",
    "    20,\n",
    "    0,\n",
    "    1,\n",
    "    color=\"k\",\n",
    "    alpha=0.5,\n",
    "    linewidth=1.0,\n",
    "    linestyle=\"dashed\",\n",
    "    label=\"Upper Saddle Node\",\n",
    ")\n",
    "ax0.axhline(\n",
    "    15.5,\n",
    "    0,\n",
    "    1,\n",
    "    color=\"k\",\n",
    "    alpha=0.5,\n",
    "    linewidth=1.0,\n",
    "    linestyle=\"dotted\",\n",
    "    label=\"Lower Saddle Node\",\n",
    ")\n",
    "\n",
    "ax0.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "ax0.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "# ax0.grid()\n",
    "# ax0.set_title(\"Volcanic Test\", size=16, loc='center')\n",
    "ax0.set_title(\"\", size=16, loc=\"center\")\n",
    "ax0.set_ylim(ax0ymin, ax0ymax)\n",
    "setaxesfsize(ax0, axts)\n",
    "ax0.legend(loc=3, fontsize=10, ncol=1)\n",
    "\n",
    "##### Top Right  ##################################################################################\n",
    "\n",
    "xinmin = 0\n",
    "xinmax = 18\n",
    "#ealp = 0.8 # transparency of ensemble runs\n",
    "################## Main Plot\n",
    "lw = 1.5\n",
    "lalp = 0.7\n",
    "ax0ymax = 32\n",
    "ax0ymin = 0\n",
    "axlbls = 14  # axes label size\n",
    "axts = 12  # axes tick label size\n",
    "# v2601 = cftime.num2date(2601*365, 'days since 0000-01-01 00:00:00', calendar='noleap')\n",
    "# ax0.vlines(v2601,ax0ymin,ax0ymax, color='red', alpha=0.7, linewidths=2.0, zorder=1)\n",
    "vgrid0_dup = [\n",
    "    cftime.num2date(yr * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\")\n",
    "    for yr in [8501, 8701, 8901]\n",
    "]\n",
    "\n",
    "[\n",
    "    ax0_dup.vlines(\n",
    "        mdates.date2num(vgrid0_dup[i]), ax0ymin, ax0ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1\n",
    "    )\n",
    "    for i in [0,1,2]\n",
    "]\n",
    "\n",
    "\n",
    "# put in boxes for warm to cold and cold to warm transition probability study\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "\n",
    "################ these are the insets that we will zoom in on\n",
    "xzoom1_dup_start = 8401\n",
    "xzoom1_dup_end = 8801\n",
    "xzoom2_dup_start = 8601\n",
    "xzoom2_dup_end = 9001\n",
    "xzoom3_dup_start = 8801\n",
    "xzoom3_dup_end = 9201\n",
    "##############################\n",
    "endmon = -1\n",
    "##############################\n",
    "\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "boxdates1_dup_start = mdates.date2num(cftime.num2date(xzoom1_dup_start * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "boxdates1_dup_end = mdates.date2num(cftime.num2date(xzoom1_dup_end * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "rect_width = boxdates1_dup_end - boxdates1_dup_start\n",
    "rectxmin = 8\n",
    "rectxmax = 26\n",
    "rect_height = rectxmax - rectxmin\n",
    "rect = Rectangle(\n",
    "    (boxdates1_dup_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    color=zoomcolor,\n",
    ")\n",
    "border = Rectangle(\n",
    "    (boxdates1_dup_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    facecolor=\"none\",\n",
    "    edgecolor=\"black\",\n",
    "    linewidth=0.7,\n",
    ")\n",
    "ax0_dup.add_patch(rect)\n",
    "ax0_dup.add_patch(border)\n",
    "\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "boxdates2_dup_start = mdates.date2num(cftime.num2date(xzoom2_dup_start * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "boxdates2_dup_end = mdates.date2num(cftime.num2date(xzoom2_dup_end * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "rect_width = boxdates2_dup_end - boxdates2_dup_start\n",
    "rectxmin = 10\n",
    "rectxmax = 28\n",
    "rect_height = rectxmax - rectxmin\n",
    "rect = Rectangle(\n",
    "    (boxdates2_dup_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    color=zoomcolor,\n",
    ")\n",
    "border = Rectangle(\n",
    "    (boxdates2_dup_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    facecolor=\"none\",\n",
    "    edgecolor=\"black\",\n",
    "    linewidth=0.7,\n",
    ")\n",
    "ax0_dup.add_patch(rect)\n",
    "ax0_dup.add_patch(border)\n",
    "\n",
    "# Create a rectangle starting from the 1st date to the 2nd date\n",
    "boxdates3_dup_start = mdates.date2num(cftime.num2date(xzoom3_dup_start * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "boxdates3_dup_end = mdates.date2num(cftime.num2date(xzoom3_dup_end * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"))\n",
    "rect_width = boxdates3_dup_end - boxdates3_dup_start\n",
    "rectxmin = 12\n",
    "rectxmax = 30\n",
    "rect_height = rectxmax - rectxmin\n",
    "rect = Rectangle(\n",
    "    (boxdates3_dup_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    color=zoomcolor,\n",
    ")\n",
    "border = Rectangle(\n",
    "    (boxdates3_dup_start, rectxmin),\n",
    "    width=rect_width,\n",
    "    height=rect_height,\n",
    "    facecolor=\"none\",\n",
    "    edgecolor=\"black\",\n",
    "    linewidth=0.7,\n",
    ")\n",
    "ax0_dup.add_patch(rect)\n",
    "ax0_dup.add_patch(border)\n",
    "\n",
    "\n",
    "cntcol = \"black\"\n",
    "cntalp = 0.5\n",
    "# volcanic run colors\n",
    "vcols = [\"purple\", \"blue\", \"red\"]\n",
    "# eruption colors\n",
    "ecols = [\"magenta\", \"cyan\", \"orange\"]\n",
    "\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax0_dup,\n",
    "    color=\"blue\",\n",
    "    linewidth=1.0,\n",
    "    alpha=0.7,\n",
    "    label=\"Control\",\n",
    ")\n",
    "#zoom in to later half of control simulation\n",
    "ax0_dup.set_xlim(cftime.num2date(\n",
    "        7500 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    "    ),\n",
    "                 cftime.num2date(\n",
    "        10000 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    "    ))\n",
    "ax0_dup.axhline(\n",
    "    20,\n",
    "    0,\n",
    "    1,\n",
    "    color=\"k\",\n",
    "    alpha=0.5,\n",
    "    linewidth=1.0,\n",
    "    linestyle=\"dashed\",\n",
    "    label=\"Upper Saddle Node\",\n",
    ")\n",
    "ax0_dup.axhline(\n",
    "    15.5,\n",
    "    0,\n",
    "    1,\n",
    "    color=\"k\",\n",
    "    alpha=0.5,\n",
    "    linewidth=1.0,\n",
    "    linestyle=\"dotted\",\n",
    "    label=\"Lower Saddle Node\",\n",
    ")\n",
    "\n",
    "#ax0_dup.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "#ax0_dup.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "ax0_dup.set_ylabel(\"\", size=axlbls)\n",
    "ax0_dup.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "# ax0_dup.grid()\n",
    "# ax0_dup.set_title(\"Volcanic Test\", size=16, loc='center')\n",
    "ax0_dup.set_title(\"\", size=16, loc=\"center\")\n",
    "ax0_dup.set_ylim(ax0ymin, ax0ymax)\n",
    "setaxesfsize(ax0_dup, axts)\n",
    "ax0_dup.legend(loc=3, fontsize=10, ncol=1)\n",
    "\n",
    "# But turn off the y-axis labels and ticks\n",
    "#ax0_dup.set_yticklabels([])\n",
    "#ax0_dup.set_yticks([])\n",
    "#ax0_dup.set_ylabel(\"\")\n",
    "\n",
    "\n",
    "################################## Zoom Plots ##################################################\n",
    "# for position of \"Volcanic Event\" text label \n",
    "vitfont = 12\n",
    "nitfont = 12\n",
    "\n",
    "\n",
    "vlabxy = (0.254, 0.54)\n",
    "vlabxyt = (0.02, 0.54)\n",
    "\n",
    "######################################################\n",
    "################## Zoom Plot1 ########################\n",
    "sns.set_style(\"darkgrid\")\n",
    "\n",
    "esalp = 0.6  # ensemble line transparency\n",
    "esmalp = 0.9  # ensemble mean line tranparency\n",
    "eslw = 0.5  # ensemble line width\n",
    "esmlw = 2.0  # ensemble mean line width\n",
    "origlw = 2.0  # original runs line width\n",
    "\n",
    "# original lines\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax1,\n",
    "    color=\"blue\",\n",
    "    linewidth=origlw,\n",
    "    alpha=esalp,\n",
    "    label=\"Control\",\n",
    ")\n",
    "#eq_v1000_moc_max_rm10a.plot(\n",
    "#    ax=ax1,\n",
    "#    color=\"red\",\n",
    "#    linewidth=origlw,\n",
    "#    alpha=esalp,\n",
    "#    label=\"CO$_2$ = 230ppm (Periodic EQ 500Tg SO$_4$)\",\n",
    "#)\n",
    "# esemble member zero , for figure label\n",
    "ens_eq_no3701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax1,\n",
    "    color=\"blue\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Noise Ensemble\",\n",
    ")\n",
    "ens_eq_3701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax1,\n",
    "    color=\"red\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Volcanic Ensemble\"\n",
    ")\n",
    "# other esemble members\n",
    "[\n",
    "    ens_eq_no3701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax1, color=\"blue\", linewidth=eslw, alpha=esalp\n",
    "    )\n",
    "    for i in range(1, nruns_3701)\n",
    "]\n",
    "[\n",
    "    ens_eq_3701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax1, color=\"red\", linewidth=eslw, alpha=esalp)\n",
    "    for i in range(1, nruns_3701)\n",
    "]\n",
    "\n",
    "#ensmean_eq_no3701_moc_max_rm10a.plot(\n",
    "#    ax=ax1,\n",
    "#    color=\"cyan\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Noise (Ens Mean)\",\n",
    "#)\n",
    "#ensmean_eq_3701_moc_max_rm10a.plot(\n",
    "#    ax=ax1,\n",
    "#    color=\"yellow\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Volc (Ens Mean)\",\n",
    "#)\n",
    "\n",
    "# Volcano at 3701 line\n",
    "eq_v3701 = cftime.num2date(\n",
    "    3701 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ymin = 8\n",
    "ymax = 26\n",
    "ax1.vlines(eq_v3701, ymin, ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1)\n",
    "setaxesfsize(ax1, axts)\n",
    "ax1.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "ax1.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "ax1.set_ylim(ymin, ymax)\n",
    "#xmin = cftime.num2date(\n",
    "#    xzoom1_start * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    "#)\n",
    "#xmax = cftime.num2date(\n",
    "#    xzoom1_end * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    "#)\n",
    "# need to offset so it ends at 6000\n",
    "xminval = 3600\n",
    "xmaxval = 4000\n",
    "xmin = cftime.num2date(\n",
    "    xminval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "xmax = cftime.num2date(\n",
    "    xmaxval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ax1.set_xlim(xmin, xmax)\n",
    "\n",
    "# Fix the ticks and tick labels to make it look like inset is a separate plot\n",
    "# Fetch current x-tick labels from outer plot\n",
    "xticks = ax1.get_xticks()\n",
    "labels = [item.get_text() for item in ax1.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    # labels[0] = ''\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "ax1.set_xticks(xticks)\n",
    "ax1.set_xticklabels(labels)\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "ax1.grid()\n",
    "\n",
    "vefc='#ffffff'\n",
    "vealp=1.0\n",
    "ax1.legend(loc=3, fontsize=10, facecolor=vefc, framealpha=vealp)\n",
    "ax1.annotate('Volcanic\\nEvent',\n",
    "             xy=vlabxy, xycoords='axes fraction',\n",
    "             xytext=vlabxyt, textcoords='axes fraction',\n",
    "             arrowprops=dict(facecolor='red', shrink=0.05),\n",
    "             horizontalalignment='left', verticalalignment='top',\n",
    "             size=12,\n",
    "             bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))\n",
    "\n",
    "ax1.tick_params(bottom=False, right=False)\n",
    "\n",
    "ax1.set_title(\"\")\n",
    "\n",
    "\n",
    "####### INSET\n",
    "# Create the inset axes, choosing a position and size manually\n",
    "# The parameters are [left, bottom, width, height] in the figure’s coordinate system\n",
    "# inset axes....\n",
    "# change this to its own gridpsect now\n",
    "#axins1 = ax1.inset_axes([0.75, 0.0, 0.25, 1.0])\n",
    "x1, x2, y1, y2 = 0, 20, ymin, ymax\n",
    "axins1.set_xlim(x1, x2)\n",
    "axins1.set_ylim(y1, y2)\n",
    "\n",
    "bwidth = 0.5\n",
    "sns.histplot(\n",
    "    ax=axins1,\n",
    "    y=ens_eq_no3701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"blue\",\n",
    "    alpha=0.3,\n",
    "    label=\"No Volc\",\n",
    ")\n",
    "# sm.ProbPlot(beforedata[0][0]).qqplot(line='s', ax=ax[1])\n",
    "sns.histplot(\n",
    "    ax=axins1,\n",
    "    y=ens_eq_3701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"red\",\n",
    "    alpha=0.3,\n",
    "    label=\"Volc\",\n",
    ")\n",
    "\n",
    "\n",
    "# If you want the ticks of the inset to match the main axis\n",
    "# axins.set_yticks(ax.get_yticks())\n",
    "# Adjust tick parameters to put the labels inside\n",
    "axins1.tick_params(axis=\"x\", direction=\"in\", pad=-15)\n",
    "axins1.set_ylabel(\"Final AMOC Max (Sv)\", size=axlbls)\n",
    "axins1.set_xlabel(\"Member Count\", size=axlbls)\n",
    "axins1.xaxis.set_label_coords(0.5, -0.04)  # Halfway along the x-axis, slightly below\n",
    "axins1.yaxis.set_label_coords(0.73, 0.5)  # Halfway along the x-axis, slightly below\n",
    "# Remove the yticklabels (or move right as below)\n",
    "# axins.set_yticklabels([])\n",
    "# Move the y-tick labels to the right side of the plot\n",
    "axins1.yaxis.set_label_position(\"right\")\n",
    "axins1.yaxis.tick_right()\n",
    "axins1.tick_params(axis=\"y\", direction=\"in\", pad=-20)\n",
    "setaxesfsize(axins1, 12)\n",
    "\n",
    "# Fetch current x-tick labels from inner plot\n",
    "xticks = axins1.get_xticks()\n",
    "labels = [item.get_text() for item in axins1.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins1.set_xticks(xticks)\n",
    "axins1.set_xticklabels(labels)\n",
    "\n",
    "# Fetch current y-tick labels from inner plot\n",
    "yticks = axins1.get_yticks()\n",
    "labels = [item.get_text() for item in axins1.get_yticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins1.set_yticks(yticks)\n",
    "axins1.set_yticklabels(labels)\n",
    "\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "axins1.grid()\n",
    "# lower saddle node line\n",
    "#ax1.axhline(\n",
    "#    15.5,\n",
    "#    0,\n",
    "#    1,\n",
    "#    color=\"k\",\n",
    "#    alpha=0.5,\n",
    "#    linewidth=1.0,\n",
    "#    linestyle=\"dotted\",\n",
    "#    label=\"Lower Saddle Tipping Point\",\n",
    "#)\n",
    "\n",
    "threshold_val = 18\n",
    "def calc_tprob(data, threshold_val):\n",
    "    # Convert data to a Pandas Series if it's not already one\n",
    "    data_series = pd.Series(data)\n",
    "    # Count elements above and below the threshold\n",
    "    count_above = data_series[data_series > threshold_val].count()\n",
    "    count_below = data_series[data_series < threshold_val].count()\n",
    "    # Total number of runs\n",
    "    total_runs = 50\n",
    "    # Calculate the transition probability\n",
    "    # The percentage of runs that are above the threshold\n",
    "    tprob = 100.0 * count_below / total_runs\n",
    "    return tprob\n",
    "\n",
    "# Optionally, add a line at the threshold and saddle nodes\n",
    "#axins1.axhline(y=threshold_val, color='green', linestyle='--')\n",
    "#plt.axhline(y=20.5, color='black', linestyle='--')\n",
    "#plt.axhline(y=15.5, color='black', linestyle='--')\n",
    "\n",
    "# Add transition text\n",
    "x_coord = -1.5\n",
    "y_coord = 0.93\n",
    "tprob_volc = calc_tprob(ens_eq_3701_moc_max_mean_last10yr, threshold_val)\n",
    "axins1.text(x_coord, y_coord, f'Volcanic Transitions: {tprob_volc:.0f}%',\n",
    "            transform=axins1.transAxes,  # Use axes fraction for coordinates\n",
    "            size=vitfont,  # Font size\n",
    "            color='r',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "x_coord = -0.15\n",
    "tprob_novolc = calc_tprob(ens_eq_no3701_moc_max_mean_last10yr, threshold_val)\n",
    "axins1.text(x_coord, y_coord, f'Noise Transitions: {tprob_novolc:.0f}%',\n",
    "            transform=axins1.transAxes,  # Use axes fraction for coordinates\n",
    "            size=nitfont,  # Font size\n",
    "            color='b',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "\n",
    "######################################################\n",
    "################## Zoom Plot2 ########################\n",
    "\n",
    "# original lines\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax2,\n",
    "    color=\"blue\",\n",
    "    linewidth=origlw,\n",
    "    alpha=esalp,\n",
    "    label=\"Control\",\n",
    ")\n",
    "#eq_v1000_moc_max_rm10a.plot(\n",
    "#    ax=ax2,\n",
    "#    color=\"red\",\n",
    "#    linewidth=origlw,\n",
    "#    alpha=esalp,\n",
    "#    label=\"CO$_2$ = 230ppm (Periodic EQ 500Tg SO$_4$)\",\n",
    "#)\n",
    "# esemble member zero , for figure label\n",
    "# only use first 225 years\n",
    "# [0:2700]\n",
    "#endmon = 2400\n",
    "ens_eq_no6701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax2,\n",
    "    color=\"blue\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Noise Ensemble\",\n",
    ")\n",
    "ens_eq_6701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax2,\n",
    "    color=\"red\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Volcanic Ensemble\"\n",
    ")\n",
    "# other esemble members\n",
    "[\n",
    "    ens_eq_no6701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax2, color=\"blue\", linewidth=eslw, alpha=esalp\n",
    "    )\n",
    "    for i in range(1, nruns_6701)\n",
    "]\n",
    "[\n",
    "    ens_eq_6701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax2, color=\"red\", linewidth=eslw, alpha=esalp)\n",
    "    for i in range(1, nruns_6701)\n",
    "]\n",
    "\n",
    "#ensmean_eq_no6701_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax2,\n",
    "#    color=\"cyan\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Noise (Ens Mean)\",\n",
    "#)\n",
    "#ensmean_eq_6701_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax2,\n",
    "#    color=\"yellow\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Volc (Ens Mean)\",\n",
    "#)\n",
    "\n",
    "# Volcano at 6701 line\n",
    "eq_v6701 = cftime.num2date(\n",
    "    6701 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ymin = 8\n",
    "ymax = 26\n",
    "ax2.vlines(eq_v6701, ymin, ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1)\n",
    "setaxesfsize(ax2, axts)\n",
    "ax2.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "ax2.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "ax2.set_ylim(ymin, ymax)\n",
    "\n",
    "xminval = 6600\n",
    "xmaxval = 7000\n",
    "xmin = cftime.num2date(\n",
    "    xminval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "xmax = cftime.num2date(\n",
    "    xmaxval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ax2.set_xlim(xmin, xmax)\n",
    "\n",
    "# Fix the ticks and tick labels to make it look like inset is a separate plot\n",
    "# Fetch current x-tick labels from outer plot\n",
    "xticks = ax2.get_xticks()\n",
    "labels = [item.get_text() for item in ax2.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    # labels[0] = ''\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "ax2.set_xticks(xticks)\n",
    "ax2.set_xticklabels(labels)\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "ax2.grid()\n",
    "\n",
    "vefc='#ffffff'\n",
    "vealp=1.0\n",
    "ax2.legend(loc=3, fontsize=10, facecolor=vefc, framealpha=vealp)\n",
    "ax2.annotate('Volcanic\\nEvent',\n",
    "             xy=vlabxy, xycoords='axes fraction',\n",
    "             xytext=vlabxyt, textcoords='axes fraction',\n",
    "             arrowprops=dict(facecolor='red', shrink=0.05),\n",
    "             horizontalalignment='left', verticalalignment='top',\n",
    "             size=12,\n",
    "             bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))\n",
    "\n",
    "ax2.tick_params(bottom=False, right=False)\n",
    "\n",
    "ax2.set_title(\"\")\n",
    "\n",
    "####### INSET\n",
    "# Create the inset axes, choosing a position and size manually\n",
    "# The parameters are [left, bottom, width, height] in the figure’s coordinate system\n",
    "# inset axes....\n",
    "#axins2 = ax3.inset_axes([0.75, 0.0, 0.25, 1.0])\n",
    "x1, x2, y1, y2 = 0, 20, ymin, ymax\n",
    "axins2.set_xlim(x1, x2)\n",
    "axins2.set_ylim(y1, y2)\n",
    "\n",
    "bwidth = 0.5\n",
    "sns.histplot(\n",
    "    ax=axins2,\n",
    "    y=ens_eq_no6701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"blue\",\n",
    "    alpha=0.3,\n",
    "    label=\"No Volc\",\n",
    ")\n",
    "\n",
    "# sm.ProbPlot(beforedata[0][0]).qqplot(line='s', ax=ax[1])\n",
    "sns.histplot(\n",
    "    ax=axins2,\n",
    "    y=ens_eq_6701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"red\",\n",
    "    alpha=0.3,\n",
    "    label=\"Volc\",\n",
    ")\n",
    "\n",
    "# If you want the ticks of the inset to match the main axis\n",
    "# axins.set_yticks(ax.get_yticks())\n",
    "# Adjust tick parameters to put the labels inside\n",
    "axins2.tick_params(axis=\"x\", direction=\"in\", pad=-15)\n",
    "axins2.set_ylabel(\"Final AMOC Max (Sv)\", size=axlbls)\n",
    "axins2.set_xlabel(\"Member Count\", size=axlbls)\n",
    "axins2.xaxis.set_label_coords(0.5, -0.04)  # Halfway along the x-axis, slightly below\n",
    "axins2.yaxis.set_label_coords(0.73, 0.5)  # Halfway along the x-axis, slightly below\n",
    "# Remove the yticklabels (or move right as below)\n",
    "# axins.set_yticklabels([])\n",
    "# Move the y-tick labels to the right side of the plot\n",
    "axins2.yaxis.set_label_position(\"right\")\n",
    "axins2.yaxis.tick_right()\n",
    "axins2.tick_params(axis=\"y\", direction=\"in\", pad=-20)\n",
    "setaxesfsize(axins2, 12)\n",
    "\n",
    "# Fetch current x-tick labels from inner plot\n",
    "xticks = axins2.get_xticks()\n",
    "labels = [item.get_text() for item in axins2.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins2.set_xticks(xticks)\n",
    "axins2.set_xticklabels(labels)\n",
    "\n",
    "# Fetch current y-tick labels from inner plot\n",
    "yticks = axins2.get_yticks()\n",
    "labels = [item.get_text() for item in axins2.get_yticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins2.set_yticks(yticks)\n",
    "axins2.set_yticklabels(labels)\n",
    "\n",
    "threshold_val = 18\n",
    "def calc_tprob(data, threshold_val):\n",
    "    # Convert data to a Pandas Series if it's not already one\n",
    "    data_series = pd.Series(data)\n",
    "    # Count elements above and below the threshold\n",
    "    count_above = data_series[data_series > threshold_val].count()\n",
    "    count_below = data_series[data_series < threshold_val].count()\n",
    "    # Total number of runs\n",
    "    total_runs = 50\n",
    "    # Calculate the transition probability\n",
    "    # The percentage of runs that are above the threshold\n",
    "    tprob = 100.0 * count_below / total_runs\n",
    "    return tprob\n",
    "\n",
    "# Optionally, add a line at the threshold and saddle nodes\n",
    "#axins2.axhline(y=threshold_val, color='green', linestyle='--')\n",
    "#plt.axhline(y=20.5, color='black', linestyle='--')\n",
    "#plt.axhline(y=15.5, color='black', linestyle='--')\n",
    "\n",
    "# Add transition text\n",
    "x_coord = -1.5\n",
    "y_coord = 0.93\n",
    "tprob_volc = calc_tprob(ens_eq_6701_moc_max_mean_last10yr, threshold_val)\n",
    "axins2.text(x_coord, y_coord, f'Volcanic Transitions: {tprob_volc:.0f}%',\n",
    "            transform=axins2.transAxes,  # Use axes fraction for coordinates\n",
    "            size=vitfont,  # Font size\n",
    "            color='r',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "x_coord = -0.15\n",
    "tprob_novolc = calc_tprob(ens_eq_no6701_moc_max_mean_last10yr, threshold_val)\n",
    "axins2.text(x_coord, y_coord, f'Noise Transitions: {tprob_novolc:.0f}%',\n",
    "            transform=axins2.transAxes,  # Use axes fraction for coordinates\n",
    "            size=nitfont,  # Font size\n",
    "            color='b',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "axins2.grid()\n",
    "\n",
    "######################################################\n",
    "################## Zoom Plot3 ########################\n",
    "\n",
    "# original lines\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax3,\n",
    "    color=\"blue\",\n",
    "    linewidth=origlw,\n",
    "    alpha=esalp,\n",
    "    label=\"Control\",\n",
    ")\n",
    "#eq_v1000_moc_max_rm10a.plot(\n",
    "#    ax=ax3,\n",
    "#    color=\"red\",\n",
    "#    linewidth=origlw,\n",
    "#    alpha=esalp,\n",
    "#    label=\"CO$_2$ = 230ppm (Periodic EQ 500Tg SO$_4$)\",\n",
    "#)\n",
    "# esemble member zero , for figure label\n",
    "# only use first 225 years\n",
    "# [0:2700]\n",
    "#endmon = 2400\n",
    "ens_eq_no7701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax3,\n",
    "    color=\"blue\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Noise Ensemble\",\n",
    ")\n",
    "ens_eq_7701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax3,\n",
    "    color=\"red\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Volcanic Ensemble\"\n",
    ")\n",
    "# other esemble members\n",
    "[\n",
    "    ens_eq_no7701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax3, color=\"blue\", linewidth=eslw, alpha=esalp\n",
    "    )\n",
    "    for i in range(1, nruns_6701)\n",
    "]\n",
    "[\n",
    "    ens_eq_7701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax3, color=\"red\", linewidth=eslw, alpha=esalp)\n",
    "    for i in range(1, nruns_6701)\n",
    "]\n",
    "\n",
    "#ensmean_eq_no7701_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax3,\n",
    "#    color=\"cyan\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Noise (Ens Mean)\",\n",
    "#)\n",
    "#ensmean_eq_7701_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax3,\n",
    "#    color=\"yellow\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Volc (Ens Mean)\",\n",
    "#)\n",
    "\n",
    "# Volcano at 6701 line\n",
    "eq_v7701 = cftime.num2date(\n",
    "    7701 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ymin = 8\n",
    "ymax = 26\n",
    "ax3.vlines(eq_v7701, ymin, ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1)\n",
    "setaxesfsize(ax3, axts)\n",
    "ax3.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "ax3.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "ax3.set_ylim(ymin, ymax)\n",
    "\n",
    "xminval = 7600\n",
    "xmaxval = 8000\n",
    "xmin = cftime.num2date(\n",
    "    xminval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "xmax = cftime.num2date(\n",
    "    xmaxval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ax3.set_xlim(xmin, xmax)\n",
    "\n",
    "# Fix the ticks and tick labels to make it look like inset is a separate plot\n",
    "# Fetch current x-tick labels from outer plot\n",
    "xticks = ax3.get_xticks()\n",
    "labels = [item.get_text() for item in ax3.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    # labels[0] = ''\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "ax3.set_xticks(xticks)\n",
    "ax3.set_xticklabels(labels)\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "ax3.grid()\n",
    "\n",
    "#vefc='white'\n",
    "#vealp=0.8\n",
    "vefc='#ffffff'\n",
    "vealp=1.0\n",
    "ax3.legend(loc=3, fontsize=10, facecolor=vefc, framealpha=vealp)\n",
    "ax3.annotate('Volcanic\\nEvent',\n",
    "             xy=vlabxy, xycoords='axes fraction',\n",
    "             xytext=vlabxyt, textcoords='axes fraction',\n",
    "             arrowprops=dict(facecolor='red', shrink=0.05),\n",
    "             horizontalalignment='left', verticalalignment='top',\n",
    "             size=12,\n",
    "             bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))\n",
    "\n",
    "ax3.tick_params(bottom=False, right=False)\n",
    "\n",
    "ax3.set_title(\"\")\n",
    "\n",
    "####### INSET\n",
    "# Create the inset axes, choosing a position and size manually\n",
    "# The parameters are [left, bottom, width, height] in the figure’s coordinate system\n",
    "# inset axes....\n",
    "#axins3 = ax3.inset_axes([0.75, 0.0, 0.25, 1.0])\n",
    "x1, x2, y1, y2 = 0, 20, ymin, ymax\n",
    "axins3.set_xlim(x1, x2)\n",
    "axins3.set_ylim(y1, y2)\n",
    "\n",
    "bwidth = 0.5\n",
    "sns.histplot(\n",
    "    ax=axins3,\n",
    "    y=ens_eq_no7701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"blue\",\n",
    "    alpha=0.3,\n",
    "    label=\"No Volc\",\n",
    ")\n",
    "\n",
    "# sm.ProbPlot(beforedata[0][0]).qqplot(line='s', ax=ax[1])\n",
    "sns.histplot(\n",
    "    ax=axins3,\n",
    "    y=ens_eq_7701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"red\",\n",
    "    alpha=0.3,\n",
    "    label=\"Volc\",\n",
    ")\n",
    "\n",
    "# If you want the ticks of the inset to match the main axis\n",
    "# axins.set_yticks(ax.get_yticks())\n",
    "# Adjust tick parameters to put the labels inside\n",
    "axins3.tick_params(axis=\"x\", direction=\"in\", pad=-15)\n",
    "axins3.set_ylabel(\"Final AMOC Max (Sv)\", size=axlbls)\n",
    "axins3.set_xlabel(\"Member Count\", size=axlbls)\n",
    "axins3.xaxis.set_label_coords(0.5, -0.04)  # Halfway along the x-axis, slightly below\n",
    "axins3.yaxis.set_label_coords(0.73, 0.5)  # Halfway along the x-axis, slightly below\n",
    "# Remove the yticklabels (or move right as below)\n",
    "# axins.set_yticklabels([])\n",
    "# Move the y-tick labels to the right side of the plot\n",
    "axins3.yaxis.set_label_position(\"right\")\n",
    "axins3.yaxis.tick_right()\n",
    "axins3.tick_params(axis=\"y\", direction=\"in\", pad=-20)\n",
    "setaxesfsize(axins3, 12)\n",
    "\n",
    "# Fetch current x-tick labels from inner plot\n",
    "xticks = axins3.get_xticks()\n",
    "labels = [item.get_text() for item in axins3.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins3.set_xticks(xticks)\n",
    "axins3.set_xticklabels(labels)\n",
    "\n",
    "# Fetch current y-tick labels from inner plot\n",
    "yticks = axins3.get_yticks()\n",
    "labels = [item.get_text() for item in axins3.get_yticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins3.set_yticks(yticks)\n",
    "axins3.set_yticklabels(labels)\n",
    "\n",
    "threshold_val = 18\n",
    "def calc_tprob(data, threshold_val):\n",
    "    # Convert data to a Pandas Series if it's not already one\n",
    "    data_series = pd.Series(data)\n",
    "    # Count elements above and below the threshold\n",
    "    count_above = data_series[data_series > threshold_val].count()\n",
    "    count_below = data_series[data_series < threshold_val].count()\n",
    "    # Total number of runs\n",
    "    total_runs = 50\n",
    "    # Calculate the transition probability\n",
    "    # The percentage of runs that are above the threshold\n",
    "    tprob = 100.0 * count_below / total_runs\n",
    "    return tprob\n",
    "\n",
    "# Optionally, add a line at the threshold and saddle nodes\n",
    "#axins3.axhline(y=threshold_val, color='green', linestyle='--')\n",
    "#plt.axhline(y=20.5, color='black', linestyle='--')\n",
    "#plt.axhline(y=15.5, color='black', linestyle='--')\n",
    "\n",
    "# Add transition text\n",
    "x_coord = -1.5\n",
    "y_coord = 0.93\n",
    "tprob_volc = calc_tprob(ens_eq_7701_moc_max_mean_last10yr, threshold_val)\n",
    "axins3.text(x_coord, y_coord, f'Volcanic Transitions: {tprob_volc:.0f}%',\n",
    "            transform=axins3.transAxes,  # Use axes fraction for coordinates\n",
    "            size=vitfont,  # Font size\n",
    "            color='r',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "x_coord = -0.15\n",
    "tprob_novolc = calc_tprob(ens_eq_no7701_moc_max_mean_last10yr, threshold_val)\n",
    "axins3.text(x_coord, y_coord, f'Noise Transitions: {tprob_novolc:.0f}%',\n",
    "            transform=axins3.transAxes,  # Use axes fraction for coordinates\n",
    "            size=nitfont,  # Font size\n",
    "            color='b',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "axins3.grid()\n",
    "\n",
    "######################################################\n",
    "################## Zoom Plot4 ########################\n",
    "sns.set_style(\"darkgrid\")\n",
    "\n",
    "esalp = 0.6  # ensemble line transparency\n",
    "esmalp = 0.9  # ensemble mean line tranparency\n",
    "eslw = 0.5  # ensemble line width\n",
    "esmlw = 2.0  # ensemble mean line width\n",
    "origlw = 2.0  # original runs line width\n",
    "\n",
    "# original lines\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax1_dup,\n",
    "    color=\"blue\",\n",
    "    linewidth=origlw,\n",
    "    alpha=esalp,\n",
    "    label=\"Control\",\n",
    ")\n",
    "# esemble member zero , for figure label\n",
    "ens_eq_no8501_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax1_dup,\n",
    "    color=\"blue\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Noise Ensemble\",\n",
    ")\n",
    "ens_eq_8501_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax1_dup,\n",
    "    color=\"red\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Volcanic Ensemble\"\n",
    ")\n",
    "# other esemble members\n",
    "[\n",
    "    ens_eq_no8501_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax1_dup, color=\"blue\", linewidth=eslw, alpha=esalp\n",
    "    )\n",
    "    for i in range(1, nruns_8501)\n",
    "]\n",
    "[\n",
    "    ens_eq_8501_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax1_dup, color=\"red\", linewidth=eslw, alpha=esalp)\n",
    "    for i in range(1, nruns_8501)\n",
    "]\n",
    "\n",
    "#ensmean_eq_no8501_moc_max_rm10a.plot(\n",
    "#    ax=ax1_dup,\n",
    "#    color=\"cyan\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Noise (Ens Mean)\",\n",
    "#)\n",
    "#ensmean_eq_8501_moc_max_rm10a.plot(\n",
    "#    ax=ax1_dup,\n",
    "#    color=\"yellow\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Volc (Ens Mean)\",\n",
    "#)\n",
    "\n",
    "# Volcano at 5701 line\n",
    "eq_v8501 = cftime.num2date(\n",
    "    8501 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ymin = 8\n",
    "ymax = 26\n",
    "\n",
    "ax1_dup.vlines(eq_v8501, ymin, ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1)\n",
    "setaxesfsize(ax1_dup, axts)\n",
    "ax1_dup.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "#ax1_dup.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "ax1_dup.set_ylabel(\"\", size=axlbls)\n",
    "ax1_dup.set_ylim(ymin, ymax)\n",
    "\n",
    "xminval = 8400\n",
    "xmaxval = 8800\n",
    "xmin = cftime.num2date(\n",
    "    xminval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "xmax = cftime.num2date(\n",
    "    xmaxval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ax1_dup.set_xlim(xmin, xmax)\n",
    "\n",
    "# Fix the ticks and tick labels to make it look like inset is a separate plot\n",
    "# Fetch current x-tick labels from outer plot\n",
    "xticks = ax1_dup.get_xticks()\n",
    "labels = [item.get_text() for item in ax1_dup.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    # labels[0] = ''\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "ax1_dup.set_xticks(xticks)\n",
    "ax1_dup.set_xticklabels(labels)\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "ax1_dup.grid()\n",
    "\n",
    "vefc='#ffffff'\n",
    "vealp=1.0\n",
    "#ax1.legend(loc=3, facecolor=vefc, framealpha=vealp)\n",
    "ax1_dup.annotate('Volcanic\\nEvent',\n",
    "             xy=vlabxy, xycoords='axes fraction',\n",
    "             xytext=vlabxyt, textcoords='axes fraction',\n",
    "             arrowprops=dict(facecolor='red', shrink=0.05),\n",
    "             horizontalalignment='left', verticalalignment='top',\n",
    "             size=12,\n",
    "             bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))\n",
    "\n",
    "ax1_dup.tick_params(bottom=False, right=False)\n",
    "\n",
    "ax1_dup.set_title(\"\")\n",
    "\n",
    "\n",
    "####### INSET\n",
    "# Create the inset axes, choosing a position and size manually\n",
    "# The parameters are [left, bottom, width, height] in the figure’s coordinate system\n",
    "# inset axes....\n",
    "# change this to its own gridpsect now\n",
    "#axins1 = ax1.inset_axes([0.75, 0.0, 0.25, 1.0])\n",
    "x1, x2, y1, y2 = xinmin, xinmax, ymin, ymax\n",
    "axins1_dup.set_xlim(x1, x2)\n",
    "axins1_dup.set_ylim(y1, y2)\n",
    "\n",
    "bwidth = 0.5\n",
    "sns.histplot(\n",
    "    ax=axins1_dup,\n",
    "    y=ens_eq_no8501_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"blue\",\n",
    "    alpha=0.3,\n",
    "    label=\"No Volc\",\n",
    ")\n",
    "# sm.ProbPlot(beforedata[0][0]).qqplot(line='s', ax=ax[1])\n",
    "sns.histplot(\n",
    "    ax=axins1_dup,\n",
    "    y=ens_eq_8501_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"red\",\n",
    "    alpha=0.3,\n",
    "    label=\"Volc\",\n",
    ")\n",
    "\n",
    "# If you want the ticks of the inset to match the main axis\n",
    "# axins.set_yticks(ax.get_yticks())\n",
    "# Adjust tick parameters to put the labels inside\n",
    "axins1_dup.tick_params(axis=\"x\", direction=\"in\", pad=-15)\n",
    "axins1_dup.set_ylabel(\"Final AMOC Max (Sv)\", size=axlbls)\n",
    "axins1_dup.set_xlabel(\"Member Count\", size=axlbls)\n",
    "axins1_dup.xaxis.set_label_coords(0.5, -0.04)  # Halfway along the x-axis, slightly below\n",
    "axins1_dup.yaxis.set_label_coords(0.73, 0.5)  # Halfway along the x-axis, slightly below\n",
    "# Remove the yticklabels (or move right as below)\n",
    "# axins_dup.set_yticklabels([])\n",
    "# Move the y-tick labels to the right side of the plot\n",
    "axins1_dup.yaxis.set_label_position(\"right\")\n",
    "axins1_dup.yaxis.tick_right()\n",
    "axins1_dup.tick_params(axis=\"y\", direction=\"in\", pad=-20)\n",
    "setaxesfsize(axins1_dup, 12)\n",
    "\n",
    "# Fetch current x-tick labels from inner plot\n",
    "xticks = axins1_dup.get_xticks()\n",
    "labels = [item.get_text() for item in axins1_dup.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins1_dup.set_xticks(xticks)\n",
    "axins1_dup.set_xticklabels(labels)\n",
    "\n",
    "# Fetch current y-tick labels from inner plot\n",
    "yticks = axins1_dup.get_yticks()\n",
    "labels = [item.get_text() for item in axins1_dup.get_yticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins1_dup.set_yticks(yticks)\n",
    "axins1_dup.set_yticklabels(labels)\n",
    "\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "axins1_dup.grid()\n",
    "# lower saddle node line\n",
    "#ax1_dup.axhline(\n",
    "#    15.5,\n",
    "#    0,\n",
    "#    1,\n",
    "#    color=\"k\",\n",
    "#    alpha=0.5,\n",
    "#    linewidth=1.0,\n",
    "#    linestyle=\"dotted\",\n",
    "#    label=\"Lower Saddle Tipping Point\",\n",
    "#)\n",
    "\n",
    "threshold_val = 18\n",
    "def calc_tprob(data, threshold_val):\n",
    "    # Convert data to a Pandas Series if it's not already one\n",
    "    data_series = pd.Series(data)\n",
    "    # Count elements above and below the threshold\n",
    "    count_above = data_series[data_series < threshold_val].count()\n",
    "    count_below = data_series[data_series > threshold_val].count()\n",
    "    # Total number of runs\n",
    "    total_runs = 50\n",
    "    # Calculate the transition probability\n",
    "    # The percentage of runs that are above the threshold\n",
    "    tprob = 100.0 * count_below / total_runs\n",
    "    return tprob\n",
    "\n",
    "# Optionally, add a line at the threshold and saddle nodes\n",
    "#axins1_dup.axhline(y=threshold_val, color='green', linestyle='--')\n",
    "#plt.axhline(y=20.5, color='black', linestyle='--')\n",
    "#plt.axhline(y=15.5, color='black', linestyle='--')\n",
    "\n",
    "# Add transition text\n",
    "x_coord = -1.5\n",
    "y_coord = 0.93\n",
    "tprob_volc = calc_tprob(ens_eq_8501_moc_max_mean_last10yr, threshold_val)\n",
    "axins1_dup.text(x_coord, y_coord, f'Volcanic Transitions: {tprob_volc:.0f}%',\n",
    "            transform=axins1_dup.transAxes,  # Use axes fraction for coordinates\n",
    "            size=vitfont,  # Font size\n",
    "            color='r',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "x_coord = -0.15\n",
    "tprob_novolc = calc_tprob(ens_eq_no8501_moc_max_mean_last10yr, threshold_val)\n",
    "axins1_dup.text(x_coord, y_coord, f'Noise Transitions: {tprob_novolc:.0f}%',\n",
    "            transform=axins1_dup.transAxes,  # Use axes fraction for coordinates\n",
    "            size=nitfont,  # Font size\n",
    "            color='b',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "\n",
    "######################################################\n",
    "################## Zoom Plot5 ########################\n",
    "\n",
    "# original lines\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax2_dup,\n",
    "    color=\"blue\",\n",
    "    linewidth=origlw,\n",
    "    alpha=esalp,\n",
    "    label=\"Control\",\n",
    ")\n",
    "ens_eq_no8701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax2_dup,\n",
    "    color=\"blue\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"NoVol Ens\",\n",
    ")\n",
    "ens_eq_8701_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax2_dup,\n",
    "    color=\"red\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Volcanic Ensemble\"\n",
    ")\n",
    "# other esemble members\n",
    "[\n",
    "    ens_eq_no8701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax2_dup, color=\"blue\", linewidth=eslw, alpha=esalp\n",
    "    )\n",
    "    for i in range(1, nruns_8701)\n",
    "]\n",
    "[\n",
    "    ens_eq_8701_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax2_dup, color=\"red\", linewidth=eslw, alpha=esalp)\n",
    "    for i in range(1, nruns_8701)\n",
    "]\n",
    "\n",
    "#ensmean_eq_no8701_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax2_dup,\n",
    "#    color=\"cyan\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Noise (Ens Mean)\",\n",
    "#)\n",
    "#ensmean_eq_8701_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax2_dup,\n",
    "#    color=\"yellow\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Volc (Ens Mean)\",\n",
    "#)\n",
    "\n",
    "# Volcano at 8701 line\n",
    "eq_v8701 = cftime.num2date(\n",
    "    8701 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ymin = 10\n",
    "ymax = 28\n",
    "ax2_dup.vlines(eq_v8701, ymin, ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1)\n",
    "setaxesfsize(ax2_dup, axts)\n",
    "ax2_dup.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "#ax2_dup.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "ax2_dup.set_ylabel(\"\", size=axlbls)\n",
    "ax2_dup.set_ylim(ymin, ymax)\n",
    "\n",
    "xminval = 8600\n",
    "xmaxval = 9000\n",
    "xmin = cftime.num2date(\n",
    "    xminval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "xmax = cftime.num2date(\n",
    "    xmaxval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ax2_dup.set_xlim(xmin, xmax)\n",
    "\n",
    "# Fix the ticks and tick labels to make it look like inset is a separate plot\n",
    "# Fetch current x-tick labels from outer plot\n",
    "xticks = ax2_dup.get_xticks()\n",
    "labels = [item.get_text() for item in ax2_dup.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    # labels[0] = ''\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "ax2_dup.set_xticks(xticks)\n",
    "ax2_dup.set_xticklabels(labels)\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "ax2_dup.grid()\n",
    "\n",
    "vefc='#ffffff'\n",
    "vealp=1.0\n",
    "#ax2_dup.legend(loc=3, facecolor=vefc, framealpha=vealp)\n",
    "ax2_dup.annotate('Volcanic\\nEvent',\n",
    "             xy=vlabxy, xycoords='axes fraction',\n",
    "             xytext=vlabxyt, textcoords='axes fraction',\n",
    "             arrowprops=dict(facecolor='red', shrink=0.05),\n",
    "             horizontalalignment='left', verticalalignment='top',\n",
    "             size=12,\n",
    "             bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))\n",
    "\n",
    "ax2_dup.tick_params(bottom=False, right=False)\n",
    "\n",
    "ax2_dup.set_title(\"\")\n",
    "\n",
    "####### INSET\n",
    "# Create the inset axes, choosing a position and size manually\n",
    "# The parameters are [left, bottom, width, height] in the figure’s coordinate system\n",
    "# inset axes....\n",
    "#axins2 = ax3.inset_axes([0.75, 0.0, 0.25, 1.0])\n",
    "x1, x2, y1, y2 = xinmin, xinmax, ymin, ymax\n",
    "axins2_dup.set_xlim(x1, x2)\n",
    "axins2_dup.set_ylim(y1, y2)\n",
    "\n",
    "bwidth = 0.5\n",
    "sns.histplot(\n",
    "    ax=axins2_dup,\n",
    "    y=ens_eq_no8701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"blue\",\n",
    "    alpha=0.3,\n",
    "    label=\"No Volc\",\n",
    ")\n",
    "\n",
    "# sm.ProbPlot(beforedata[0][0]).qqplot(line='s', ax=ax[1])\n",
    "sns.histplot(\n",
    "    ax=axins2_dup,\n",
    "    y=ens_eq_8701_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"red\",\n",
    "    alpha=0.3,\n",
    "    label=\"Volc\",\n",
    ")\n",
    "\n",
    "# If you want the ticks of the inset to match the main axis\n",
    "# axins.set_yticks(ax.get_yticks())\n",
    "# Adjust tick parameters to put the labels inside\n",
    "axins2_dup.tick_params(axis=\"x\", direction=\"in\", pad=-15)\n",
    "axins2_dup.set_ylabel(\"Final AMOC Max (Sv)\", size=axlbls)\n",
    "axins2_dup.set_xlabel(\"Member Count\", size=axlbls)\n",
    "axins2_dup.xaxis.set_label_coords(0.5, -0.04)  # Halfway along the x-axis, slightly below\n",
    "axins2_dup.yaxis.set_label_coords(0.73, 0.5)  # Halfway along the x-axis, slightly below\n",
    "# Remove the yticklabels (or move right as below)\n",
    "# axins.set_yticklabels([])\n",
    "# Move the y-tick labels to the right side of the plot\n",
    "axins2_dup.yaxis.set_label_position(\"right\")\n",
    "axins2_dup.yaxis.tick_right()\n",
    "axins2_dup.tick_params(axis=\"y\", direction=\"in\", pad=-20)\n",
    "setaxesfsize(axins2_dup, 12)\n",
    "\n",
    "# Fetch current x-tick labels from inner plot\n",
    "xticks = axins2_dup.get_xticks()\n",
    "labels = [item.get_text() for item in axins2_dup.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins2_dup.set_xticks(xticks)\n",
    "axins2_dup.set_xticklabels(labels)\n",
    "\n",
    "# Fetch current y-tick labels from inner plot\n",
    "yticks = axins2_dup.get_yticks()\n",
    "labels = [item.get_text() for item in axins2_dup.get_yticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins2_dup.set_yticks(yticks)\n",
    "axins2_dup.set_yticklabels(labels)\n",
    "\n",
    "threshold_val = 18\n",
    "def calc_tprob(data, threshold_val):\n",
    "    # Convert data to a Pandas Series if it's not already one\n",
    "    data_series = pd.Series(data)\n",
    "    # Count elements above and below the threshold\n",
    "    count_above = data_series[data_series < threshold_val].count()\n",
    "    count_below = data_series[data_series > threshold_val].count()\n",
    "    # Total number of runs\n",
    "    total_runs = 50\n",
    "    # Calculate the transition probability\n",
    "    # The percentage of runs that are above the threshold\n",
    "    tprob = 100.0 * count_below / total_runs\n",
    "    return tprob\n",
    "\n",
    "# Optionally, add a line at the threshold and saddle nodes\n",
    "#axins2_dup.axhline(y=threshold_val, color='green', linestyle='--')\n",
    "#plt.axhline(y=20.5, color='black', linestyle='--')\n",
    "#plt.axhline(y=15.5, color='black', linestyle='--')\n",
    "\n",
    "# Add transition text\n",
    "x_coord = -1.5\n",
    "y_coord = 0.93\n",
    "tprob_volc = calc_tprob(ens_eq_8701_moc_max_mean_last10yr, threshold_val)\n",
    "axins2_dup.text(x_coord, y_coord, f'Volcanic Transitions: {tprob_volc:.0f}%',\n",
    "            transform=axins2_dup.transAxes,  # Use axes fraction for coordinates\n",
    "            color='r',\n",
    "            size=vitfont,  # Font size\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "x_coord = -0.15\n",
    "tprob_novolc = calc_tprob(ens_eq_no8701_moc_max_mean_last10yr, threshold_val)\n",
    "axins2_dup.text(x_coord, y_coord, f'Noise Transitions: {tprob_novolc:.0f}%',\n",
    "            transform=axins2_dup.transAxes,  # Use axes fraction for coordinates\n",
    "            size=nitfont,  # Font size\n",
    "            color='b',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "axins2_dup.grid()\n",
    "\n",
    "######################################################\n",
    "################## Zoom Plot6 ########################\n",
    "\n",
    "# original lines\n",
    "ctrl230i_moc_max_rm10a.plot(\n",
    "    ax=ax3_dup,\n",
    "    color=\"blue\",\n",
    "    linewidth=origlw,\n",
    "    alpha=esalp,\n",
    "    label=\"Control\",\n",
    ")\n",
    "ens_eq_no8901_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax3_dup,\n",
    "    color=\"blue\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Noise Ensemble\",\n",
    ")\n",
    "ens_eq_8901_moc_max_rm10a[0][0:endmon].plot(\n",
    "    ax=ax3_dup,\n",
    "    color=\"red\",\n",
    "    linewidth=eslw,\n",
    "    alpha=esalp,\n",
    "    label=\"Volcanic Ensemble\"\n",
    ")\n",
    "# other esemble members\n",
    "[\n",
    "    ens_eq_no8901_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax3_dup, color=\"blue\", linewidth=eslw, alpha=esalp\n",
    "    )\n",
    "    for i in range(1, nruns_8901)\n",
    "]\n",
    "[\n",
    "    ens_eq_8901_moc_max_rm10a[i][0:endmon].plot(\n",
    "        ax=ax3_dup, color=\"red\", linewidth=eslw, alpha=esalp)\n",
    "    for i in range(1, nruns_8901)\n",
    "]\n",
    "\n",
    "#ensmean_eq_no8901_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax3_dup,\n",
    "#    color=\"cyan\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Noise (Ens Mean)\",\n",
    "#)\n",
    "#ensmean_eq_8901_moc_max_rm10a[0:endmon].plot(\n",
    "#    ax=ax3_dup,\n",
    "#    color=\"yellow\",\n",
    "#    linewidth=esmlw,\n",
    "#    alpha=esmalp,\n",
    "#    label=\"Volc (Ens Mean)\",\n",
    "#)\n",
    "\n",
    "# Volcano at 8901 line\n",
    "eq_v8901 = cftime.num2date(\n",
    "    8901 * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ymin = 12\n",
    "ymax = 30\n",
    "ax3_dup.vlines(eq_v8901, ymin, ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1)\n",
    "setaxesfsize(ax3_dup, axts)\n",
    "ax3_dup.set_xlabel(\"Model Time (Years)\", size=axlbls)\n",
    "#ax3_dup.set_ylabel(\"AMOC Max (Sv)\", size=axlbls)\n",
    "ax3_dup.set_ylabel(\"\", size=axlbls)\n",
    "ax3_dup.set_ylim(ymin, ymax)\n",
    "\n",
    "xminval = 8800\n",
    "xmaxval = 9200\n",
    "xmin = cftime.num2date(\n",
    "    xminval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "xmax = cftime.num2date(\n",
    "    xmaxval * 365, \"days since 0000-01-01 00:00:00\", calendar=\"noleap\"\n",
    ")\n",
    "ax3_dup.set_xlim(xmin, xmax)\n",
    "\n",
    "# Fix the ticks and tick labels to make it look like inset is a separate plot\n",
    "# Fetch current x-tick labels from outer plot\n",
    "xticks = ax3_dup.get_xticks()\n",
    "labels = [item.get_text() for item in ax3_dup.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    # labels[0] = ''\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "ax3_dup.set_xticks(xticks)\n",
    "ax3_dup.set_xticklabels(labels)\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "ax3_dup.grid()\n",
    "\n",
    "#vefc='white'\n",
    "#vealp=0.8\n",
    "vefc='#ffffff'\n",
    "vealp=1.0\n",
    "#ax3.legend(loc=3, facecolor=vefc, framealpha=vealp)\n",
    "ax3_dup.annotate('Volcanic\\nEvent',\n",
    "             xy=vlabxy, xycoords='axes fraction',\n",
    "             xytext=vlabxyt, textcoords='axes fraction',\n",
    "             arrowprops=dict(facecolor='red', shrink=0.05),\n",
    "             horizontalalignment='left', verticalalignment='top',\n",
    "             size=12,\n",
    "             bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))\n",
    "\n",
    "ax3_dup.tick_params(bottom=False, right=False)\n",
    "\n",
    "ax3_dup.set_title(\"\")\n",
    "\n",
    "####### INSET\n",
    "# Create the inset axes, choosing a position and size manually\n",
    "# The parameters are [left, bottom, width, height] in the figure’s coordinate system\n",
    "# inset axes....\n",
    "#axins3 = ax3.inset_axes([0.75, 0.0, 0.25, 1.0])\n",
    "x1, x2, y1, y2 = xinmin, xinmax, ymin, ymax\n",
    "axins3_dup.set_xlim(x1, x2)\n",
    "axins3_dup.set_ylim(y1, y2)\n",
    "\n",
    "bwidth = 0.5\n",
    "sns.histplot(\n",
    "    ax=axins3_dup,\n",
    "    y=ens_eq_no8901_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"blue\",\n",
    "    alpha=0.3,\n",
    "    label=\"No Volc\",\n",
    ")\n",
    "\n",
    "# sm.ProbPlot(beforedata[0][0]).qqplot(line='s', ax=ax[1])\n",
    "sns.histplot(\n",
    "    ax=axins3_dup,\n",
    "    y=ens_eq_8901_moc_max_mean_last10yr,\n",
    "    kde=True,\n",
    "    orientation=\"vertical\",\n",
    "    binwidth=bwidth,\n",
    "    binrange=[ymin, ymax],\n",
    "    color=\"red\",\n",
    "    alpha=0.3,\n",
    "    label=\"Volc\",\n",
    ")\n",
    "\n",
    "# If you want the ticks of the inset to match the main axis\n",
    "# axins.set_yticks(ax.get_yticks())\n",
    "# Adjust tick parameters to put the labels inside\n",
    "axins3_dup.tick_params(axis=\"x\", direction=\"in\", pad=-15)\n",
    "axins3_dup.set_ylabel(\"Final AMOC Max (Sv)\", size=axlbls)\n",
    "axins3_dup.set_xlabel(\"Member Count\", size=axlbls)\n",
    "axins3_dup.xaxis.set_label_coords(0.5, -0.04)  # Halfway along the x-axis, slightly below\n",
    "axins3_dup.yaxis.set_label_coords(0.73, 0.5)  # Halfway along the x-axis, slightly below\n",
    "# Remove the yticklabels (or move right as below)\n",
    "# axins3_dup.set_yticklabels([])\n",
    "# Move the y-tick labels to the right side of the plot\n",
    "axins3_dup.yaxis.set_label_position(\"right\")\n",
    "axins3_dup.yaxis.tick_right()\n",
    "axins3_dup.tick_params(axis=\"y\", direction=\"in\", pad=-20)\n",
    "setaxesfsize(axins3_dup, 12)\n",
    "\n",
    "# Fetch current x-tick labels from inner plot\n",
    "xticks = axins3_dup.get_xticks()\n",
    "labels = [item.get_text() for item in axins3_dup.get_xticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins3_dup.set_xticks(xticks)\n",
    "axins3_dup.set_xticklabels(labels)\n",
    "\n",
    "# Fetch current y-tick labels from inner plot\n",
    "yticks = axins3_dup.get_yticks()\n",
    "labels = [item.get_text() for item in axins3_dup.get_yticklabels()]\n",
    "# Check if there are enough tick labels to modify\n",
    "if len(labels) > 1:\n",
    "    labels[0] = \"\"\n",
    "    # Remove first label\n",
    "    labels[-1] = \"\"\n",
    "    # Remove last label\n",
    "# Set new labels\n",
    "axins3_dup.set_yticks(yticks)\n",
    "axins3_dup.set_yticklabels(labels)\n",
    "\n",
    "threshold_val = 18\n",
    "def calc_tprob(data, threshold_val):\n",
    "    # Convert data to a Pandas Series if it's not already one\n",
    "    data_series = pd.Series(data)\n",
    "    # Count elements above and below the threshold\n",
    "    count_above = data_series[data_series < threshold_val].count()\n",
    "    count_below = data_series[data_series > threshold_val].count()\n",
    "    # Total number of runs\n",
    "    total_runs = 50\n",
    "    # Calculate the transition probability\n",
    "    # The percentage of runs that are above the threshold\n",
    "    tprob = 100.0 * count_below / total_runs\n",
    "    return tprob\n",
    "\n",
    "# Optionally, add a line at the threshold and saddle nodes\n",
    "#axins3_dup.axhline(y=threshold_val, color='green', linestyle='--')\n",
    "#plt.axhline(y=20.5, color='black', linestyle='--')\n",
    "#plt.axhline(y=15.5, color='black', linestyle='--')\n",
    "\n",
    "# Add transition text\n",
    "x_coord = -1.5\n",
    "y_coord = 0.93\n",
    "tprob_volc = calc_tprob(ens_eq_8901_moc_max_mean_last10yr, threshold_val)\n",
    "axins3_dup.text(x_coord, y_coord, f'Volcanic Transitions: {tprob_volc:.0f}%',\n",
    "            transform=axins3_dup.transAxes,  # Use axes fraction for coordinates\n",
    "            color='r',\n",
    "            size=vitfont,  # Font size\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "x_coord = -0.15\n",
    "tprob_novolc = calc_tprob(ens_eq_no8901_moc_max_mean_last10yr, threshold_val)\n",
    "axins3_dup.text(x_coord, y_coord, f'Noise Transitions: {tprob_novolc:.0f}%',\n",
    "            transform=axins3_dup.transAxes,  # Use axes fraction for coordinates\n",
    "            size=nitfont,  # Font size\n",
    "            color='b',\n",
    "            bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor=vefc, alpha=vealp))  # Optional: bounding box\n",
    "# just to compare the same yaxis and make sure the histogram is right\n",
    "axins3_dup.grid()\n",
    "\n",
    "###########################END ZOOM PLOTS ###############################################\n",
    "\n",
    "# Fix up the labels\n",
    "#ax0.set_ylabel(\"AMOC Max (non-dim)\", size=axlbls)\n",
    "#ax0.set_xlabel(\"Model Time (non-dim)\", size=axlbls)\n",
    "ax0.set_title(\"Complex Model Interstadial Age Test\", size=14, loc=\"center\")\n",
    "#ax0.set_ylim(ax0ymin, ax0ymax)\n",
    "#setaxesfsize(ax0, axts)\n",
    "#ax0.legend(loc=3, fontsize=10, ncol=4)\n",
    "\n",
    "#ax0_dup.set_ylabel(\"AMOC Max (non-dim)\", size=axlbls)\n",
    "#ax0_dup.set_xlabel(\"Model Time (non-dim)\", size=axlbls)\n",
    "ax0_dup.set_title(\"Complex Model Stadial Age Test\", size=14, loc=\"center\")\n",
    "#ax0_dup.set_ylim(ax0ymin, ax0ymax)\n",
    "#setaxesfsize(ax0_dup, axts)\n",
    "#ax0_dup.legend(loc=3, fontsize=10, ncol=4)\n",
    "\n",
    "\n",
    "# Annotating my plots:\n",
    "flsize=20\n",
    "xannleft = 0.01\n",
    "yannleft = 0.86\n",
    "ax0.annotate('A', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "yannleft = 0.91\n",
    "ax1.annotate('B', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "ax2.annotate('C', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ") \n",
    "ax3.annotate('D', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ")\n",
    "\n",
    "## annotate ax0\n",
    "flax0size=10\n",
    "xb = 0.235\n",
    "yb = 0.8\n",
    "ax0.annotate('See B', \n",
    "             xy=(xb, yb), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "xc = 0.58\n",
    "yc = 0.8\n",
    "ax0.annotate('See C', \n",
    "             xy=(xc, yc), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold' \n",
    ")\n",
    "xd = 0.69\n",
    "yd = 0.8\n",
    "ax0.annotate('See D', \n",
    "             xy=(xd, yd), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold' \n",
    ") \n",
    "\n",
    "flsize=20\n",
    "xannleft = 0.01\n",
    "yannleft = 0.86\n",
    "ax0_dup.annotate('E', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "yannleft = 0.91\n",
    "ax1_dup.annotate('F', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "ax2_dup.annotate('G', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ") \n",
    "ax3_dup.annotate('H', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ") \n",
    "## annotate ax0_dup\n",
    "flax0size=10\n",
    "xb = 0.41\n",
    "yb = 0.732\n",
    "ax0_dup.annotate('See F', \n",
    "             xy=(xb, yb), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "xc = 0.49\n",
    "yc = 0.812\n",
    "ax0_dup.annotate('See G', \n",
    "             xy=(xc, yc), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold' \n",
    ")\n",
    "xd = 0.57\n",
    "yd = 0.888\n",
    "ax0_dup.annotate('See H', \n",
    "             xy=(xd, yd), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold' \n",
    ") \n",
    "# Save figure\n",
    "icepaperdir = \"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot = False\n",
    "if saveplot:\n",
    "    plt.savefig(\n",
    "        icepaperdir + \"Figure3_volc_perturb_complex_model.png\",\n",
    "        dpi=300,\n",
    "        bbox_inches=\"tight\",\n",
    "        facecolor=fig.get_facecolor(),\n",
    "        transparent=False,\n",
    "    )\n",
    "    plt.savefig(\n",
    "        icepaperdir + \"Figure3_volc_perturb_complex_model.pdf\",\n",
    "        dpi=300,\n",
    "        bbox_inches=\"tight\",\n",
    "        facecolor=fig.get_facecolor(),\n",
    "        transparent=False,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb6a4d6f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0322929",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (geosci)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
