{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1bb1e8a2-fce3-479f-842d-1a7b0db2e363",
   "metadata": {},
   "source": [
    "# Periodic Volcanic Forcing Experiments - Figure 2\n",
    "\n",
    "This notebook generates Figure 2, which demonstrates the impact of periodic equatorial volcanic forcing on Atlantic Meridional Overturning Circulation (AMOC) strength in a glacial climate state (230 ppm CO₂).\n",
    "\n",
    "\n",
    "\n",
    "## Experimental Design- Long-term climate system memory and response\n",
    "\n",
    "- Magnitude thresholds for triggering\n",
    "\n",
    "The analysis compares a control run (no volcanic forcing) with three periodic volcanic forcing experiments:- The repeatability of volcanic-triggered transitions\n",
    "\n",
    "This figure tests whether volcanic forcing can act as a pacemaker for millennial-scale climate oscillations similar to Dansgaard-Oeschger events. By applying forcing at regular intervals, we can assess:\n",
    "\n",
    "**Control Run:**\n",
    "\n",
    "- CO₂ concentration: 230 ppm (glacial levels)## Scientific Context\n",
    "\n",
    "- No volcanic forcing applied\n",
    "\n",
    "- Represents natural AMOC variability- Run length: ~7000+ model years\n",
    "\n",
    "- Forcing frequency: Every 1000 years\n",
    "\n",
    "**Volcanic Forcing Experiments:**- Volcanic forcing location: Equatorial\n",
    "\n",
    "All experiments apply equatorial volcanic eruptions with 1000-year periodicity:- Climate state: Glacial conditions (230 ppm CO₂)\n",
    "\n",
    "1. **High magnitude (230 Tg SO₄)**: Tests response to very large eruptions- Ocean component: POP2 (Parallel Ocean Program)\n",
    "\n",
    "2. **Medium magnitude (115 Tg SO₄)**: Tests response to moderate eruptions  - Model: CESM1.3 (iCESM)\n",
    "\n",
    "3. **Low magnitude (58 Tg SO₄)**: Tests response to smaller eruptions\n",
    "\n",
    "## Model Configuration\n",
    "\n",
    "## Key Questions Addressed\n",
    "\n",
    "- Red vertical lines: Timing of volcanic events\n",
    "\n",
    "1. Can periodic volcanic forcing trigger AMOC transitions in a glacial climate?- Colored line: Volcanic forcing experiment\n",
    "\n",
    "2. How does the magnitude of volcanic forcing affect AMOC response?- Blue line: Control run (no volcanism)\n",
    "\n",
    "3. Is there a threshold eruption size needed to induce transitions?\n",
    "\n",
    "4. How does the climate system respond over multiple forcing cycles?\n",
    "\n",
    "**Panel C**: AMOC response to 58 Tg eruptions (every 1000 years)\n",
    "\n",
    "## Figure Structure**Panel B**: AMOC response to 115 Tg eruptions (every 1000 years)\n",
    "\n",
    "**Panel A**: AMOC response to 230 Tg eruptions (every 1000 years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c4b04fec-df11-4321-86f8-090aac8f47bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "    <style>\n",
       "        .bk-notebook-logo {\n",
       "            display: block;\n",
       "            width: 20px;\n",
       "            height: 20px;\n",
       "            background-image: url(data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABQAAAAUCAYAAACNiR0NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAABx0RVh0U29mdHdhcmUAQWRvYmUgRmlyZXdvcmtzIENTNui8sowAAAOkSURBVDiNjZRtaJVlGMd/1/08zzln5zjP1LWcU9N0NkN8m2CYjpgQYQXqSs0I84OLIC0hkEKoPtiH3gmKoiJDU7QpLgoLjLIQCpEsNJ1vqUOdO7ppbuec5+V+rj4ctwzd8IIbbi6u+8f1539dt3A78eXC7QizUF7gyV1fD1Yqg4JWz84yffhm0qkFqBogB9rM8tZdtwVsPUhWhGcFJngGeWrPzHm5oaMmkfEg1usvLFyc8jLRqDOMru7AyC8saQr7GG7f5fvDeH7Ej8CM66nIF+8yngt6HWaKh7k49Soy9nXurCi1o3qUbS3zWfrYeQDTB/Qj6kX6Ybhw4B+bOYoLKCC9H3Nu/leUTZ1JdRWkkn2ldcCamzrcf47KKXdAJllSlxAOkRgyHsGC/zRday5Qld9DyoM4/q/rUoy/CXh3jzOu3bHUVZeU+DEn8FInkPBFlu3+nW3Nw0mk6vCDiWg8CeJaxEwuHS3+z5RgY+YBR6V1Z1nxSOfoaPa4LASWxxdNp+VWTk7+4vzaou8v8PN+xo+KY2xsw6une2frhw05CTYOmQvsEhjhWjn0bmXPjpE1+kplmmkP3suftwTubK9Vq22qKmrBhpY4jvd5afdRA3wGjFAgcnTK2s4hY0/GPNIb0nErGMCRxWOOX64Z8RAC4oCXdklmEvcL8o0BfkNK4lUg9HTl+oPlQxdNo3Mg4Nv175e/1LDGzZen30MEjRUtmXSfiTVu1kK8W4txyV6BMKlbgk3lMwYCiusNy9fVfvvwMxv8Ynl6vxoByANLTWplvuj/nF9m2+PDtt1eiHPBr1oIfhCChQMBw6Aw0UulqTKZdfVvfG7VcfIqLG9bcldL/+pdWTLxLUy8Qq38heUIjh4XlzZxzQm19lLFlr8vdQ97rjZVOLf8nclzckbcD4wxXMidpX30sFd37Fv/GtwwhzhxGVAprjbg0gCAEeIgwCZyTV2Z1REEW8O4py0wsjeloKoMr6iCY6dP92H6Vw/oTyICIthibxjm/DfN9lVz8IqtqKYLUXfoKVMVQVVJOElGjrnnUt9T9wbgp8AyYKaGlqingHZU/uG2NTZSVqwHQTWkx9hxjkpWDaCg6Ckj5qebgBVbT3V3NNXMSiWSDdGV3hrtzla7J+duwPOToIg42ChPQOQjspnSlp1V+Gjdged7+8UN5CRAV7a5EdFNwCjEaBR27b3W890TE7g24NAP/mMDXRWrGoFPQI9ls/MWO2dWFAar/xcOIImbbpA3zgAAAABJRU5ErkJggg==);\n",
       "        }\n",
       "    </style>\n",
       "    <div>\n",
       "        <a href=\"https://bokeh.org\" target=\"_blank\" class=\"bk-notebook-logo\"></a>\n",
       "        <span id=\"f90465bf-9e7b-4649-a886-f3521d68b0a5\">Loading BokehJS ...</span>\n",
       "    </div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": "'use strict';\n(function(root) {\n  function now() {\n    return new Date();\n  }\n\n  const force = true;\n\n  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n    root._bokeh_onload_callbacks = [];\n    root._bokeh_is_loading = undefined;\n  }\n\nconst JS_MIME_TYPE = 'application/javascript';\n  const HTML_MIME_TYPE = 'text/html';\n  const EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n  const CLASS_NAME = 'output_bokeh rendered_html';\n\n  /**\n   * Render data to the DOM node\n   */\n  function render(props, node) {\n    const script = document.createElement(\"script\");\n    node.appendChild(script);\n  }\n\n  /**\n   * Handle when an output is cleared or removed\n   */\n  function handleClearOutput(event, handle) {\n    function drop(id) {\n      const view = Bokeh.index.get_by_id(id)\n      if (view != null) {\n        view.model.document.clear()\n        Bokeh.index.delete(view)\n      }\n    }\n\n    const cell = handle.cell;\n\n    const id = cell.output_area._bokeh_element_id;\n    const server_id = cell.output_area._bokeh_server_id;\n\n    // Clean up Bokeh references\n    if (id != null) {\n      drop(id)\n    }\n\n    if (server_id !== undefined) {\n      // Clean up Bokeh references\n      const cmd_clean = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n      cell.notebook.kernel.execute(cmd_clean, {\n        iopub: {\n          output: function(msg) {\n            const id = msg.content.text.trim()\n            drop(id)\n          }\n        }\n      });\n      // Destroy server and session\n      const cmd_destroy = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n      cell.notebook.kernel.execute(cmd_destroy);\n    }\n  }\n\n  /**\n   * Handle when a new output is added\n   */\n  function handleAddOutput(event, handle) {\n    const output_area = handle.output_area;\n    const output = handle.output;\n\n    // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n    if ((output.output_type != \"display_data\") || (!Object.prototype.hasOwnProperty.call(output.data, EXEC_MIME_TYPE))) {\n      return\n    }\n\n    const toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n\n    if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n      toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n      // store reference to embed id on output_area\n      output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n    }\n    if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n      const bk_div = document.createElement(\"div\");\n      bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n      const script_attrs = bk_div.children[0].attributes;\n      for (let i = 0; i < script_attrs.length; i++) {\n        toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n        toinsert[toinsert.length - 1].firstChild.textContent = bk_div.children[0].textContent\n      }\n      // store reference to server id on output_area\n      output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n    }\n  }\n\n  function register_renderer(events, OutputArea) {\n\n    function append_mime(data, metadata, element) {\n      // create a DOM node to render to\n      const toinsert = this.create_output_subarea(\n        metadata,\n        CLASS_NAME,\n        EXEC_MIME_TYPE\n      );\n      this.keyboard_manager.register_events(toinsert);\n      // Render to node\n      const props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n      render(props, toinsert[toinsert.length - 1]);\n      element.append(toinsert);\n      return toinsert\n    }\n\n    /* Handle when an output is cleared or removed */\n    events.on('clear_output.CodeCell', handleClearOutput);\n    events.on('delete.Cell', handleClearOutput);\n\n    /* Handle when a new output is added */\n    events.on('output_added.OutputArea', handleAddOutput);\n\n    /**\n     * Register the mime type and append_mime function with output_area\n     */\n    OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n      /* Is output safe? */\n      safe: true,\n      /* Index of renderer in `output_area.display_order` */\n      index: 0\n    });\n  }\n\n  // register the mime type if in Jupyter Notebook environment and previously unregistered\n  if (root.Jupyter !== undefined) {\n    const events = require('base/js/events');\n    const OutputArea = require('notebook/js/outputarea').OutputArea;\n\n    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n      register_renderer(events, OutputArea);\n    }\n  }\n  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n    root._bokeh_timeout = Date.now() + 5000;\n    root._bokeh_failed_load = false;\n  }\n\n  const NB_LOAD_WARNING = {'data': {'text/html':\n     \"<div style='background-color: #fdd'>\\n\"+\n     \"<p>\\n\"+\n     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n     \"</p>\\n\"+\n     \"<ul>\\n\"+\n     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n     \"</ul>\\n\"+\n     \"<code>\\n\"+\n     \"from bokeh.resources import INLINE\\n\"+\n     \"output_notebook(resources=INLINE)\\n\"+\n     \"</code>\\n\"+\n     \"</div>\"}};\n\n  function display_loaded(error = null) {\n    const el = document.getElementById(\"f90465bf-9e7b-4649-a886-f3521d68b0a5\");\n    if (el != null) {\n      const html = (() => {\n        if (typeof root.Bokeh === \"undefined\") {\n          if (error == null) {\n            return \"BokehJS is loading ...\";\n          } else {\n            return \"BokehJS failed to load.\";\n          }\n        } else {\n          const prefix = `BokehJS ${root.Bokeh.version}`;\n          if (error == null) {\n            return `${prefix} successfully loaded.`;\n          } else {\n            return `${prefix} <b>encountered errors</b> while loading and may not function as expected.`;\n          }\n        }\n      })();\n      el.innerHTML = html;\n\n      if (error != null) {\n        const wrapper = document.createElement(\"div\");\n        wrapper.style.overflow = \"auto\";\n        wrapper.style.height = \"5em\";\n        wrapper.style.resize = \"vertical\";\n        const content = document.createElement(\"div\");\n        content.style.fontFamily = \"monospace\";\n        content.style.whiteSpace = \"pre-wrap\";\n        content.style.backgroundColor = \"rgb(255, 221, 221)\";\n        content.textContent = error.stack ?? error.toString();\n        wrapper.append(content);\n        el.append(wrapper);\n      }\n    } else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(() => display_loaded(error), 100);\n    }\n  }\n\n  function run_callbacks() {\n    try {\n      root._bokeh_onload_callbacks.forEach(function(callback) {\n        if (callback != null)\n          callback();\n      });\n    } finally {\n      delete root._bokeh_onload_callbacks\n    }\n    console.debug(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(css_urls, js_urls, callback) {\n    if (css_urls == null) css_urls = [];\n    if (js_urls == null) js_urls = [];\n\n    root._bokeh_onload_callbacks.push(callback);\n    if (root._bokeh_is_loading > 0) {\n      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls == null || js_urls.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n    root._bokeh_is_loading = css_urls.length + js_urls.length;\n\n    function on_load() {\n      root._bokeh_is_loading--;\n      if (root._bokeh_is_loading === 0) {\n        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n        run_callbacks()\n      }\n    }\n\n    function on_error(url) {\n      console.error(\"failed to load \" + url);\n    }\n\n    for (let i = 0; i < css_urls.length; i++) {\n      const url = css_urls[i];\n      const element = document.createElement(\"link\");\n      element.onload = on_load;\n      element.onerror = on_error.bind(null, url);\n      element.rel = \"stylesheet\";\n      element.type = \"text/css\";\n      element.href = url;\n      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n      document.body.appendChild(element);\n    }\n\n    for (let i = 0; i < js_urls.length; i++) {\n      const url = js_urls[i];\n      const element = document.createElement('script');\n      element.onload = on_load;\n      element.onerror = on_error.bind(null, url);\n      element.async = false;\n      element.src = url;\n      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.head.appendChild(element);\n    }\n  };\n\n  function inject_raw_css(css) {\n    const element = document.createElement(\"style\");\n    element.appendChild(document.createTextNode(css));\n    document.body.appendChild(element);\n  }\n\n  const js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.7.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.7.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.7.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.7.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.7.3.min.js\"];\n  const css_urls = [];\n\n  const inline_js = [    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\nfunction(Bokeh) {\n    }\n  ];\n\n  function run_inline_js() {\n    if (root.Bokeh !== undefined || force === true) {\n      try {\n            for (let i = 0; i < inline_js.length; i++) {\n      inline_js[i].call(root, root.Bokeh);\n    }\n\n      } catch (error) {display_loaded(error);throw error;\n      }if (force === true) {\n        display_loaded();\n      }} else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(run_inline_js, 100);\n    } else if (!root._bokeh_failed_load) {\n      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n      root._bokeh_failed_load = true;\n    } else if (force !== true) {\n      const cell = $(document.getElementById(\"f90465bf-9e7b-4649-a886-f3521d68b0a5\")).parents('.cell').data().cell;\n      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n    }\n  }\n\n  if (root._bokeh_is_loading === 0) {\n    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(css_urls, js_urls, function() {\n      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(window));",
      "application/vnd.bokehjs_load.v0+json": ""
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# numpy/scipy stuff\n",
    "import numpy as np\n",
    "from scipy.signal import argrelextrema\n",
    "\n",
    "# matplotlib stuff\n",
    "import matplotlib.pyplot as plt\n",
    "#plt.style.use('seaborn-talk')\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.ticker as ticker\n",
    "import matplotlib.patches as patches\n",
    "from matplotlib import colors\n",
    "import matplotlib.gridspec as gridspec\n",
    "from matplotlib.patches import FancyBboxPatch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.lines as mlines\n",
    "\n",
    "mpl.rcParams['figure.figsize'] = (16.0, 10.0)\n",
    "\n",
    "\n",
    "# data sci stuff\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "\n",
    "# system stuff\n",
    "import os\n",
    "import sys\n",
    "import glob\n",
    "import fnmatch\n",
    "\n",
    "# netcdf stuff\n",
    "from netCDF4 import Dataset\n",
    "import datetime\n",
    "import cftime\n",
    "import nc_time_axis\n",
    "\n",
    "# widgets stuff\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact,interactive\n",
    "\n",
    "# dask stuff\n",
    "from dask.diagnostics import Profiler, ResourceProfiler, CacheProfiler\n",
    "from dask.diagnostics import visualize\n",
    "\n",
    "# bokeh for dask extension\n",
    "import bokeh\n",
    "from bokeh.io import output_notebook\n",
    "output_notebook()\n",
    "\n",
    "from scipy.io import netcdf\n",
    "\n",
    "#import nsim\n",
    "# getting more fonts in matplotlib\n",
    "# sudo apt install msttcorefonts\n",
    "#\n",
    "# cd /usr/share/fonts/truetype/msttcorefonts\n",
    "# cp * /home/guido/anaconda3/lib/python3.8/site-packages/matplotlib/mpl-data/fonts/ttf\n",
    "#\n",
    "# rebuild fonts\n",
    "from matplotlib import rcParams\n",
    "#import matplotlib.font_manager as font_manager; font_manager._rebuild()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "da9b135a-74fe-479c-9132-a22edb559f94",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as feature\n",
    "\n",
    "# try geoviews later\n",
    "#import geoviews as gv\n",
    "#import geoviews.feature as gf\n",
    "#from bokeh.models import HoverTool\n",
    "from bokeh.models import CustomJS, HoverTool\n",
    "from bokeh.plotting import figure, show\n",
    "\n",
    "# widgets stuff\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact,interactive\n",
    "from ipywidgets import interact,FloatSlider,IntSlider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ac3dbc76-524f-495d-bcb0-6fbb4f7398c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def setaxesfsize(axl,fontsize):\n",
    "    ticklabelbot = axl.get_xticklabels()\n",
    "    ticklabelleft = axl.get_yticklabels()\n",
    "    for labelx in ticklabelbot:\n",
    "        labelx.set_fontsize(fontsize)\n",
    "#    i=0\n",
    "    for labely in ticklabelleft:\n",
    "#        if i != 0:\n",
    "            labely.set_fontsize(fontsize)\n",
    "#        else:\n",
    "#            labely.set_visible(False)\n",
    "#        i+=1\n",
    "#    # remove top label on y-axis\n",
    "#    labely.set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "08590df9-3604-4102-a565-e5cffd011103",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_amoc_max(mocarray):\n",
    "    latmin=20.\n",
    "    latmax=50.\n",
    "    zmin=50000.\n",
    "    zmax=200000.\n",
    "    region = 1\n",
    "    moc_max_em = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=0, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max_ei = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=1, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max_sm = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=2, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max = (moc_max_em +  moc_max_ei +  moc_max_sm).max(dim=('lat_aux_grid', 'moc_z'))\n",
    " \n",
    "    return moc_max\n",
    "\n",
    "def get_aabw_max(mocarray):\n",
    "    latmin=-80.\n",
    "    latmax=0.\n",
    "    zmin=300000.\n",
    "    zmax=400000.\n",
    "    region = 0\n",
    "    moc_max_em = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=0, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max_ei = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=1, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max_sm = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=2, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max = (moc_max_em +  moc_max_ei +  moc_max_sm).min(dim=('lat_aux_grid', 'moc_z'))\n",
    " \n",
    "    return moc_max\n",
    "\n",
    "def get_atl_aabw_max(mocarray):\n",
    "    latmin=-40.\n",
    "    latmax=0.\n",
    "    zmin=200000.\n",
    "    zmax=400000.\n",
    "    region = 1\n",
    "    moc_max_em = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=0, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max_ei = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=1, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max_sm = (mocarray\n",
    "               .sel(transport_reg=region, moc_comp=2, lat_aux_grid=slice(latmin, latmax), moc_z=slice(zmin, zmax))\n",
    "               .load())\n",
    "    moc_max = (moc_max_em +  moc_max_ei +  moc_max_sm).min(dim=('lat_aux_grid', 'moc_z'))\n",
    "\n",
    "    return moc_max\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a68cefbd-6e2d-4755-a357-e72fc14c6857",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.coding.times.CFDatetimeCoder at 0x76e23735bb90>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create the CFDatetimeCoder instance\n",
    "time_coder = xr.coders.CFDatetimeCoder(use_cftime=True)\n",
    "time_coder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "705d738d-99a1-42b0-b5af-f2e50a096935",
   "metadata": {},
   "outputs": [],
   "source": [
    "#run_230i = \"/cdata/climo/cesmi6gat31rblc230i/ocn/hist/ts/cesmi6gat31rblc230i_ANN_210912_998911_pop_MOC.nc\"\n",
    "run_230i = \"/cdata/climo/cesmi6gat31rblc230i/ocn/hist/ts/cesmi6gat31rblc230i_ANN_210912_998911_pop_decclimots.nc\"\n",
    "\n",
    "xrdata_230i = xr.open_dataset(run_230i, engine='netcdf4', chunks={'time': 20}, decode_times=time_coder, decode_timedelta=True)\n",
    "\n",
    "\n",
    "#run_1000p = \"/cdata/climo/cesmi6gat31rblc230iv1000p/ocn/hist/ts/cesmi6gat31rblc230iv1000p_ANN_260912_998911_pop_MOC.nc\"\n",
    "run_1000p = \"/cdata/climo/cesmi6gat31rblc230iv1000p/ocn/hist/ts/cesmi6gat31rblc230iv1000p_ANN_260912_998911_pop_decclimots.nc\"\n",
    "xrdata_1000p = xr.open_dataset(run_1000p, engine='netcdf4', chunks={'time': 20}, decode_times=time_coder, decode_timedelta=True)\n",
    "#run_1000ph = \"/cdata/climo/cesmi6gat31rblc230iv1000ph/ocn/hist/ts/cesmi6gat31rblc230iv1000ph_ANN_260912_998911_pop_MOC.nc\"\n",
    "run_1000ph = \"/cdata/climo/cesmi6gat31rblc230iv1000ph/ocn/hist/ts/cesmi6gat31rblc230iv1000ph_ANN_260912_998911_pop_decclimots.nc\"\n",
    "xrdata_1000ph = xr.open_dataset(run_1000ph, engine='netcdf4', chunks={'time': 20}, decode_times=time_coder, decode_timedelta=True)\n",
    "#run_1000pq = \"/cdata/climo/cesmi6gat31rblc230iv1000pq/ocn/hist/ts/cesmi6gat31rblc230iv1000pq_ANN_260912_998911_pop_MOC.nc\"\n",
    "run_1000pq = \"/cdata/climo/cesmi6gat31rblc230iv1000pq/ocn/hist/ts/cesmi6gat31rblc230iv1000pq_ANN_260912_998911_pop_decclimots.nc\"\n",
    "xrdata_1000pq = xr.open_dataset(run_1000pq, engine='netcdf4', chunks={'time': 20}, decode_times=time_coder, decode_timedelta=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "10f4b256-ff81-410d-8e27-ba3a5981dd07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# control\n",
    "amoc_max_230i = get_amoc_max(xrdata_230i.MOC)\n",
    "# volc experiments\n",
    "amoc_max_1000p = get_amoc_max(xrdata_1000p.MOC)\n",
    "amoc_max_1000ph = get_amoc_max(xrdata_1000ph.MOC)\n",
    "amoc_max_1000pq = get_amoc_max(xrdata_1000pq.MOC)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "4c1f2ced-72f2-4455-9f94-bbdd27529c2c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 301536 characters]",
      "text/plain": [
       "<Figure size 1600x900 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(16,9))\n",
    "\n",
    "\n",
    "lw=1.5\n",
    "lalp=0.7\n",
    "\n",
    "ax0=plt.subplot(311)\n",
    "ax1=plt.subplot(312)\n",
    "ax2=plt.subplot(313)\n",
    "#ax3=plt.subplot(514)\n",
    "#ax4=plt.subplot(515)\n",
    "ax = [ax0, ax1, ax2]\n",
    "\n",
    "cntcol = 'blue'\n",
    "cntalp = 0.8\n",
    "# volcanic run colors\n",
    "vcols = [\"red\", \"green\", \"purple\"]\n",
    "ecols = vcols\n",
    "for i in range(3):\n",
    "    amoc_max_230i.plot(ax=ax[i], color=cntcol, alpha=cntalp,linewidth=lw, label=\"Control: No Volcanism (CO$_2$ = 230 ppm)\")\n",
    "\n",
    "amoc_max_1000p.plot(ax=ax[0], color=vcols[0], alpha=lalp,linewidth=lw, label=\"Volcanism: 1000 Year Periodic EQ 230 Tg (CO$_2$ = 230 ppm)\")\n",
    "amoc_max_1000ph.plot(ax=ax[1], color=vcols[1], alpha=lalp,linewidth=lw, label=\"Volcanism: 1000 Year Periodic EQ 115 Tg (CO$_2$ = 230 ppm)\")\n",
    "amoc_max_1000pq.plot(ax=ax[2], color=vcols[2], alpha=lalp,linewidth=lw, label=\"Volcanism: 1000 Year Periodic EQ 58 Tg (CO$_2$ = 230 ppm)\")\n",
    "\n",
    "for i in range(3):\n",
    "\n",
    "    ax[i].set_ylabel(\"AMOC Max (Sv)\",size=16)\n",
    "    ax[i].set_xlabel(\"\",size=16)\n",
    "    #ax[i].grid()\n",
    "    ax[i].set_title(\"\", size=16, loc='center')\n",
    "\n",
    "\n",
    "ax[2].set_ylabel(\"AMOC Max (Sv)\",size=16)\n",
    "ax[2].set_xlabel(\"Model Time (Years)\",size=16)\n",
    "#ax[2].grid()\n",
    "#ax[2].set_title(\"Volcanic Test\", size=16, loc='center')\n",
    "ax[2].set_title(\"\", size=16, loc='center')\n",
    "\n",
    "# do background volcanic lines \n",
    "\n",
    "\n",
    "v2601 = cftime.num2date(2601*365, 'days since 0000-01-01 00:00:00', calendar='noleap')\n",
    "vgrid0 = [ cftime.num2date(yr*365, 'days since 0000-01-01 00:00:00', calendar='noleap') for yr in range(2701,10000,1000) ]\n",
    "# These are for fading lines\n",
    "vgrid1 = [ cftime.num2date((yr+20)*365, 'days since 0000-01-01 00:00:00', calendar='noleap') for yr in range(2701,10000,1000) ]\n",
    "vgrid2 = [ cftime.num2date((yr+40)*365, 'days since 0000-01-01 00:00:00', calendar='noleap') for yr in range(2701,10000,1000) ]\n",
    "for j in range(3):\n",
    "    ax[j].set_ylim(0,30)\n",
    "    setaxesfsize(ax[j],14)\n",
    "    \n",
    "    if j < 3:\n",
    "        volcol = ecols[0]\n",
    "    elif j == 3:\n",
    "        volcol = ecols[1]\n",
    "    else:\n",
    "        volcol = ecols[2]\n",
    "    \n",
    "    # volcanic pulse linewidths\n",
    "    vlw = 5.0\n",
    "    if j == 0 or j == 3 or j == 4:\n",
    "        alpfac = 1.0\n",
    "    elif j == 1:\n",
    "        alpfac = 0.5\n",
    "    elif j == 2:\n",
    "        alpfac = 0.25\n",
    "        \n",
    "    # try to move to back with zorder\n",
    "    zo = 0\n",
    "    # make thin black line at start\n",
    "    [ ax[j].vlines(vgrid0[i],0,30, color=vcols[j], alpha=0.7, linewidths=1.5, zorder=1) for i in range(0,8) ]\n",
    "\n",
    "\n",
    "    ax[j].legend(loc=4,fontsize=12,ncol=2)\n",
    "\n",
    "\n",
    "# add a,b,c,d\n",
    "lx = 0.006\n",
    "ly = 0.88\n",
    "ax[0].text(lx, ly, 'A',\n",
    "         transform=ax[0].transAxes, size=18, weight='bold')\n",
    "ax[1].text(lx, ly, 'B',\n",
    "         transform=ax[1].transAxes, size=18, weight='bold')\n",
    "ax[2].text(lx, ly, 'C',\n",
    "         transform=ax[2].transAxes, size=18, weight='bold')\n",
    "\n",
    "ax[0].grid()\n",
    "ax[1].grid()\n",
    "ax[2].grid()\n",
    "\n",
    "# volcanic forcing at 7950\n",
    "icepaperdir = \"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot=False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Figure2_volc_230_periodic_EQ.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=True)\n",
    "    plt.savefig(icepaperdir+\"Figure2_volc_230_periodic_EQ.pdf\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba07aab0-70fc-49e2-b5b7-3500fd7ee1ba",
   "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
}
