{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate averaged surface temperature and precipitation anomalies\n",
    "\n",
    "Goals:\n",
    "- Compute the averaged surface temperature anomaly and the same for precipitation anomaly over the period 2081-2100 compared to the period 1981-2000\n",
    "- Display delta-T delta-P diagram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example calculates the averaged temperature anomaly (using the **TG**) indicator) vs the precipitation anomaly (using the **PRCPTOT** indicator) for the period 2081-2100 compared to the reference 1981-2000 for SSP 585 and several climate models.\n",
    "\n",
    "We assume to have the **tas** and the **pr** variables in netCDF files in a `./data` folder.\n",
    "The data can be dowloaded using the [metalink](data/cmcc_gfdl_pr_and_tas.metalink) provided with this notebook.\n",
    "The data described in a `.metalink` file can be dowloaded with tools such as [aria2](https://aria2.github.io/) or a browser plugin such as [DownThemAll!](https://addons.mozilla.org/en-US/firefox/addon/downthemall/)\n",
    "If you wish to use a different dataset, you can use the [climate 4 impact portal](https://www.climate4impact.eu/c4i-frontend/) to search and select the data you wish to use and a metalink file to the [ESGF](https://esgf.llnl.gov/) data will be provided.\n",
    "\n",
    "\n",
    "The data is read using xarray and a plot of the a deltaT-deltaP diagram with data averaged over Europe is generated.\n",
    "\n",
    "The datasets that are expected for this notebook are tas and pr parameters (needed to calculate the TG and PRCPTOT indicators respectively) for several climate models, for the historical (1981-2000) and future (2081-2100) SSP 585 experiment and for several climate models and one member. Daily data is used. In C4I, you can find all of the data needed in the CMIP6 project, at the **esgf-data3.ceda.ac.uk** and **esgf.nci.org.au** mirrors.\n",
    "\n",
    "### Preparation of the modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: icclim in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (6.6.0)\n",
      "Requirement already satisfied: matplotlib in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (3.8.2)\n",
      "Requirement already satisfied: nc_time_axis in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (1.4.1)\n",
      "Requirement already satisfied: numpy>=1.16 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (1.26.2)\n",
      "Requirement already satisfied: xarray>=2022.6 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (2023.10.1)\n",
      "Requirement already satisfied: xclim<=0.47,>=0.45 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (0.47.0)\n",
      "Requirement already satisfied: cf_xarray>=0.7.4 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (0.8.6)\n",
      "Requirement already satisfied: cftime>=1.4.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (1.6.3)\n",
      "Requirement already satisfied: dask[array] in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (2023.12.1)\n",
      "Requirement already satisfied: netCDF4>=1.5.7 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (1.6.5)\n",
      "Requirement already satisfied: psutil in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (5.9.7)\n",
      "Requirement already satisfied: zarr in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (2.16.1)\n",
      "Requirement already satisfied: rechunker!=0.4,>=0.3 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (0.5.2)\n",
      "Requirement already satisfied: fsspec in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (2023.12.2)\n",
      "Requirement already satisfied: pandas>=1.3 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (2.1.4)\n",
      "Requirement already satisfied: dateparser in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (1.2.0)\n",
      "Requirement already satisfied: pint in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (0.19.2)\n",
      "Requirement already satisfied: jinja2 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from icclim) (3.1.2)\n",
      "Requirement already satisfied: contourpy>=1.0.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (1.2.0)\n",
      "Requirement already satisfied: cycler>=0.10 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (0.12.1)\n",
      "Requirement already satisfied: fonttools>=4.22.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (4.46.0)\n",
      "Requirement already satisfied: kiwisolver>=1.3.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (1.4.5)\n",
      "Requirement already satisfied: packaging>=20.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (23.2)\n",
      "Requirement already satisfied: pillow>=8 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (10.1.0)\n",
      "Requirement already satisfied: pyparsing>=2.3.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (3.1.1)\n",
      "Requirement already satisfied: python-dateutil>=2.7 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from matplotlib) (2.8.2)\n",
      "Requirement already satisfied: certifi in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from netCDF4>=1.5.7->icclim) (2023.11.17)\n",
      "Requirement already satisfied: pytz>=2020.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from pandas>=1.3->icclim) (2023.3.post1)\n",
      "Requirement already satisfied: tzdata>=2022.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from pandas>=1.3->icclim) (2023.3)\n",
      "Requirement already satisfied: six>=1.5 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)\n",
      "Requirement already satisfied: mypy-extensions in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from rechunker!=0.4,>=0.3->icclim) (1.0.0)\n",
      "Requirement already satisfied: boltons>=20.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (23.0.0)\n",
      "Requirement already satisfied: bottleneck>=1.3.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (1.3.7)\n",
      "Requirement already satisfied: Click>=8.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (8.1.7)\n",
      "Requirement already satisfied: jsonpickle in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (3.0.2)\n",
      "Requirement already satisfied: lmoments3>=1.0.5 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (1.0.6)\n",
      "Requirement already satisfied: numba in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (0.58.1)\n",
      "Requirement already satisfied: pyyaml in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (6.0.1)\n",
      "Requirement already satisfied: scikit-learn>=0.21.3 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (1.3.2)\n",
      "Requirement already satisfied: scipy>=1.2 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (1.11.4)\n",
      "Requirement already satisfied: statsmodels in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from xclim<=0.47,>=0.45->icclim) (0.14.1)\n",
      "Requirement already satisfied: cloudpickle>=1.5.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from dask[array]->icclim) (3.0.0)\n",
      "Requirement already satisfied: partd>=1.2.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from dask[array]->icclim) (1.4.1)\n",
      "Requirement already satisfied: toolz>=0.10.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from dask[array]->icclim) (0.12.0)\n",
      "Requirement already satisfied: importlib-metadata>=4.13.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from dask[array]->icclim) (7.0.0)\n",
      "Requirement already satisfied: asciitree in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from zarr->icclim) (0.3.3)\n",
      "Requirement already satisfied: fasteners in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from zarr->icclim) (0.17.3)\n",
      "Requirement already satisfied: numcodecs>=0.10.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from zarr->icclim) (0.12.1)\n",
      "Requirement already satisfied: regex!=2019.02.19,!=2021.8.27 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from dateparser->icclim) (2023.10.3)\n",
      "Requirement already satisfied: tzlocal in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from dateparser->icclim) (5.2)\n",
      "Requirement already satisfied: MarkupSafe>=2.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from jinja2->icclim) (2.1.3)\n",
      "Requirement already satisfied: zipp>=0.5 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from importlib-metadata>=4.13.0->dask[array]->icclim) (3.17.0)\n",
      "Requirement already satisfied: locket in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from partd>=1.2.0->dask[array]->icclim) (1.0.0)\n",
      "Requirement already satisfied: joblib>=1.1.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from scikit-learn>=0.21.3->xclim<=0.47,>=0.45->icclim) (1.3.2)\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from scikit-learn>=0.21.3->xclim<=0.47,>=0.45->icclim) (3.2.0)\n",
      "Requirement already satisfied: bokeh>=2.4.2 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from dask[array,diagnostics]->rechunker!=0.4,>=0.3->icclim) (3.3.2)\n",
      "Requirement already satisfied: llvmlite<0.42,>=0.41.0dev0 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from numba->xclim<=0.47,>=0.45->icclim) (0.41.1)\n",
      "Requirement already satisfied: patsy>=0.5.4 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from statsmodels->xclim<=0.47,>=0.45->icclim) (0.5.4)\n",
      "Requirement already satisfied: tornado>=5.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from bokeh>=2.4.2->dask[array,diagnostics]->rechunker!=0.4,>=0.3->icclim) (6.3.3)\n",
      "Requirement already satisfied: xyzservices>=2021.09.1 in /home/bzah/micromamba/envs/icclim-dev/lib/python3.11/site-packages (from bokeh>=2.4.2->dask[array,diagnostics]->rechunker!=0.4,>=0.3->icclim) (2023.10.1)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install icclim matplotlib nc_time_axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "python:  3.11.7 | packaged by conda-forge | (main, Dec 15 2023, 08:38:37) [GCC 12.3.0]\n",
      "numpy:  1.26.2\n",
      "xarray:  2023.10.1\n",
      "pandas:  2.1.4\n",
      "icclim:  6.6.0\n",
      "cftime:  1.6.3\n"
     ]
    }
   ],
   "source": [
    "import datetime\n",
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "import cftime\n",
    "import icclim\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "\n",
    "print(\"python: \", sys.version)\n",
    "print(\"numpy: \", np.__version__)\n",
    "print(\"xarray: \", xr.__version__)\n",
    "print(\"pandas: \", pd.__version__)\n",
    "print(\"icclim: \", icclim.__version__)\n",
    "print(\"cftime: \", cftime.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters Setup  \n",
    "The time period of interest as well as the reference period are defined here.\n",
    "A list of models is listed here as an example.\n",
    "Here we used Monthly data (Amon) but daily data could also be used.\n",
    "The corresponding datafiles must have been selected by the user, containing both the studied and referenced periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# studied period\n",
    "yearb = 2081\n",
    "yeare = 2100\n",
    "dt1 = datetime.datetime(yearb, 1, 1, tzinfo=datetime.timezone.utc)\n",
    "dt2 = datetime.datetime(yeare, 12, 31, tzinfo=datetime.timezone.utc)\n",
    "\n",
    "# reference period\n",
    "yearrefb = 1981\n",
    "yearrefe = 2000\n",
    "dtr1 = datetime.datetime(yearrefb, 1, 1, tzinfo=datetime.timezone.utc)\n",
    "dtr2 = datetime.datetime(yearrefe, 12, 31, tzinfo=datetime.timezone.utc)\n",
    "\n",
    "# studied domain\n",
    "# Western Europe\n",
    "minlat = 30\n",
    "maxlat = 56\n",
    "minlon = -30\n",
    "maxlon = 30\n",
    "\n",
    "models = [\"CMCC-ESM2\", \"GFDL-ESM4\"]\n",
    "out_f = {}\n",
    "out_hist_f = {}\n",
    "out_f_pr = {}\n",
    "out_hist_f_pr = {}\n",
    "for model in models:\n",
    "    out_f[model] = f\"data/dtdp_tg_icclim_{model}.nc\"\n",
    "    out_hist_f[model] = f\"data/dtdp_tg_icclim_{model}_hist.nc\"\n",
    "    out_f_pr[model] = f\"data/dtdp_pr_icclim_{model}.nc\"\n",
    "    out_hist_f_pr[model] = f\"data/dtdp_pr_icclim_{model}_hist.nc\"\n",
    "    data_dir = Path(\"data\")\n",
    "    filenames_hist = data_dir.glob(f\"tas_day_{model}_historical_*.nc\")\n",
    "    filenames_585 = data_dir.glob(f\"tas_day_{model}_ssp585_*.nc\")\n",
    "    filenames_hist_pr = data_dir.glob(f\"pr_day_{model}_historical_*.nc\")\n",
    "    filenames_585_pr = data_dir.glob(f\"pr_day_{model}_ssp585_*.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "icclim is then executed for both periods for each climate model separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for model in models:\n",
    "    icclim.index(\n",
    "        index_name=\"TG\",\n",
    "        in_files=[str(f) for f in filenames_585],\n",
    "        var_name=\"tas\",\n",
    "        slice_mode=\"year\",\n",
    "        time_range=[dt1, dt2],\n",
    "        out_file=out_f[model],\n",
    "        logs_verbosity=\"LOW\",\n",
    "    )\n",
    "    icclim.index(\n",
    "        index_name=\"TG\",\n",
    "        in_files=[str(f) for f in filenames_hist],\n",
    "        var_name=\"tas\",\n",
    "        slice_mode=\"year\",\n",
    "        time_range=[dtr1, dtr2],\n",
    "        out_file=out_hist_f[model],\n",
    "        logs_verbosity=\"LOW\",\n",
    "    )\n",
    "    icclim.index(\n",
    "        index_name=\"PRCPTOT\",\n",
    "        in_files=[str(f) for f in filenames_585_pr],\n",
    "        var_name=\"pr\",\n",
    "        slice_mode=\"year\",\n",
    "        time_range=[dt1, dt2],\n",
    "        out_file=out_f_pr[model],\n",
    "        logs_verbosity=\"LOW\",\n",
    "    )\n",
    "    icclim.index(\n",
    "        index_name=\"PRCPTOT\",\n",
    "        in_files=[str(f) for f in filenames_hist_pr],\n",
    "        var_name=\"pr\",\n",
    "        slice_mode=\"year\",\n",
    "        time_range=[dtr1, dtr2],\n",
    "        out_file=out_hist_f_pr[model],\n",
    "        logs_verbosity=\"LOW\",\n",
    "    )\n",
    "\n",
    "print(\"Processing finished.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data preparation\n",
    "\n",
    "Here all data is loaded in 2 separate variables, one containing all the historical periods for all the models, and the same for the future time period.\n",
    "\n",
    "Extract geographical domain and deal with longitudes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open datasets\n",
    "tg = []\n",
    "tg_hist = []\n",
    "pr = []\n",
    "pr_hist = []\n",
    "for model in models:\n",
    "    dsl = xr.open_dataset(out_f[model], decode_times=False)\n",
    "    dsl[\"time\"] = xr.decode_cf(dsl).time\n",
    "    dsl = dsl.assign_coords({\"model_id\": model})\n",
    "\n",
    "    dslp = xr.open_dataset(out_f_pr[model], decode_times=False)\n",
    "    dslp[\"time\"] = xr.decode_cf(dslp).time\n",
    "    dslp = dslp.assign_coords({\"model_id\": model})\n",
    "\n",
    "    dshl = xr.open_dataset(out_hist_f[model], decode_times=False)\n",
    "    dshl[\"time\"] = xr.decode_cf(dshl).time\n",
    "    dshl = dshl.assign_coords({\"model_id\": model})\n",
    "\n",
    "    dshlp = xr.open_dataset(out_hist_f_pr[model], decode_times=False)\n",
    "    dshlp[\"time\"] = xr.decode_cf(dshlp).time\n",
    "    dshlp = dshlp.assign_coords({\"model_id\": model})\n",
    "\n",
    "    # Reorder longitudes if needed and extract domain\n",
    "    if minlon > maxlon or minlon < 0:\n",
    "        dsl = dsl.assign_coords(lon=(((dsl.lon + 180) % 360) - 180)).roll(\n",
    "            lon=(dsl.dims[\"lon\"] // 2), roll_coords=True\n",
    "        )\n",
    "        dslp = dslp.assign_coords(lon=(((dslp.lon + 180) % 360) - 180)).roll(\n",
    "            lon=(dslp.dims[\"lon\"] // 2), roll_coords=True\n",
    "        )\n",
    "        dshl = dshl.assign_coords(lon=(((dshl.lon + 180) % 360) - 180)).roll(\n",
    "            lon=(dshl.dims[\"lon\"] // 2), roll_coords=True\n",
    "        )\n",
    "        dshlp = dshlp.assign_coords(lon=(((dshlp.lon + 180) % 360) - 180)).roll(\n",
    "            lon=(dshlp.dims[\"lon\"] // 2), roll_coords=True\n",
    "        )\n",
    "\n",
    "        if minlon >= 180:\n",
    "            minlon = minlon - 360\n",
    "        if maxlon >= 180:\n",
    "            maxlon = maxlon - 360\n",
    "\n",
    "    dsl = dsl.sel(lon=slice(minlon, maxlon), lat=slice(minlat, maxlat))\n",
    "    dslp = dslp.sel(lon=slice(minlon, maxlon), lat=slice(minlat, maxlat))\n",
    "    dshl = dshl.sel(lon=slice(minlon, maxlon), lat=slice(minlat, maxlat))\n",
    "    dshlp = dshlp.sel(lon=slice(minlon, maxlon), lat=slice(minlat, maxlat))\n",
    "\n",
    "    tg.append(dsl[\"TG\"])\n",
    "    pr.append(dslp[\"PRCPTOT\"])\n",
    "    tg_hist.append(dshl[\"TG\"])\n",
    "    pr_hist.append(dshlp[\"PRCPTOT\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform spatial average on the extracted geographical domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average different grids\n",
    "for ii in range(len(tg)):\n",
    "    tg[ii] = tg[ii].mean(dim=[\"lon\", \"lat\"])\n",
    "for ii in range(len(pr)):\n",
    "    pr[ii] = pr[ii].mean(dim=[\"lon\", \"lat\"])\n",
    "for ii in range(len(tg_hist)):\n",
    "    tg_hist[ii] = tg_hist[ii].mean(dim=[\"lon\", \"lat\"])\n",
    "for ii in range(len(pr_hist)):\n",
    "    pr_hist[ii] = pr_hist[ii].mean(dim=[\"lon\", \"lat\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Define a function to align all different calendar types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function to align different calendars using annual data\n",
    "def to_pandas_dt(da):\n",
    "    \"\"\"Takes an annual DataArray. Change the calendar to a pandas datetime.\"\"\"\n",
    "    val = da.copy()\n",
    "    # val.resample(time='Y').mean('time')\n",
    "    timev = []\n",
    "    years = [int(val) for val in da.time.dt.strftime(\"%Y\")]\n",
    "    for itime in range(val.sizes[\"time\"]):\n",
    "        timev.append(years[itime])\n",
    "\n",
    "    timevp = pd.to_datetime(timev, format=\"%Y\")\n",
    "    time1 = xr.DataArray(data=timevp, dims=[\"time\"])\n",
    "    time1.name = \"time\"\n",
    "    # We rename the time dimension and coordinate to time360 to make it clear it isn't\n",
    "    # the original time coordinate.\n",
    "    val = val.rename({\"time\": \"timepd\"})\n",
    "    time1 = time1.rename({\"time\": \"timepd\"})\n",
    "    val = val.assign_coords({\"timepd\": time1})\n",
    "\n",
    "    return val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Align calendars of all input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert all calendars to annual precision (we have configured icclim to output yearly data)\n",
    "ll = [to_pandas_dt(da) for da in tg]\n",
    "llp = [to_pandas_dt(da) for da in pr]\n",
    "ll_hist = [to_pandas_dt(da) for da in tg_hist]\n",
    "llp_hist = [to_pandas_dt(da) for da in pr_hist]\n",
    "\n",
    "# Concatenate all models into one\n",
    "full_tg = xr.concat(ll, \"model_id\", join=\"outer\")\n",
    "full_pr = xr.concat(llp, \"model_id\", join=\"outer\")\n",
    "full_tg_hist = xr.concat(ll_hist, \"model_id\", join=\"outer\")\n",
    "full_pr_hist = xr.concat(llp_hist, \"model_id\", join=\"outer\")\n",
    "full_tg_anomaly = full_tg.mean(dim=\"timepd\") - full_tg_hist.mean(dim=\"timepd\")\n",
    "full_pr_anomaly = (\n",
    "    (full_pr.mean(dim=\"timepd\") - full_pr_hist.mean(dim=\"timepd\"))\n",
    "    / full_pr.mean(dim=\"timepd\")\n",
    "    * 100.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Plot a multi-model scatter plot of the anomalies of temperature and precipitation of the future time period 2081-2100 SSP 585 compared to the historical period\n",
    "\n",
    "Temperature and Precipitation anomalies for SSP585 of the period 2080-2100 vs 1981-2000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Temperature Anomaly K')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn0AAAHgCAYAAADQY9vNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB64ElEQVR4nO3dd1gUV9sG8HvoHaUXFayAFXuPWLGCPWqsYEnU2GLsBbsxaizRqAmCXaNiSewFLFEj1lixgaJCFFTAAlLO94cf87rSdtdFJHv/rmuusGfOnPOcnc3u45wpkhBCgIiIiIj+03QKOgAiIiIiyn9M+oiIiIi0AJM+IiIiIi3ApI+IiIhICzDpIyIiItICTPqIiIiItACTPiIiIiItwKSPiIiISAsw6SMiIiLSAkz6CpHg4GBIkiQvRkZGcHBwQOPGjTFnzhw8efJE7bbDwsIgSRLCwsLksr179yIgIOCj2lNmyYmXl5dcR0dHB+bm5ihTpgy6dOmCbdu2ISMjQ63YAKBv375wdXVVKJs9ezZ27typdpsAsozN0tISXl5e2LNnj1Ixvb+tqakpXF1d4ePjg6CgIKSkpGTZxsvLC15eXh8V83/N0aNH4efnB3d3d5iamsLZ2Rm+vr44f/58tvUvXLiAZs2awczMDEWKFEHHjh1x7969LPViY2MxdOhQlCpVCsbGxnBxcYG/vz8ePHigUO/hw4cYMWIEGjVqhCJFikCSJAQHBysdf3p6OhYuXIiWLVuiWLFiMDExgYeHB8aNG4cXL15kqb9o0SJ07NgRJUuWhCRJuX4enjx5gr59+8LGxgYmJiaoW7cujhw5km3dw4cPo27dujAxMYGNjQ369u37Ud8xmqTKmA8cOID69evD2NgYlpaWaNeuHa5du5alXkpKCn788UdUrFgRpqamsLe3R6tWrXDq1KksdSdNmoS2bdvC2dkZkiShb9++KsUfEhKC7t27o0yZMjA2Noarqyu++uor3L59O9v6yu6L1NRUTJs2Da6urjA0NIS7uzuWLl2abZv37t1Dx44dUaRIEZiZmaF58+a4cOGCSuOgQkhQoREUFCQAiKCgIHH69Glx/PhxsW3bNjFixAhhaWkprKysxKFDh9RqOzQ0VAAQoaGhctmQIUOEuh+RhIQEcfr0aYXFwcFB1K9fP0t5Tho1aiRKlSol1zt8+LD49ddfRZs2bQQA0bBhQ/HixQu14uvTp49wcXFRKDM1NRV9+vRRq71MAETnzp3F6dOnxV9//SXWrVsn3NzchCRJ4s8//8wzJmNjY3m8R48eFWvWrBHdunUTurq6okKFCiI6Olphm2vXrolr1659VMz/NZ07dxaNGzcWy5cvF2FhYWLr1q2iTp06Qk9PTxw5ckSh7o0bN4S5ublo2LCh2LNnj9i+fbuoUKGCcHJyEk+ePJHrJScni7JlywobGxuxbNkyERoaKlasWCHs7e2Fs7OzSExMlOuGhoYKGxsb0axZM9G9e3f5/1llJSUlCXNzczFw4ECxdetWERoaKhYsWCCKFi0qypcvL16/fq1Q383NTVSrVk34+fkJW1tb0ahRo2zbTU5OFhUrVhTFihUT69evFwcPHhS+vr5CT09PhIWFKdQNCwsTenp6wtfXVxw8eFCsX79eODs7i4oVK4rk5GSlx5JflB3zzp07hSRJon379mLPnj1i48aNws3NTRQtWlTcuXNHoW6vXr2Ejo6OmDhxojhy5IjYunWrqF69utDT0xN///23Ql0TExNRp04d8fXXXwsDAwOVvzdq1aolfHx8xOrVq0VYWJhYt26d8PDwEGZmZuLq1asKdVXZF/379xeGhoZi3rx5IjQ0VIwbN05IkiRmzZqlUO/JkyfCyclJVKhQQWzfvl3s2bNHNGjQQJibm4ubN2+qNBYqXJj0FSKZSV94eHiWdffv3xfFixcX5ubmIjY2VuW2NZ30ZcfFxUW0adNG6fqNGjUSFSpUyHbd6tWrBQDRtWtXtWLJz6RvyJAhCmV37twRAESzZs3yjMnU1DTbdQcOHBD6+vqidu3aHxWfprx+/VpkZGQUdBjZ+vfff7OUJSUlCXt7e9G0aVOF8i5duggbGxuRkJAgl0VFRQl9fX0xZswYuezQoUMCgPjtt98Utt+4caMAIEJCQuSy9PR0+e/w8HCVk760tDQRFxeXpXzr1q0CgFi3bp1C+fv9VahQIccEaNmyZQKAOHXqlFyWmpoqypcvL2rVqqVQt2bNmqJ8+fIiNTVVLvvrr78EALF8+XKlx5JflB2zm5ubqFy5ssJnNSoqShgYGIgePXrIZcnJyUJXV1f07NlTYfvHjx8LAGLYsGE59q/O90Z2n9FHjx4JfX194e/vr1Cu7L64evWqkCRJzJ49W2H7AQMGCGNjYxEfHy+Xff/990JfX19ERUXJZQkJCcLGxkbt71QqHDi9+x9RokQJLFiwAElJSVi5cqXCunPnzsHHxwdWVlYwMjJC1apV8fvvv+faXt++fbFs2TIAilOWUVFRAIBly5bhiy++gJ2dHUxNTVGpUiXMmzcPqamp+TK+D/Xr1w+tW7fG1q1bcf/+fblcCIHly5fD09MTxsbGKFq0KDp37pztdN37JEnCq1evsGbNGnmsmVNGT58+xeDBg1G+fHmYmZnBzs4OTZo0wYkTJ5SKtXTp0rC1tVWIU1UtWrTAgAED8Pfff+P48eNyeXbTu9OmTUPt2rVhZWUFCwsLVKtWDYGBgRBCKNRLSUnBd999BwcHB5iYmOCLL77A+fPn4erqqjBdlXlawcGDB+Hn5wdbW1uYmJggJSUFd+7cQb9+/VC2bFmYmJjA2dkZ7dq1w5UrVxT6ypzu37hxI8aOHQtHR0eYmZmhXbt2+Pfff5GUlISBAwfCxsYGNjY26NevH16+fKnWe2VnZ5elzMzMDOXLl0d0dLRclpaWhj///BOdOnWChYWFXO7i4oLGjRtjx44dcpm+vj4AwNLSUqHdIkWKAACMjIzkMh2dj/ta1dXVhbW1dZbyWrVqAYDCGFTpb8eOHXBzc0PdunXlMj09PfTs2RNnz57Fo0ePAACPHj1CeHg4evXqBT09PbluvXr1UK5cOYX35UOpqamws7NDr169sqx78eIFjI2NMWrUKABARkYGZs6cCTc3NxgbG6NIkSKoXLkyFi9enOdYlBlzfHw8IiIi0KpVK4XTSFxcXFCxYkXs3LkT6enpcns6OjpZ9q+FhQV0dHQU9q+y/ecmu8+ok5MTihUrprB/VdkXO3fuhBAC/fr1U2i3X79+ePPmDfbv3y+X7dixA02aNIGLi4tcZmFhgY4dO+KPP/5AWlraR42PPl9M+v5DWrduDV1dXYWkIDQ0FPXr18eLFy+wYsUK7Nq1C56envjyyy9zPc9o8uTJ6Ny5MwDg9OnT8uLo6AgAuHv3Lnr06IF169bhzz//hL+/P3788UcMGjQoX8f4Ph8fHwghFJKvQYMGYcSIEWjWrBl27tyJ5cuX49q1a6hXrx7+/fffHNs6ffo0jI2N0bp1a3msy5cvBwA8e/YMADB16lTs2bMHQUFBKFWqFLy8vBTOgczJ8+fPER8fD1tb248eLwCF/ZudqKgoDBo0CL///jtCQkLQsWNHfPvtt5gxY4ZCvX79+mHRokXo168fdu3ahU6dOqFDhw7ZnjcGAH5+ftDX18e6deuwbds26Ovr4/Hjx7C2tsbcuXOxf/9+LFu2DHp6eqhduzYiIiKytDFhwgQ8efIEwcHBWLBgAcLCwtC9e3d06tQJlpaW2LRpE8aMGYN169ZhwoQJ6r1R2UhISMCFCxdQoUIFuezu3bt48+YNKleunKV+5cqVcefOHSQnJwMA6tevj+rVqyMgIADh4eF4+fIlLly4gAkTJqBatWpo1qyZxmLNydGjRwFAYQyquHr1ao5jBSCf53b16lWF8g/rZq7Pjr6+Pnr27Int27cjMTFRYd2mTZuQnJwsJyXz5s1DQEAAunfvjj179mDLli3w9/fP8fOnqrdv3wIADA0Ns6wzNDTE69evcffuXTnuwYMHY82aNdi5cycSExMRFRWFAQMGwNLSEgMGDNBITLm5d+8e7t+/r7B/VdkXV69eha2tLRwcHLLUe7+tN2/e4O7duzm2+ebNmzz/kUyFWMEeaCRV5Da9m8ne3l54eHjIr93d3UXVqlUVpgaEEKJt27bC0dFRnqb4mOnd9PR0kZqaKtauXSt0dXXFs2fPsq2nyeldIYTYt2+fACB++OEHIYQQp0+fFgDEggULFOpFR0cLY2Njhem6j5neTUtLE6mpqaJp06aiQ4cOCusAiMGDB4vU1FTx9u1bcePGDdGqVSsBQCxbtizXdnOb3hXi3flnAMQ333wjlzVq1CjHqS0h/rdvpk+fLqytreVprmvXrgkAYuzYsQr1N23aJAAovA+Zn7vevXvnGr8Q796bt2/firJly4qRI0fK5Zmfr3bt2inUHzFiRLbTZ+3btxdWVlZ59qesr776Sujp6Ylz587JZZlTZJs2bcpSf/bs2QKAePz4sVyWmJgo2rVrJwDIi5eXl8K02YfUmd7NzsOHD4W9vb2oUaOGwtTih3Kb6tTX1xeDBg3KUn7q1CkBQGzcuFEIIcSGDRsEgGzPtx04cKAwMDDINdZ//vlHABCrVq1SKK9Vq5aoXr26/Lpt27bC09Mz17aUkdOY09PThZWVVZYp/efPnwtzc/MsU90ZGRliypQpQkdHR96/JUqUEBcvXsy1f02cFpKamiq8vLyEhYWFePDggVyuyr5o3ry5cHNzy7Z9AwMDMXDgQCHEu2lkAGLOnDlZ6mWervD++0L/LTzS9x8j3pvCu3PnDm7evImvvvoKwLvprMyldevWiImJyfZojDIuXrwIHx8fWFtbQ1dXF/r6+ujduzfS09Nx69YtpdvJyMhQiCtzukUZ4oPpyj///BOSJKFnz54KbTo4OKBKlSpKHZXLyYoVK1CtWjUYGRlBT08P+vr6OHLkCG7cuJGl7vLly6Gvrw8DAwN4eHjg1KlTmD59OgYPHqx2/0DW8ebk6NGjaNasGSwtLeV9M2XKFMTHx8tX/B07dgwA0LVrV4VtO3furDCN9L5OnTplKUtLS8Ps2bNRvnx5GBgYQE9PDwYGBrh9+3a2703btm0VXnt4eAAA2rRpk6X82bNnak/xvm/y5MnYsGEDfvrpJ1SvXj3L+tyuIM9cl5qaii+//BKXLl3Cr7/+iuPHj2PNmjV49OgRmjdvjoSEBJXjUvaz/+zZM7Ru3RpCCGzZsuWjphaVGWtedXNrAwAqVaqE6tWrIygoSC67ceMGzp49Cz8/P7msVq1auHz5MgYPHowDBw5kOTL4sXR0dDBkyBAcOXIEM2bMwJMnT3Dnzh307NkTr1+/lutkmjVrFubPn4+AgACEhoZi165dcHNzQ/PmzXHx4kWV+xdCKOzfnKZMhRDw9/fHiRMnsHbtWhQvXjxLHWX3hSb2b17rqHBj0vcf8urVK8THx8PJyQkA5OnM0aNHQ19fX2HJTEDi4uJU7ufBgwdo2LAhHj16hMWLF+PEiRMIDw+XzwF88+aN0m1lThlmLk2bNlV628xz5N4frxAC9vb2WcZ75swZtcYKAAsXLsQ333yD2rVrY/v27Thz5gzCw8PRsmXLbMfatWtXhIeH49y5c4iIiEB8fDwmT56sVt/v+3C82Tl79ixatGgBAPj111/x119/ITw8HBMnTgTwv30THx8PALC3t1fYXk9PL9vzyQDIU/vvGzVqFCZPnoz27dvjjz/+wN9//43w8HBUqVIl2/fGyspK4bWBgUGu5ZnTq+qaNm0aZs6ciVmzZmHo0KEK6zLHmflevO/Zs2eQJEk+Zy8wMBD79u1DSEgI+vfvj4YNG6J3797Yv38/Lly4gEWLFqkcmzKf/efPn6N58+Z49OgRDh06hFKlSqncTyZra+scxwr8bx/k9b58uK+y4+fnh9OnT+PmzZsAgKCgIBgaGqJ79+5ynfHjx2P+/Pk4c+YMWrVqBWtrazRt2hTnzp1TfXA5mDJlCkaOHImZM2fC3t4eZcuWBQB5itnZ2RnAu6R0ypQpmDZtGiZPngwvLy/4+Phgz549KFKkiHweoirWrFmT5XvoQ0II9O/fH+vXr0dwcDB8fX0V1quyL3Lav69evcLbt2/lukWLFoUkSUp9Fui/J/t/0lOhtGfPHqSnp8sn9tvY2AB49+XasWPHbLdxc3NTuZ+dO3fi1atXCAkJUTgR+NKlSyq3FRAQoPBjbG5urvS2u3fvhiRJ+OKLLwC8G68kSThx4kSO5/GoY/369fDy8sIvv/yiUJ6UlJRtfVtbW9SoUUOtvnKze/duAMj1nmSbN2+Gvr4+/vzzT4WTzz+8/2Dmj8m///4r//AB747cZfdjAGT/r//169ejd+/emD17tkJ5XFycnDAVlGnTpiEgIAABAQHZnh9YunRpGBsbZ7noBACuXLmCMmXKyO/hpUuXoKuri2rVqinUK1WqFKytrXM9zy0neX32nz9/jmbNmiEyMhJHjhzJ9hwsVVSqVCnHsQJAxYoVFf575coVtG7dOkvdzPW56d69O0aNGoXg4GDMmjUL69atQ/v27VG0aFG5jp6eHkaNGoVRo0bhxYsXOHz4MCZMmABvb29ER0fDxMRE7bG+38fChQsxffp0REZGwsbGBo6OjvD29kbJkiVRrFgxAMDly5chhEDNmjUVttfX10eVKlXkI+OqaNeuHcLDw3Ncn5nwBQUFITAwED179sxSR5V9UalSJWzevBmxsbEK5/V9uH+NjY1RpkyZHD8LxsbGH/WPC/q88Ujff8SDBw8wevRoWFpayhdTuLm5oWzZsrh8+TJq1KiR7ZJbkpWZJH14xCbzx//9JEoIgV9//VXluF1dXRXiUTYJDQoKwr59+9C9e3eUKFECwLupQyEEHj16lO1YK1WqlGubhoaG2R6dkiQpS8L4zz//4PTp00qO8uMdOnQIv/32G+rVq4cGDRrkWE+SJOjp6UFXV1cue/PmDdatW6dQLzNR3rJli0L5tm3bVLpyL7v3Zs+ePfKVoAVlxowZCAgIwKRJkzB16tRs6+jp6aFdu3YICQlRSOAfPHiA0NBQhX8oOTk5IT09PcuP+K1btxAfHy8nD6rI7bOfmfDdu3cPBw8eRNWqVVVu/0MdOnTAzZs38ffff8tlaWlpWL9+PWrXri0fQXZ2dkatWrWwfv16hSnnM2fOICIiIsd/QL6vaNGiaN++PdauXYs///wTsbGxClO7HypSpAg6d+6MIUOG4NmzZ/JdAjTFzMwMlSpVgqOjIy5cuIAjR45g+PDh8vrMsZ85c0Zhu5SUFFy4cEGt/WttbZ3lOyiTEAIDBgxAUFAQVq5cmeWK20yq7AtfX19IkoQ1a9YotBEcHAxjY2O0bNlSLuvQoQOOHj2qcKVwUlISQkJC4OPjk+MpHlT4cc8WQlevXpXPEXny5AlOnDiBoKAg6OrqYseOHQpXia5cuRKtWrWCt7c3+vbtC2dnZzx79gw3btzAhQsXsHXr1hz7yUySfvjhB7Rq1Qq6urqoXLkymjdvDgMDA3Tv3h1jxoxBcnIyfvnlFzx//lzjY33z5o38RZx5VdnOnTvx559/olGjRlixYoVct379+hg4cCD69euHc+fO4YsvvoCpqSliYmJw8uRJVKpUCd98802u4w0LC8Mff/wBR0dHmJubw83NDW3btsWMGTMwdepUNGrUCBEREZg+fTpKliyp8VsbZGRkyONNSUnBgwcPsG/fPvz+++/w8PDI81Y7bdq0wcKFC9GjRw8MHDgQ8fHxmD9/fpbErEKFCujevTsWLFgAXV1dNGnSBNeuXcOCBQtgaWmp9Hljbdu2RXBwMNzd3VG5cmWcP38eP/74o1o/kpqyYMECTJkyBS1btkSbNm2y/JDXqVNH/nvatGmoWbMm2rZti3HjxiE5ORlTpkyBjY0NvvvuO7lev3798NNPP6FTp06YNGkS3NzccO/ePcyePRumpqb4+uuvFfrYtm0bAMhXQZ47dw5mZmYAIF8Vn5M3b97A29sbFy9exKJFi5CWlqYwBltbW5QuXVp+fe7cOTlJSkxMhBBC7r9mzZry0Xg/Pz8sW7YMXbp0wdy5c2FnZ4fly5cjIiIChw8fVojhhx9+QPPmzdGlSxcMHjwYT548wbhx41CxYsUcE5QP+fn5YcuWLRg6dCiKFSuW5Qrndu3aoWLFiqhRo4Z8S6NFixbBxcVFnobNibJjDgsLQ3h4OCpXrgwhBM6ePYsffvgBLVu2VDjK2qBBA9SsWRMBAQF4/fo1vvjiCyQkJGDp0qWIjIzM8o+mY8eO4enTpwDePUHl/v37cv+NGjXK80r9YcOGITAwEH5+fqhUqZLC/jU0NFRI8pXdFxUqVIC/vz+mTp0KXV1d1KxZEwcPHsSqVaswc+ZMhSnb0aNHY926dWjTpg2mT58OQ0NDzJ07F8nJyWo/hYkKiU9/7QipK/MqyszFwMBA2NnZiUaNGonZs2crPEHgfZcvXxZdu3YVdnZ2Ql9fXzg4OIgmTZqIFStWyHWyu3o3JSVF9O/fX9ja2gpJkgQAERkZKYQQ4o8//hBVqlQRRkZGwtnZWXz//ffy1bTvt/E+da7efX+8pqamolSpUqJz585i69atOV7FuHr1alG7dm1hamoqjI2NRenSpUXv3r0VrtzM7urdS5cuifr16wsTExMBQL4iMCUlRYwePVo4OzsLIyMjUa1aNbFz585s20A2N2dWVp8+fRTGa2xsLEqUKCHatWsnVq9eLVJSUrJ9jz68cnH16tXCzc1NGBoailKlSok5c+aIwMBAhf0nxLsb0o4aNUrY2dkJIyMjUadOHXH69GlhaWmpcOVtbleNP3/+XPj7+ws7OzthYmIiGjRoIE6cOJElrszP19atWxW2z6ntqVOnCgDi6dOnKryD/3tP3n8fP1w+dO7cOdG0aVNhYmIiLCwsRPv27bM8rUEIIW7fvi169eolXF1dhaGhoShRooT48ssvs30iiir9fygyMjLX7T+8UvTDz837y4dXDcfGxorevXsLKysreZ/n9BSfgwcPijp16ggjIyNhZWUlevfune1NhXOSnp4uihcvLgCIiRMnZlm/YMECUa9ePWFjYyMMDAxEiRIlhL+/v8INg3Oi7Jj/+usvUbt2bWFhYSEMDQ1FxYoVxfz588Xbt2+ztPnixQsxceJE4eHhIUxMTISdnZ3w8vISe/fuzVI3t89YTt9/73Nxcclx+w+/U4RQfl+8fftWTJ06VZQoUUIYGBiIcuXKiSVLlmQbw507d0T79u2FhYWFMDExEU2bNhXnz5/PM3Yq3CQhlLwkkIj+806dOoX69etjw4YN6NGjR0GHQ0REGsSkj0hLHTp0CKdPn0b16tVhbGyMy5cvY+7cubC0tMQ///yT5SkERERUuPGcPiItZWFhgYMHD2LRokVISkqCjY0NWrVqhTlz5jDhIyL6D+KRPiIiIiItwFu2EBEREWkBJn1EREREWoBJHxEREZEWYNJHREREpAWY9BERERFpASZ9RERERFqASR8RERGRFmDSR0RERKQFmPQRERERaQEmfURERERagEkfERERkRZg0kdERESkBZj0aaHg4GBIkiQvRkZGcHBwQOPGjTFnzhw8efKkoENUiqurK/r27Zvv/fTt2xeurq5qbbtx40YsWrRIY7F8uO/i4uLkdZs2bcIXX3wBe3t7GBoawsnJCe3atcOpU6eybWvz5s3w9PSEkZERnJycMGLECLx8+VLt2Pr27asQW+bi7u6ebf379+/Dz88PTk5OMDQ0hLOzMzp06KBQZ9GiRTmON6cY1N1X+e3996dixYpKbbNkyRLUqVMHNjY2MDQ0RIkSJdCtWzdcu3Yt2/pLly6Fu7s7DA0NUbJkSUybNg2pqalqxzxx4kRUrVoVVlZWMDIyQqlSpTBw4EDcv38/S93U1FRMmzYNrq6uMDQ0hLu7O5YuXZql3ogRI+T3wczMTO3YiEh1egUdABWcoKAguLu7IzU1FU+ePMHJkyfxww8/YP78+diyZQuaNWtW0CHmaseOHbCwsCjoMHK1ceNGXL16FSNGjNBouyEhIXB0dESRIkXksvj4eNSvXx/Dhw+HjY0NYmJisHDhQnzxxRc4cuQIGjVqJNfdsGEDevbsif79++Onn37CrVu3MHbsWFy/fh0HDx5UOy5jY2McPXo0S9mHrl69Ci8vL5QqVQrz589HsWLFEBMTgwMHDijU69atG+rUqYPffvsNgYGBasf1uXBwcMCOHTtgYmKiVP34+Hi0atUKVapUQdGiRXHv3j3MnTsXtWvXxvnz5+Hm5ibXnTVrFiZPnoxx48ahRYsWCA8Px6RJk/Do0SOsWrVKrXhfvHiB7t27w8PDA+bm5rh+/TpmzpyJ3bt349q1a7C2tpbrDh48GOvWrcOMGTNQs2ZNHDhwAMOHD0dSUhImTJgg1xs5ciS6deuGGTNm4NixY2rFRURqEqR1goKCBAARHh6eZd39+/dF8eLFhbm5uYiNjS2A6D4/ffr0ES4uLmpt26ZNG7W3zU7mvouMjFSq/osXL4S+vr7o1auXXJaWliYcHR1FixYtFOpu2LBBABB79+5VK7Y+ffoIU1PTPOtlZGQIT09P4enpKZKTk5Vqe+rUqQKAePr0aZ4xaPL91iRNxXb9+nUBQEyePFkui4uLE0ZGRmLgwIEKdWfNmiUkSRLXrl376H4z7d27VwAQgYGBctnVq1eFJEli9uzZCnUHDBggjI2NRXx8fJZ2lP28EJHmcHqXFJQoUQILFixAUlISVq5cqbBu9+7dqFu3LkxMTGBubo7mzZvj9OnTCnUCAgIgSRL++ecfdOnSBZaWlrCyssKoUaOQlpaGiIgItGzZEubm5nB1dcW8efMUtk9OTsZ3330HT09Pedu6deti165dWWL9cHo3LCwMkiRh06ZNmDhxIpycnGBhYYFmzZohIiJCc28SgGXLluGLL76AnZ0dTE1NUalSJcybN09hKs3Lywt79uzB/fv3FaYoPyVzc3MYGRlBT+9/B/XPnDmDmJgY9OvXT6Fuly5dYGZmhh07duRrTMePH8elS5cwYsQIGBoa5mtfwLvP1Pjx41GyZEkYGBjA2dkZQ4YMwYsXLxTqubq6om3btti/fz+qVasGY2NjuLu7Y/Xq1fkeoypsbW0BQGGf7t+/H8nJyVn2ab9+/SCEwM6dO/O1/507d0IIkW3/b968wf79+zXWPxGpj0kfZdG6dWvo6uri+PHjctnGjRvh6+sLCwsLbNq0CYGBgXj+/Dm8vLxw8uTJLG107doVVapUwfbt2zFgwAD89NNPGDlyJNq3b482bdpgx44daNKkCcaOHYuQkBB5u5SUFDx79gyjR4/Gzp07sWnTJjRo0AAdO3bE2rVrlYp/woQJuH//Pn777TesWrUKt2/fRrt27ZCeni7XyUwQAwIC1HqP7t69ix49emDdunX4888/4e/vjx9//BGDBg2S6yxfvhz169eHg4MDTp8+LS+Z0tPTkZaWlueSkZGhUmzp6elITU1FVFQUvvnmGwghMGTIEHn91atXAQCVK1dW2E5fXx/u7u7yenW8efMGDg4O0NXVRbFixTB06FA8e/ZMoU7m58rc3BytW7eGkZERzMzM0LZtW9y8eVPtvrMjhED79u0xf/589OrVC3v27MGoUaOwZs0aNGnSBCkpKQr1L1++jO+++w4jR47Erl27ULlyZfj7+yv8vwBAqf2WlpYGIYRGxpGeno6UlBTcvHkT/fv3h52dnUKClbnPKlWqpLCdo6MjbGxsPmqfAu/G++bNG1y8eBEjRoxAuXLl0LFjR4X+bW1t4eDgoLBd5mfsY/snIg0p0OOMVCBym97NZG9vLzw8PIQQQqSnpwsnJydRqVIlkZ6eLtdJSkoSdnZ2ol69enJZ5jTcggULFNrz9PQUAERISIhclpqaKmxtbUXHjh1zjCMtLU2kpqYKf39/UbVqVYV1Li4uok+fPvLr0NBQAUC0bt1aod7vv/8uAIjTp0/LZWFhYUJXV1dMmzYtx74z5TUtl56eLlJTU8XatWuFrq6uePbsmbwut+ldFxcXASDPZerUqfI2ykzvurm5yds6OjqKkydPKqyfNWuWACBiYmKybNuiRQtRrly5HNvOzcKFC8XChQvFwYMHxcGDB8XEiROFiYmJcHd3F0lJSXK9QYMGCQDCwsJC+Pv7i8OHD4t169YJFxcXYWNjIx4/fpylbXWnd/fv3y8AiHnz5inU27JliwAgVq1aJZe5uLgIIyMjcf/+fbnszZs3wsrKSgwaNEgui4yMVGq/ARChoaE5xqYKQ0NDuc1y5cqJ69evK6wfMGCAMDQ0zHbbcuXKZZnKV0VMTIzCmGrXri0ePXqkUKd58+bCzc0t2+0NDAyyTDsLweldooLACzkoW+K9IxQRERF4/PgxRowYAR2d/x0cNjMzQ6dOnbBy5Uq8fv1a4eT0tm3bKrTn4eGBy5cvo1WrVnKZnp4eypQpk+VKwK1bt2LRokW4fPkyXr16JZcbGRkpFbuPj4/C68yjDffv30edOnUAAI0aNUJaWppS7WXn4sWLmDp1Kv76668sR7Ju3bqF2rVr59nGH3/8keVIU3acnJxUim379u149eoVHjx4gBUrVqBVq1bYvXs3vLy8FOrlNNWs7hT0yJEjFV43b94cVatWRefOnfHrr7/K6zOPXNatWxe//fabXL9ixYqoWrUqli1bhpkzZ6oVw4cyLyr58CrvLl26wM/PD0eOHMGAAQPkck9PT5QoUUJ+bWRkhHLlyil8Rp2cnBAeHq5U/+9faJGTjIwMhaO5kiRBV1dXoc6pU6fw9u1b3L17Fz/99BMaN26MI0eOoEKFCgrb5eRjTiuwsbFBeHg4UlJScOPGDcybNw+NGzdGWFgYHB0d871/ItIcJn2UxatXrxAfHy9PFcXHxwOAwhd8JicnJ2RkZOD58+cKSZ+VlZVCPQMDA5iYmGRJ3AwMDJCYmCi/DgkJQdeuXdGlSxd8//33cHBwgJ6eHn755Relz616/4pCAPJ5Y2/evFFq+7w8ePAADRs2hJubGxYvXgxXV1cYGRnh7NmzGDJkiNL9lC9fXqnpv/cTbWVkJgK1atVC+/btUbVqVQwfPhyXL18G8L/3Jz4+Hvb29grbPnv2LMu++xgdOnSAqakpzpw5I5dl9u/t7a1Q19PTE46Ojrhw4YLG+o+Pj4eenp58HlomSZLg4OAgf7Y/jO19hoaGCvvUwMAAnp6eSvX/YfKWHT8/P6xZs0Z+3ahRI4SFhSnUqVatGgCgTp068PHxQZkyZTBhwgT5XFdra2skJydn+ccX8G6fVq9eXal4s6Onp4caNWoAAOrXr4+WLVuiZMmSmDt3LhYvXiz3f+nSpSzbvnr1Cm/fvtXoZ4qI1Mdz+iiLPXv2ID09XT4ylPlDGBMTk6Xu48ePoaOjg6JFi2qk7/Xr16NkyZLYsmUL2rdvjzp16qBGjRpKHRH7VHbu3IlXr14hJCQEPXv2RIMGDVCjRg0YGBio1E7p0qWhr6+f5zJ9+nS1Y9XT00O1atVw69YtuSwzmb9y5YpC3bS0NNy8eVPpe8gpSwihkLh+eC5hbnU/lrW1NdLS0vD06dMs/cTGxsLGxkblNqOiopTab/r6+krdkiQgIADh4eHy8uEFVB8yNzeHu7u7Uvs0NjYWcXFxGt2nxYoVg5OTU5b+nz59itjYWIW6mfFo+jNFROrhkT5S8ODBA4wePRqWlpbyRQlubm5wdnbGxo0bMXr0aHmq5tWrV9i+fbt8Ra8mSJIEAwMDhemg2NjYbK/eLSiZsb1/5akQAr/++muWuh8eJXpffk3vvi85ORlnzpxBmTJl5LLatWvD0dERwcHB+PLLL+Xybdu24eXLlwon6H+sbdu24fXr1/K0OgC0atUKJiYm2Ldvn8KU8IULFxAbG6tQ92M1bdoU8+bNw/r16xX6ypwCb9q0qcptanp619XVVaUbSsfFxeHKlSuoX7++XNayZUsYGRkhODhY4dSCzJt5t2/fXun283Lnzh08fPhQ4TQKX19fTJo0CWvWrMHYsWMV+jc2NkbLli011j8RqY9Jnxa7evWqfJXhkydPcOLECQQFBUFXVxc7duyQp8R0dHQwb948fPXVV2jbti0GDRqElJQU/Pjjj3jx4gXmzp2rsZjatm2LkJAQDB48GJ07d0Z0dDRmzJgBR0dH3L59W2P9HDt2DE2bNsWUKVMwZcoUlbZt3rw5DAwM0L17d4wZMwbJycn45Zdf8Pz58yx1K1WqhJCQEPzyyy+oXr06dHR05KmyD6+0/Fj16tWDj48PPDw8YGlpiaioKPzyyy+4e/euwm1YdHV1MW/ePPTq1QuDBg1C9+7dcfv2bYwZMwbNmzfP8gMtSVK2U47vu3//Pnr06IFu3bqhTJkykCQJx44dw6JFi1ChQgX0799frlukSBFMnz4do0ePRt++fdG9e3fExsZi8uTJKFGiBAYPHqyx96R58+bw9vbG2LFjkZiYiPr16+Off/7B1KlTUbVqVfTq1UvlNg0MDOR9mJ8SEhLQvHlz9OjRA2XLloWxsTFu3bqFxYsXIyUlBVOnTpXrWllZYdKkSZg8eTKsrKzkmzMHBASgf//+KF++vFw3KioKJUuWRJ8+fRAcHJxj///88w9GjhyJzp07o1SpUtDR0cGVK1fw008/wdraGqNHj5brVqhQAf7+/pg6dSp0dXVRs2ZNHDx4EKtWrcLMmTM5vUv0mWDSp8Uyb/lgYGCAIkWKwMPDA2PHjkX//v2znAPVo0cPmJqaYs6cOfjyyy+hq6uLOnXqIDQ0FPXq1dNoTE+ePMGKFSuwevVqlCpVCuPGjcPDhw8xbdo0jfUjhEB6errKt0MBAHd3d2zfvh2TJk1Cx44dYW1tjR49emDUqFEKF6oAwPDhw3Ht2jVMmDABCQkJEEJo7DYeH6pXrx42b96MqKgovHr1CjY2Nqhbty5++umnLPuoZ8+e0NXVxdy5cxEcHAwrKyv07t0bs2bNUqiX+Vi27M7nfJ+FhQXs7e2xcOFC/Pvvv0hPT4eLiwuGDRuGCRMmwNTUVKH+d999B0tLSyxevBibNm2Cubk5WrZsiblz52o0QZAkCTt37kRAQACCgoIwa9Ys2NjYoFevXpg9e/YnuU+guoyMjFClShWsWrUK0dHRSE5OhoODA7y8vLB9+3aFRA5498g0c3NzLFu2DPPnz4eDgwPGjRuHiRMnKtRTdp/a29vDyckJCxYsQExMDNLS0lCsWDG0bdsWEyZMQPHixRXqL1++HM7Ozli6dCliY2Ph6uqKxYsX49tvv9XAu0FEmiCJ/PoFIiKNCw4ORr9+/XDnzh24uLgo3CA3P+zduxdt27bF5cuXNX5kUhmZyfn06dMxY8YMPH36VK3z8D4Hffv2RVhYGO7cuZPtFbqfyvLlyzFmzBjcvXs3y4U8n0Lm1cr+/v7Yvn37Rz3vmYhUwws5iAqhMmXKQF9fH3FxcfnaT2hoKLp161YgCR8ALF68GPr6+pgxY0aB9K9p9+/fh76+PqpUqVJgMYSGhmLYsGEFkvABwKhRo6Cvr6/0zdaJSHN4pI+oEImPj0dkZKT82tPTM9+P9hWkJ0+e4MGDB/LrwjzeqKgoOUk3NjZWuMeeNomOjsa///4L4N35pVWrVi3giIi0B5M+IiIiIi3A6V0iIiIiLcCkj4iIiEgLFM6TYz6hjIwMPH78GObm5nx+JBER5UoIgaSkJDg5OWn06TIfSk5Oxtu3bz+6HQMDA6Wfa06FH5O+PDx+/DjL/aiIiIhyEx0djWLFiuVL28nJySjpYobYJ+kf3ZaFhQUcHR2ho6ODIUOGYMiQIRqIkD5XTPryYG5uDuDd/8AWFhYFHA0REX3OEhMTUbx4cfm3Iz+8ffsWsU/Scf+8KyzM1T+amJiUAZfqUfx90yJM+vKQOaVrYWHB/ymIiEgpn+J0IAtzHViYF8xNvqlw4oUcRERERFqASR8RERGRFmDSR0RERKQFmPQRERERaQEmfURERERagEkfERERkRZg0kdERESkBZj0EREREWkBJn1EREREWoBP5CAiIq2WnpYOIQQAQE+fP4v038VPNxERaaWMjAzo6OjgxpnbCN9/EebW5mjRuxFMi5hCV5cTYfTfw6SPiIi0kgQJs79ajNBNJ+Wy4EmbMG3nWHg2rgBdPT7Xlv5b+E8ZIiLSOmmp6fhr11mFhA8AUt68xQL/5ZB0pAKKjCj/MOkjIiKto6evi792nM123dOH8bhzIfITR0SU/5j0ERGRVtI3zPkMJ31D/U8YCdGnwaSPiIi0TnpaOpr39sp2nWuF4ihZqcSnDYjoE2DSR0REWkdXTxeVGnqgd0BXhQs27ErYYNKWkUhLTS/A6IjyB6/eJSIirdVrShe0+7oFwvdfQhE7C1RvXgUZGQJ6+rxyl/57mPQREZFWK2JnicbdG0BHR4KOrg50mO/RfxSTPiIi0no8skfagOf0EREREWkBJn1EREREWoBJHxEREZEWYNJHREREpAWY9BERERFpASZ9RERERFqASR8RERGRFmDSR0RERKQFmPQRERERaQEmfURERERagEkfERERkRZg0kdERESkBZj0EREREWkBJn1EREREWoBJHxEREZEWYNJHREREpAWY9BERERFpASZ9RERERFqASR8RERGRFmDSR0RERKQFmPQRERERaQEmfURERERagEkfERERkRZg0kdERESkBZj0EREREWkBJn1EREREWoBJHxEREZEWYNJHREREpAWY9BERERFpASZ9RERERFqASR8RERGRFigUSV9YWBgkScp2CQ8Pz3ab1NRUjB07FpUqVYKpqSmcnJzQu3dvPH78+BNHT0RERFTwCkXSV69ePcTExCgs/fv3h6urK2rUqJHtNq9fv8aFCxcwefJkXLhwASEhIbh16xZ8fHw+cfREREREBU+voANQhoGBARwcHOTXqamp2L17N4YOHQpJkrLdxtLSEocOHVIoW7p0KWrVqoUHDx6gRIkS+RozERER0eekUCR9H9q9ezfi4uLQt29flbZLSEiAJEkoUqRIjnVSUlKQkpIiv05MTFQzSiIiIqLPR6GY3v1QYGAgvL29Ubx4caW3SU5Oxrhx49CjRw9YWFjkWG/OnDmwtLSUF1X6ICIiIvpcFWjSFxAQkOMFGpnLuXPnFLZ5+PAhDhw4AH9/f6X7SU1NRbdu3ZCRkYHly5fnWnf8+PFISEiQl+joaLXGRkRERPQ5KdDp3aFDh6Jbt2651nF1dVV4HRQUBGtra6UvyEhNTUXXrl0RGRmJo0eP5nqUDwAMDQ1haGioVNtEREREhUWBJn02NjawsbFRur4QAkFBQejduzf09fXzrJ+Z8N2+fRuhoaGwtrb+mHCJiIiICq1CdU7f0aNHERkZmePUrru7O3bs2AEASEtLQ+fOnXHu3Dls2LAB6enpiI2NRWxsLN6+ffspwyYiIiIqcIXq6t3AwEDUq1cPHh4e2a6PiIhAQkICgHfn/u3evRsA4OnpqVAvNDQUXl5e+RkqERER0WelUCV9GzduzHW9EEL+29XVVeE1ERERkTYrVNO7RERERKQeJn1EREREWoBJHxEREZEWYNJHREREpAWY9BERERFpASZ9RERERFqASR8RERGRFmDSR0RERKQFmPQRERERaQEmfURERERagEkfERERkRZg0kdERESkBZj0EREREWkBJn1EREREWoBJHxEREZEWYNJHREREpAWY9BERERFpASZ9RERERFqASR8RERGRFmDSR0RERKQFmPQRERERaQEmfURERERagEkfERERkRZg0kdERESkBZj0ERERkdJiY2MxfPhwlClTBkZGRrC3t0eDBg2wYsUKvH79GgDg6uoKSZIUlmLFisltvL/e2NgYrq6u6Nq1K44eParQV1RUFCRJwqVLl5SOz8vLK0vfkiTh66+/luuEhoaicePGsLKygomJCcqWLYs+ffogLS0NABAWFgZJklC0aFEkJycrtH/27Fm5zezcuXMH5ubmKFKkiNIxfypM+oiIiEgp9+7dQ9WqVXHw4EHMnj0bFy9exOHDhzFy5Ej88ccfOHz4sFx3+vTpiImJkZeLFy8qtJW5PiIiAmvXrkWRIkXQrFkzzJo166PjHDBggELfMTExmDdvHgDg2rVraNWqFWrWrInjx4/jypUrWLp0KfT19ZGRkaHQjrm5OXbs2KFQtnr1apQoUSLbflNTU9G9e3c0bNjwo8eQH/QKOgAiIiIqHAYPHgw9PT2cO3cOpqamcnmlSpXQqVMnCCHkMnNzczg4OOTY1vvrS5QogS+++AKOjo6YMmUKOnfuDDc3N7XjNDExybHvQ4cOwdHRUU4CAaB06dJo2bJllrp9+vTB6tWr0b17dwDAmzdvsHnzZgwbNgwzZszIUn/SpElwd3dH06ZNcerUKbXjzy880kdERKTFEhMTFZaUlJRs68XHx+PgwYMYMmSIQsL3vpymPJU1fPhwCCGwa9euj2onNw4ODoiJicHx48fzrNurVy+cOHECDx48AABs374drq6uqFatWpa6R48exdatW7Fs2TKNx6wpTPqIiIi0WPHixWFpaSkvc+bMybbenTt3IITIcgTOxsYGZmZmMDMzw9ixY+XysWPHyuVmZmZYsmRJnrFYWVnBzs4OUVFRHzWm5cuXK/RtZmaGNWvWAAC6dOmC7t27o1GjRnB0dESHDh3w888/IzExMUs7dnZ2aNWqFYKDgwG8m9r18/PLUi8+Ph59+/ZFcHAwLCwsPir2/MSkj4iISItFR0cjISFBXsaPH59r/Q+P5p09exaXLl1ChQoVFI4Sfv/997h06ZK89O7dW6l4hBBKHTHcsGGDQlJ34sQJed1XX32l0PelS5fQoUMHAICuri6CgoLw8OFDzJs3D05OTpg1axYqVKiAmJiYLP34+fkhODgY9+7dw+nTp/HVV19lqTNgwAD06NEDX3zxhVJjLCg8p4+IiEiLWVhYKHV0qkyZMpAkCTdv3lQoL1WqFADA2NhYodzGxgZlypRRKZb4+Hg8ffoUJUuWzLOuj48PateuLb92dnaW/7a0tMyzb2dnZ/Tq1Qu9evXCzJkzUa5cOaxYsQLTpk1TqNe6dWsMGjQI/v7+aNeuHaytrbO0dfToUezevRvz588H8C5xzcjIgJ6eHlatWpXt0cGCwKSPiIiI8mRtbY3mzZvj559/xrfffpvjeX0fY/HixdDR0UH79u3zrGtubg5zc3ON9Fu0aFE4Ojri1atXWdbp6uqiV69emDdvHvbt25ft9qdPn0Z6err8eteuXfjhhx9w6tQphWS0oDHpIyIiIqUsX74c9evXR40aNRAQEIDKlStDR0cH4eHhuHnzJqpXr650W0lJSYiNjUVqaioiIyOxfv16/Pbbb5gzZ06Wo3QRERFZti9fvjwMDAyybfv169eIjY1VKDM0NETRokWxcuVKebq3dOnSSE5Oxtq1a3Ht2jUsXbo02/ZmzJiB77//PtujfADg4eGh8PrcuXPQ0dFBxYoVcxx/QWDSR0REREopXbo0Ll68iNmzZ2P8+PF4+PAhDA0NUb58eYwePRqDBw9Wuq0pU6ZgypQpMDAwgIODA+rUqYMjR46gcePGWep269YtS1lkZCRcXV2zbfvXX3/Fr7/+qlDm7e2N/fv3o1atWjh58iS+/vprPH78GGZmZqhQoQJ27tyJRo0aZduegYEBbGxslB7b50oS799Uh7JITEyEpaUlEhISPusrcoiIqOB9it+MzD6e3yoFC3Nd9dtJSkfRcvf4+6ZFePUuERERkRZg0kdERESkBT7qnL64uDj8/fffSE9PR82aNeHo6KipuIiIiIhIg9RO+rZv3w5/f3+UK1cOqampiIiIwLJly9CvXz9NxkdEREREGqD09O7Lly8VXk+bNg1nz57F2bNncfHiRWzduhUTJ07UeIBERERE9PGUTvqqV6+u8ABkPT09PHnyRH7977//5ni/HCIiIiIqWEpP7x44cACDBw9GcHAwli1bhsWLF+PLL79Eeno60tLSoKOjIz+QmIiIiIg+L0onfa6urti7dy82btyIRo0aYfjw4bhz5w7u3LmD9PR0uLu7w8jIKD9jJSIiIiI1qXzLlh49esjn8Xl5eSEjIwOenp5M+IiIiIg+Yypdvbtv3z5cv34dVapUQWBgIMLCwtCjRw+0bt0a06dPh7GxcX7FSUREREQfQekjfWPGjEHfvn0RHh6OQYMGYcaMGfDy8sLFixdhaGgIT09P7Nu3Lz9jJSIiIiI1Kf3sXRsbGxw4cADVq1fHs2fPUKdOHdy6dUtef+3aNQwaNAgnT57Mt2ALAp+9S0REyuKzd+lzpvSRPhMTE0RGRgIAoqOjs5zDV6FChf9cwkdERET0X6F00jdnzhz07t0bTk5OaNSoEWbMmJGfcRERERGRBil9IcdXX32Fli1b4t69eyhbtiyKFCmSj2ERERERkSapdPWutbU1rK2t8ysWIiIiIsonKt+nryCEhYVBkqRsl/DwcKXaGDRoECRJwqJFi/I3WCIiIqLPkEpH+gpKvXr1EBMTo1A2efJkHD58GDVq1Mhz+507d+Lvv/+Gk5NTfoVIRERE9FkrFEmfgYEBHBwc5NepqanYvXs3hg4dCkmSct320aNHGDp0KA4cOIA2bdrk2VdKSgpSUlLk14mJieoHTkRERPSZUHl699WrV/kRh0p2796NuLg49O3bN9d6GRkZ6NWrF77//ntUqFBBqbbnzJkDS0tLeSlevLgGIiYiIiIqWConffb29vDz8yvQe/IFBgbC29s7z4Tshx9+gJ6eHoYNG6Z02+PHj0dCQoK8REdHf2y4RERERAVO5aRv06ZNSEhIQNOmTVGuXDnMnTsXjx8/VqvzgICAHC/QyFzOnTunsM3Dhw9x4MAB+Pv759r2+fPnsXjxYgQHB+c5Bfw+Q0NDWFhYKCxEREREhZ3Sj2H7UHx8PNauXYvg4GBcv34d3t7e8PPzg4+PD/T0lDtVMC4uDnFxcbnWcXV1VXj6x4wZM7B06VI8evQI+vr6OW63aNEijBo1Cjo6/8tr09PToaOjg+LFiyMqKkqpGPkYNiIiUhYfw0afM7WTvvctXboU33//Pd6+fQsbGxt8/fXXGDduHExMTDQRo0wIgdKlS6Njx46YP39+rnXj4+OzXPHr7e2NXr16oV+/fnBzc1OqTyZ9RESkLCZ99DlT++rd2NhYrF27FkFBQXjw4AE6d+4Mf39/PH78GHPnzsWZM2dw8OBBTcaKo0ePIjIyMsepXXd3d8yZMwcdOnTI9kbS+vr6cHBwUDrhIyIiIvqvUDnpCwkJQVBQEA4cOIDy5ctjyJAh6Nmzp8Jj2Tw9PVG1alVNxgng3QUc9erVg4eHR7brIyIikJCQoPF+iYiIiAo7lZO+fv36oVu3bvjrr79Qs2bNbOuUKlUKEydO/OjgPrRx48Zc1+c1U63seXxERERE/zUqJ30xMTF5nqtnbGyMqVOnqh0UEREREWmWUknfh0+lyO0pFTwZlIiIiOjzo1TSV6RIkTzvdSeEgCRJSE9P10hgRERERKQ5SiV9oaGh+R0HEREREeUjpZK+Ro0a5XccRERERJSP1L5P3+vXr/HgwQO8fftWobxy5cofHRQRERERaZbKSd/Tp0/Rr18/7Nu3L9v1PKePiIiI6POjk3cVRSNGjMDz589x5swZGBsbY//+/VizZg3Kli2L3bt350eMRERERPSRVD7Sd/ToUezatQs1a9aEjo4OXFxc0Lx5c1hYWGDOnDlo06ZNfsRJRERERB9B5SN9r169gp2dHQDAysoKT58+BQBUqlQJFy5c0Gx0RERERKQRKid9bm5uiIiIAPDuGbsrV67Eo0ePsGLFCjg6Omo8QCIiIiL6eCpP744YMQIxMTEAgKlTp8Lb2xsbNmyAgYEBgoODNR0fEREREWmAyknfV199Jf9dtWpVREVF4ebNmyhRogRsbGw0GhwRERERaYba9+nLZGJigmrVqmkiFiIiIiLKJyonfUIIbNu2DaGhoXjy5AkyMjIU1oeEhGgsOCIiIiLSDJWTvuHDh2PVqlVo3Lgx7O3tIUlSfsRFRERERBqkctK3fv16hISEoHXr1vkRDxERERHlA5Vv2WJpaYlSpUrlRyxERERElE9UTvoCAgIwbdo0vHnzJj/iISIiIqJ8oPL0bpcuXbBp0ybY2dnB1dUV+vr6Cuv5VA4iIiKiz4/KSV/fvn1x/vx59OzZkxdyEBERERUSKid9e/bswYEDB9CgQYP8iIeIiIiI8oHK5/QVL14cFhYW+RELEREREeUTlZO+BQsWYMyYMYiKisqHcIiIiIgoP6g8vduzZ0+8fv0apUuXhomJSZYLOZ49e6ax4IiIiIhIM1RO+hYtWpQPYRARERFRflI56evTp09+xEFERERE+UjlpA8A0tPTsXPnTty4cQOSJKF8+fLw8fGBrq6upuMjIiIiIg1QOem7c+cOWrdujUePHsHNzQ1CCNy6dQvFixfHnj17ULp06fyIk4iIiIg+gspX7w4bNgylS5dGdHQ0Lly4gIsXL+LBgwcoWbIkhg0blh8xEhEREdFHUvlI37Fjx3DmzBlYWVnJZdbW1pg7dy7q16+v0eCIiIiISDNUTvoMDQ2RlJSUpfzly5cwMDDQSFBERERE2qZjx45K1w0JCVG5fZWnd9u2bYuBAwfi77//hhACQgicOXMGX3/9NXx8fFQOgIiIiIgAS0tLebGwsMCRI0dw7tw5ef358+dx5MgRWFpaqtW+ykf6lixZgj59+qBu3bryjZnT0tLg4+ODxYsXqxUEERERkbYLCgqS/x47diy6du2KFStWyHdHSU9Px+DBg9V+HK4khBDqbHj79m3cvHkTQgiUL18eZcqUUSuAz11iYiIsLS2RkJDAZw4TEVGuPsVvRmYfz2+VgoW5+rdKS0xKR9Fy9/j79pmytbXFyZMn4ebmplAeERGBevXqIT4+XuU21bpPHwCULVsWZcuWVXdzIiIiIspBWloabty4kSXpu3HjBjIyMtRqU+WkLz09HcHBwThy5AiePHmSpeOjR4+qFQgRERERvdOvXz/4+fnhzp07qFOnDgDgzJkzmDt3Lvr166dWmyonfcOHD0dwcDDatGmDihUrQpIktTomIiIiouzNnz8fDg4O+OmnnxATEwMAcHR0xJgxY/Ddd9+p1abK5/TZ2Nhg7dq1aN26tVodFjY8p4+IiJTFc/ooPyQmJgLAR+8nlY/0GRgY/Gcv2iAiIiL63GgqKVc56fvuu++wePFi/Pzzz5zaJSIiItKQatWq4ciRIyhatCiqVq2aa5514cIFldtXOek7efIkQkNDsW/fPlSoUEG+V18mde4QTURERKTtfH19YWhoCABo3769xttXOekrUqQIOnTooPFAiIiIiLTZ1KlTs/07N5s2bYKPjw9MTU3zrKty0vf+3aKJiIiIqOAMGjQItWvXRqlSpfKsq/Kzd7Pz/PlzLF26FJ6enppojoiIiIiUoMpNWNR+IgcAHD58GIGBgdi5cydsbGzQsWPHj2mOiIiIiPKJyknfgwcPEBQUhKCgILx8+RLPnz/H77//jk6dOuVHfERERESkAUpP7/7+++9o0aIFPDw8cPXqVSxevBiPHz+Gjo4OPDw88jNGIiIiIvpISh/p69GjB8aMGYPt27fD3Nw8P2MiIiIiIg1T+kifn58fli9fjpYtW2LFihV4/vx5fsZFRERERHlwcXHJcs/knCid9K1atQoxMTEYOHAgNm3aBEdHR/j6+kIIgYyMDLWDVUZYWBgkScp2CQ8Pz3XbGzduwMfHB5aWljA3N0edOnXw4MGDfI2XiIiI6FO4evUqihcvrlRdlW7ZYmxsjD59+uDYsWO4cuUKypcvD3t7e9SvXx89evTIt6dx1KtXDzExMQpL//794erqiho1auS43d27d9GgQQO4u7sjLCwMly9fxuTJk2FkZJQvcRIRERGpq2jRorCyslJqUYckVLnBSzYyMjKwZ88eBAYGYt++fUhJSfmY5pSSmpqKYsWKYejQoZg8eXKO9bp16wZ9fX2sW7dO7b4SExNhaWmJhIQEjT3wmIiI/ps+xW9GZh/Pb5WChbmu+u0kpaNouXv8ffuMrFmzRum6ffr0Ubn9j0763vfkyRPY2dlpqrkcbd++HV27dkVUVFSOhzQzMjJgaWmJMWPG4OTJk7h48SJKliyJ8ePH5/o8u5SUFIXENTExEcWLF+f/FERElCcmffQ508gTOTJ9ioQPAAIDA+Ht7Z3rHPaTJ0/w8uVLzJ07Fy1btsTBgwfRoUMHdOzYEceOHctxuzlz5sDS0lJelJ0nJyIiItKku3fvYtKkSejevTuePHkCANi/fz+uXbumVnsaTfpUFRAQkOMFGpnLuXPnFLZ5+PAhDhw4AH9//1zbzry4xNfXFyNHjoSnpyfGjRuHtm3bYsWKFTluN378eCQkJMhLdHT0xw+UiIiISAXHjh1DpUqV8PfffyMkJAQvX74EAPzzzz+YOnWqWm1+1GPYPtbQoUPRrVu3XOu4uroqvA4KCoK1tTV8fHxy3c7GxgZ6enooX768QrmHhwdOnjyZ43aGhoYwNDTMPXAiIiKifDRu3DjMnDkTo0aNUrg/cuPGjbF48WK12izQpM/GxgY2NjZK1xdCICgoCL17987znjQGBgaoWbMmIiIiFMpv3boFFxcXteIlIiIi+hSuXLmCjRs3Zim3tbVFfHy8Wm2qnfS9ffsWT548yXKPvhIlSqjbZJ6OHj2KyMjIHKd23d3dMWfOHHTo0AEA8P333+PLL7/EF198gcaNG2P//v34448/EBYWlm8xEhERfQodylWCnqTcTXmzkyZSAdxDzZo1oauriyFDhmDIkCGaC5A+SpEiRRATE4OSJUsqlF+8eBHOzs5qtaly0nf79m34+fnh1KlTCuVCCEiShPT0dLUCUUZgYCDq1auX47N+IyIikJCQIL/u0KEDVqxYgTlz5mDYsGFwc3PD9u3b0aBBg3yLkYiIqDAJDw/n1bufoR49emDs2LHYunUrJElCRkYG/vrrL4wePRq9e/dWq02Vb9lSv3596OnpYdy4cXB0dIQkSQrrq1SpolYgnyvep4+IiJT1KW/Z4gXfjz7SF4Zd/H37TKWmpqJv377YvHkzhBDQ09NDeno6evTogeDgYOjqqn67HpWP9F26dAnnz5+Hu7u7yp0RERERUd709fWxYcMGTJ8+HRcvXkRGRgaqVq2KsmXLqt2myklf+fLlERcXp3aHRERERKSc0qVLo3Tp0hppS+Wk74cffsCYMWMwe/ZsVKpUKctVtDxETERERKS6UaNGKV134cKFKrevctLXrFkzAEDTpk0Vyj/FhRxERERE/1UXL15UeH3+/Hmkp6fDzc0NwLvbzunq6qJ69epqta9y0hcaGqpWR0RERESUs/dzrIULF8Lc3Bxr1qxB0aJFAQDPnz9Hv3790LBhQ7XaV/nqXW3Dq3eJiEhZvHqXNMXZ2RkHDx5EhQoVFMqvXr2KFi1a4PHjxyq3qdbNmV+8eIHAwEDcuHEDkiShfPny8PPzg6WlpTrNEREREdF7EhMT8e+//2ZJ+p48eYKkpCS12tRRdYNz586hdOnS+Omnn/Ds2TPExcVh4cKFKF26NC5cuKBWEERERET0Px06dEC/fv2wbds2PHz4EA8fPsS2bdvg7++Pjh07qtWmykf6Ro4cCR8fH/z666/Q03u3eVpaGvr3748RI0bg+PHjagVCRERERO+sWLECo0ePRs+ePZGamgoA0NPTg7+/P3788Ue12lT5nD5jY2NcvHgxy82Zr1+/jho1auD169dqBfK54jl9RESkLJ7TR5r26tUr3L17F0IIlClTBqampmq3pfL0roWFBR48eJClPDo6Gubm5moHQkRERESKTE1NYWVlBRsbm49K+AA1kr4vv/wS/v7+2LJlC6Kjo/Hw4UNs3rwZ/fv3R/fu3T8qGCIiIiICMjIyMH36dFhaWsLFxQUlSpRAkSJFMGPGDGRkZKjVpsrn9M2fPx+SJKF3795IS0sD8O75cN988w3mzp2rVhBERERE9D8TJ05EYGAg5s6di/r160MIgb/++gsBAQFITk7GrFmzVG5T7fv0vX79WmGO2cTERJ1mPns8p4+IiJTFc/pIU5ycnLBixQr4+PgolO/atQuDBw/Go0ePVG5Trfv0AYCJiQkqVaqk7uZERERElINnz55luWgWANzd3fHs2TO12lQq6evYsSOCg4NhYWGR571hQkJC1AqEiIiIiN6pUqUKfv75ZyxZskSh/Oeff0aVKlXUalOppM/S0hKSJAF4d/Vu5t9EREREpHnz5s1DmzZtcPjwYdStWxeSJOHUqVN48OAB9u3bp1abfPZuHnhOHxERKYvn9JEmPXr0CL/88gtu3LgBIQTKly+PwYMHw8nJSa32VD6nr0mTJggJCUGRIkUUyhMTE9G+fXscPXpUrUCIiIiI6H+sra3h4+ODOnXqyLdpOXfuHABkucBDGSonfWFhYXj79m2W8uTkZJw4cULlAIiIiIhI0f79+9G7d2/Ex8fjw0lZSZKQnp6ucptKJ33//POP/Pf169cRGxsrv05PT8f+/fvh7OyscgBEREREpGjo0KHo0qULpkyZAnt7e420qXTS5+npCUmSIEkSmjRpkmW9sbExli5dqpGgiIiIiLTZkydPMGrUKI0lfIAKSV9kZCSEEChVqhTOnj0LW1tbeZ2BgQHs7Oygq6urscCIiIiItFXnzp0RFhaG0qVLa6xNpZM+FxcXAFD7eW9EREREpJyff/4ZXbp0wYkTJ1CpUiXo6yteqT1s2DCV21T7iRzXr1/HgwcPslzUoc7VJERERET0Pxs3bsSBAwdgbGyMsLAwhXskS5L0aZK+e/fuoUOHDrhy5QokSZKvKMkMRp2rSYiIiIjofyZNmoTp06dj3Lhx0NHR0UibKrcyfPhwlCxZEv/++y9MTExw7do1HD9+HDVq1EBYWJhGgiIiIiLSZm/fvsWXX36psYQPUCPpO336NKZPnw5bW1vo6OhAR0cHDRo0wJw5c9Q61EhEREREivr06YMtW7ZotE2Vp3fT09NhZmYGALCxscHjx4/h5uYGFxcXREREaDQ4IiIiIm2Unp6OefPm4cCBA6hcuXKWCzkWLlyocpsqJ30VK1bEP//8g1KlSqF27dqYN28eDAwMsGrVKpQqVUrlAIiIiIhI0ZUrV1C1alUAwNWrVxXWvX9RhypUTvomTZqEV69eAQBmzpyJtm3bomHDhrC2tsbmzZvVCoKIiIiI/ic0NFTjbaqc9Hl7e8t/lypVCtevX8ezZ89QtGhRtTNPIiIiIspfKl/I4efnh6SkJIUyKysrvH79Gn5+fhoLjIiIiIg0R+Wkb82aNXjz5k2W8jdv3mDt2rUaCYqIiIiINEvp6d3ExEQIISCEQFJSEoyMjOR16enp2Lt3L+zs7PIlSCIiIiL6OEonfUWKFIEkSZAkCeXKlcuyXpIkTJs2TaPBEREREZFmKJ30hYaGQgiBJk2aYPv27bCyspLXGRgYwMXFBU5OTvkSJBERERF9HKWTvkaNGgEAIiMjUaJECV6pS0RERFSIKJX0/fPPP6hYsSJ0dHSQkJCAK1eu5Fi3cuXKGguOiIiIiDRDqaTP09MTsbGxsLOzg6enJyRJghAiSz1JkpCenq7xIImIiIjo4yiV9EVGRsLW1lb+m4iIiIgKF6WSPhcXl2z/JiIiIqLCQeXHsAFAREQEli5dihs3bkCSJLi7u+Pbb7+Fm5ubpuMjIiIiIg1Q+Ykc27ZtQ8WKFXH+/HlUqVIFlStXxoULF1CxYkVs3bo1P2IkIiIioo+k8pG+MWPGYPz48Zg+fbpC+dSpUzF27Fh06dJFY8ERERERkWaofKQvNjYWvXv3zlLes2dPxMbGaiQoIiIiItIslZM+Ly8vnDhxIkv5yZMn0bBhQ40ERURERESapfL0ro+PD8aOHYvz58+jTp06AIAzZ85g69atmDZtGnbv3q1Ql4iIiIgKniSyu8tyLnR0lDs4+F+5UXNiYiIsLS2RkJAACwuLgg6HiIg+Y5/iNyOzDy/4Qk/SV7udNJGKMOzi75sWUflIX0ZGRn7EQURERET5SOVz+oiIiIio8FHqSN+SJUswcOBAGBkZYcmSJbnWHTZsmEYCe19YWBgaN26c7bqzZ8+iZs2a2a57+fIlxo0bh507dyI+Ph6urq4YNmwYvvnmG43HSERERPQ5U+qcvpIlS+LcuXOwtrZGyZIlc25MknDv3j2NBggAb9++xbNnzxTKJk+ejMOHD+PevXuQJCnb7QYMGIDQ0FD89ttvcHV1xcGDBzF48GBs374dvr6+SvXNc/qIiEhZPKePPmdKHemLjIzM9u9PxcDAAA4ODvLr1NRU7N69G0OHDs0x4QOA06dPo0+fPvDy8gIADBw4ECtXrsS5c+dyTPpSUlKQkpIiv05MTNTMIIiIiIgKUKE8p2/37t2Ii4tD3759c63XoEED7N69G48ePYIQAqGhobh16xa8vb1z3GbOnDmwtLSUl+LFi2s4eiIiIqJPT+Wkr3Pnzpg7d26W8h9//PGTPYItMDAQ3t7eeSZkS5YsQfny5VGsWDEYGBigZcuWWL58ORo0aJDjNuPHj0dCQoK8REdHazp8IiIiok9O5aTv2LFjaNOmTZbyli1b4vjx4yq1FRAQAEmScl3OnTunsM3Dhw9x4MAB+Pv759n+kiVLcObMGezevRvnz5/HggULMHjwYBw+fDjHbQwNDWFhYaGwEBERERV2Kt+n7+XLlzAwMMhSrq+vr/L5b0OHDkW3bt1yrePq6qrwOigoCNbW1nk+7ePNmzeYMGECduzYISeplStXxqVLlzB//nw0a9ZMpViJiIiICjOVk76KFStiy5YtmDJlikL55s2bUb58eZXasrGxgY2NjdL1hRAICgpC7969oa+f+xVLqampSE1NzfIEEV1dXd5gmoiIiLSOyknf5MmT0alTJ9y9exdNmjQBABw5cgSbNm3C1q1bNR7g+44ePYrIyMgcp3bd3d0xZ84cdOjQARYWFmjUqBG+//57GBsbw8XFBceOHcPatWuxcOHCfI2TiIiI6HOjctLn4+ODnTt3Yvbs2di2bRuMjY1RuXJlHD58GI0aNcqPGGWBgYGoV68ePDw8sl0fERGBhIQE+fXmzZsxfvx4fPXVV3j27BlcXFwwa9YsfP311/kaJxEREdHnRqmbM2sz3pyZiIiUxZsz0+dMrfv0vXjxAr/99hsmTJggPynjwoULePTokUaDIyIiIiLNUHl6959//kGzZs1gaWmJqKgo9O/fH1ZWVtixYwfu37+PtWvX5kecRERERPQRVD7SN2rUKPTt2xe3b9+GkZGRXN6qVSuV79NHRERERJ+GyklfeHg4Bg0alKXc2dkZsbGxGgmKiIiIiDRL5aTPyMgo25swR0REwNbWViNBEREREZFmqZz0+fr6Yvr06UhNTQUASJKEBw8eYNy4cejUqZPGAyQiIiKij6dy0jd//nw8ffoUdnZ2ePPmDRo1aoQyZcrA3Nwcs2bNyo8YiYiIiOgjqXz1roWFBU6ePImjR4/iwoULyMjIQLVq1fgsWyIiIqLPmEpJX1paGoyMjHDp0iU0adJEfgwbEREREX3eVJre1dPTg4uLC9LT0/MrHiIiIiLKByqf0zdp0iSMHz9efhIHEREREX3+VD6nb8mSJbhz5w6cnJzg4uICU1NThfUXLlzQWHBEREREpBkqJ32+vr6QJCk/YiEiIiKifKJy0hcQEJAPYRARERFRflL6nL7Xr19jyJAhcHZ2hp2dHXr06IG4uLj8jI2IiIiINETppG/q1KkIDg5GmzZt0K1bNxw6dAjffPNNfsZGRERERBqi9PRuSEgIAgMD0a1bNwBAz549Ub9+faSnp0NXVzffAiQiIiKij6f0kb7o6Gg0bNhQfl2rVi3o6enh8ePH+RIYEREREWmO0klfeno6DAwMFMr09PSQlpam8aCIiIiISLOUnt4VQqBv374wNDSUy5KTk/H1118r3KsvJCREsxESERER0UdTOunr06dPlrKePXtqNBgiIiIiyh9KJ31BQUH5GQcRERER5SOVn71LRERERIUPkz4iIiIiLcCkj4iIiJQWGxuLb7/9FqVKlYKhoSGKFy+Odu3a4ciRIwAAV1dXSJKEzZs3Z9m2QoUKkCQJwcHBCuUXL15Ely5dYG9vDyMjI5QrVw4DBgzArVu3FOpt374dXl5esLS0hJmZGSpXrozp06fj2bNnOcYbEBAASZKyLO7u7nKde/fuoXv37nBycoKRkRGKFSsGX19fhf4ztztz5oxC+ykpKbC2toYkSQgLCwMAREVFwd/fHyVLloSxsTFKly6NqVOn4u3bt0q9x/mFSR8REREpJSoqCtWrV8fRo0cxb948XLlyBfv370fjxo0xZMgQuV7x4sWzXAtw5swZxMbGKtzxAwD+/PNP1KlTBykpKdiwYQNu3LiBdevWwdLSEpMnT5brTZw4EV9++SVq1qyJffv24erVq1iwYAEuX76MdevW5Rp3hQoVEBMTo7CcPHkSAPD27Vs0b94ciYmJCAkJQUREBLZs2YKKFSsiISFBoZ3sxrVjxw6YmZkplN28eRMZGRlYuXIlrl27hp9++gkrVqzAhAkT8niH85ckhBAFGsFnLjExEZaWlkhISICFhUVBh0NERJ+xT/GbkdmHF3yhJ+mr3U6aSEUYdiE6OlohVkNDQ4Xbs72vdevW+OeffxAREZEleXvx4gWKFCkCV1dXdO/eHT/99BNu376N4sWLAwAGDhwIIyMjrF27FosWLULfvn3x+vVruLi4oEGDBtixY0eW/jLbPHv2LGrXro1FixZh+PDhOdbLTkBAAHbu3IlLly5lu/7SpUuoWrUqoqKi4OLikm0d4N2RvkmTJmHJkiWIjY2FsbExAKBFixaoU6cOZsyYgdDQUHh5eWW7/Y8//ohffvkF9+7dy7GP/MYjfURERFqsePHisLS0lJc5c+ZkW+/Zs2fYv38/hgwZkiXhA6CQdNnb28Pb2xtr1qwBALx+/RpbtmyBn5+fwjYHDhxAXFwcxowZk22fmW1u2LABZmZmGDx4cK711GFrawsdHR1s27YN6enpudatXr06SpYsie3btwN497Sy48ePo1evXnn2k5CQACsrK7Xj1AQmfURERFosOjoaCQkJ8jJ+/Phs6925cwdCCIVz4XLj5+eH4OBgCCGwbds2lC5dGp6engp1bt++DQB5tnn79m2UKlUK+vrqHdm8cuUKzMzMFJb+/fsDAJydnbFkyRJMmTIFRYsWRZMmTTBjxowcj8j169cPq1evBvDudnatW7eGra1trv3fvXsXS5cuxddff61W/JrCpI+IiEiLWVhYKCw5Te1mng0mSZJS7bZp0wYvX77E8ePHsXr16ixH+d5vMy9CiDz7ffDggUJSN3v2bHmdm5sbLl26pLDMmjVLXj9kyBDExsZi/fr1qFu3LrZu3YoKFSrg0KFDWfrp2bMnTp8+jXv37iE4ODjbcb3v8ePHaNmyJbp06SInmgVF6ZszExERkfYqW7YsJEnCjRs30L59+zzr6+npoVevXpg6dSr+/vvvbM/ZK1euHIB3Fz7UrVs3x7bKlSuHkydPIjU1NcejfU5OTgrn7b0/lWpgYIAyZcrkGq+5uTl8fHzg4+ODmTNnwtvbGzNnzkTz5s0V6llbW6Nt27bw9/dHcnIyWrVqhaSkpGzbfPz4MRo3boy6deti1apVufb/KfBIHxEREeXJysoK3t7eWLZsGV69epVl/YsXL7KU+fn54dixY/D19UXRokWzrG/RogVsbGwwb968bPvMbLNHjx54+fIlli9fnmM9PT09lClTRl4+5vy5zFu6ZDdO4N24wsLC0Lt3b+jq6mZb59GjR/Dy8kK1atUQFBQEHZ2CT7l4pI+IiIiUsnz5ctSrVw+1atXC9OnTUblyZaSlpeHQoUP45ZdfcOPGDYX6Hh4eiIuLg4mJSbbtmZqa4rfffkOXLl3g4+ODYcOGoUyZMoiLi8Pvv/+OBw8eYPPmzahduzbGjBmD7777Do8ePUKHDh3g5OSEO3fuYMWKFWjQoEG2V/VmSktLQ2xsrEKZJEmwt7fHpUuXMHXqVPTq1Qvly5eHgYEBjh07htWrV2Ps2LHZtteyZUs8ffo0xyu0Hz9+DC8vL5QoUQLz58/H06dP5XUODg45xpnfmPQRERGRUkqWLIkLFy5g1qxZ+O677xATEwNbW1tUr14dv/zyS7bbWFtb59qmr68vTp06hTlz5qBHjx5ITExE8eLF0aRJE8ycOVOu98MPP6B69epYtmwZVqxYgYyMDJQuXRqdO3dGnz59cu3j2rVrcHR0VCgzNDREcnIyihUrBldXV0ybNg1RUVGQJEl+PXLkyGzbkyQJNjY2OfZ38OBB3LlzB3fu3EGxYsUU1hXknfJ4n7488D59RESkrMJ4nz7+vmmPgp9gJiIiIqJ8x6SPiIiISAsw6SMiIiLSAkz6iIiIiLQAkz4iIiIiLcCkj4iIiEgLMOkjIiIi0gJM+oiIiIi0AJM+IiIiIi3ApI+IiIhICzDpIyIiItICTPqIiIiItACTPiIiIiItwKSPiIiISAsw6SMiIiLSAoUm6bt16xZ8fX1hY2MDCwsL1K9fH6GhobluI4RAQEAAnJycYGxsDC8vL1y7du0TRUxERET0+Sg0SV+bNm2QlpaGo0eP4vz58/D09ETbtm0RGxub4zbz5s3DwoUL8fPPPyM8PBwODg5o3rw5kpKSPmHkRERERAWvUCR9cXFxuHPnDsaNG4fKlSujbNmymDt3Ll6/fp3jkTshBBYtWoSJEyeiY8eOqFixItasWYPXr19j48aNn3gERERERAWrUCR91tbW8PDwwNq1a/Hq1SukpaVh5cqVsLe3R/Xq1bPdJjIyErGxsWjRooVcZmhoiEaNGuHUqVM59pWSkoLExESFhYiIiKiw0yvoAJQhSRIOHToEX19fmJubQ0dHB/b29ti/fz+KFCmS7TaZ07729vYK5fb29rh//36Ofc2ZMwfTpk3TWOxEREREn4MCPdIXEBAASZJyXc6dOwchBAYPHgw7OzucOHECZ8+eha+vL9q2bYuYmJhc+5AkSeG1ECJL2fvGjx+PhIQEeYmOjtbIWImIiIgKUoEe6Rs6dCi6deuWax1XV1ccPXoUf/75J54/fw4LCwsAwPLly3Ho0CGsWbMG48aNy7Kdg4MDgHdH/BwdHeXyJ0+eZDn69z5DQ0MYGhqqMxwiIiKiz1aBJn02NjawsbHJs97r168BADo6igcmdXR0kJGRke02JUuWhIODAw4dOoSqVasCAN6+fYtjx47hhx9++MjIiYiIiAqXQnEhR926dVG0aFH06dMHly9fxq1bt/D9998jMjISbdq0keu5u7tjx44dAN5N644YMQKzZ8/Gjh07cPXqVfTt2xcmJibo0aNHQQ2FiIiIqEAUigs5bGxssH//fkycOBFNmjRBamoqKlSogF27dqFKlSpyvYiICCQkJMivx4wZgzdv3mDw4MF4/vw5ateujYMHD8Lc3LwghkFERERUYCQhhCjoID5niYmJsLS0REJCgnw+IRERUXY+xW9GZh9e8IWepK92O2kiFWHYxd83LVIopneJiIiI6OMw6SMiIiLSAkz6iIiIiLQAkz4iIiIiLcCkj4iIiEgLMOkjIiIi0gJM+oiIiIi0AJM+IiIiIi3ApI+IiIhICzDpIyIiItICTPqIiIiItACTPiIiIiItwKSPiIiISAsw6SMiIiLSAkz6iIiIiLQAkz4iIiIiLcCkj4iIiEgLMOkjIiIi0gJM+oiIiIi0AJM+IiIiIi3ApI+IiIhICzDpIyIiItICTPqIiIiItACTPiIiIiItwKSPiIiISAsw6SMiIiLSAkz6iIiIiLQAkz4iIiIiLcCkj4iIiEgLMOkjIiIi0gJM+oiIiIi0AJM+IiIiIi3ApI+IiIhICzDpIyIiItICTPqIiIiItACTPiIiIiItwKSPiIiISAvoFXQA2kAIASDt/1/pQpKYaxMREdGnxaQvnwkhgIwY4M0uCJEKyagVoF+2oMMiIiIiLcOkLx8JkQG82QaROAVAxruyVz9DmH4DHfORBRscERERaRXOM+Yn8RIicToyEz7Zq18g0u6+SwqJiIiIPgEmfflEiDQgJRTA2+wrJO8DkP4pQyIiIiItxqQv30gA9HNZb/CpAiEiIiJi0pdfJEkXMGoKSBbZrNUFjH3BUyqJiIjoU2HSl6/0IBVZAkiW75UZQbKcA+jYQpKkAouMiIiItAsPNeUjSdKFMKgNye4kkBIGiFTAsDEgGfNefURERPRJMenLZ5KkC0AXwrDZe6+JiIiIPi0mfZ8Ikz0iIiIqSJxjJCIiItICTPqIiIiItACTPiIiIiItUGiSvlu3bsHX1xc2NjawsLBA/fr1ERoammP91NRUjB07FpUqVYKpqSmcnJzQu3dvPH78+BNGTURERPR5KDRJX5s2bZCWloajR4/i/Pnz8PT0RNu2bREbG5tt/devX+PChQuYPHkyLly4gJCQENy6dQs+Pj6fOHIiIiKigicJIURBB5GXuLg42Nra4vjx42jYsCEAICkpCRYWFjh8+DCaNm2qVDvh4eGoVasW7t+/jxIlSmRbJyUlBSkpKfLrxMREFC9eHAkJCbCwyO7pGkRERO8kJibC0tIyX38zMvvwgi/0pNwe95m7NJGKMOzi75sWKRRH+qytreHh4YG1a9fi1atXSEtLw8qVK2Fvb4/q1asr3U5CQgIkSUKRIkVyrDNnzhxYWlrKS/HixTUwAiIiIqKCVSiSPkmScOjQIVy8eBHm5uYwMjLCTz/9hP379+eawL0vOTkZ48aNQ48ePXL9F8348eORkJAgL9HR0RoaBREREVHBKdCkLyAgAJIk5bqcO3cOQggMHjwYdnZ2OHHiBM6ePQtfX1+0bdsWMTExefaTmpqKbt26ISMjA8uXL8+1rqGhISwsLBQWIiIiosKuQM/pi4uLQ1xcXK51XF1d8ddff6FFixZ4/vy5QhJWtmxZ+Pv7Y9y4cTlun5qaiq5du+LevXs4evQorK2tVYrxU5yfQURE/w08p48+ZwX6GDYbGxvY2NjkWe/169cAAB0dxQOTOjo6yMjIyHG7zITv9u3bCA0NVTnhA4DMnDgxMVHlbYmISLtk/lZ8iuMpaUgFPqKbNKRqLhgqHEQh8PTpU2FtbS06duwoLl26JCIiIsTo0aOFvr6+uHTpklzPzc1NhISECCGESE1NFT4+PqJYsWLi0qVLIiYmRl5SUlKU7js6Olrg3f9WXLhw4cKFi1JLdHS0xn8LM71580Y4ODhoJE4LCwvh5uYmPDw8xM8//5xvMdPnoUCP9CnLxsYG+/fvx8SJE9GkSROkpqaiQoUK2LVrF6pUqSLXi4iIQEJCAgDg4cOH2L17NwDA09NTob3Q0FB4eXkp1beTkxOio6Nhbm4OSZI0Mp6Clnkbmujo6P/cIf3/6tg4rsLnvzo2jit3QggkJSXByclJg9EpMjIyQmRkJN6+ffvRbRkYGMDIyEgDUVFhUCju00ea9V8+T/G/OjaOq/D5r46N4yIqvArFLVuIiIiI6OMw6SMiIiLSAkz6tJChoSGmTp0KQ0PDgg5F4/6rY+O4Cp//6tg4LqLCi+f0EREREWkBHukjIiIi0gJM+oiIiIi0AJM+IiIiIi3ApI+IiIhICzDp+4/55ZdfULlyZVhYWMDCwgJ169bFvn37lNr2r7/+gp6eXpYnmHwu1BlbSkoKJk6cCBcXFxgaGqJ06dJYvXr1J4pYOeqMa8OGDahSpQpMTEzg6OiIfv36IT4+/hNFrJ45c+ZAkiSMGDEi13rHjh1D9erVYWRkhFKlSmHFihWfJkA1KTOukJAQNG/eHLa2tvI+PnDgwKcLUk3K7rNMn/t3SCZlx1UYvj+IVMGk7z+mWLFimDt3Ls6dO4dz586hSZMm8PX1xbVr13LdLiEhAb1790bTpk0/UaSqU2dsXbt2xZEjRxAYGIiIiAhs2rQJ7u7unzDqvKk6rpMnT6J3797w9/fHtWvXsHXrVoSHh6N///6fOHLlhYeHY9WqVahcuXKu9SIjI9G6dWs0bNgQFy9exIQJEzBs2DBs3779E0WqGmXHdfz4cTRv3hx79+7F+fPn0bhxY7Rr1w4XL178RJGqTtmxZSoM3yGAauMqDN8fRCop0Cf/0idRtGhR8dtvv+Va58svvxSTJk0SU6dOFVWqVPk0gWlAbmPbt2+fsLS0FPHx8Z84qo+X27h+/PFHUapUKYWyJUuWiGLFin2K0FSWlJQkypYtKw4dOiQaNWokhg8fnmPdMWPGCHd3d4WyQYMGiTp16uRzlKpTZVzZKV++vJg2bVr+BPeR1BlbYfgOUWVchfn7gygnPNL3H5aeno7Nmzfj1atXqFu3bo71goKCcPfuXUydOvUTRvdxlBnb7t27UaNGDcybNw/Ozs4oV64cRo8ejTdv3nziaJWnzLjq1auHhw8fYu/evRBC4N9//8W2bdvQpk2bTxytcoYMGYI2bdqgWbNmedY9ffo0WrRooVDm7e2Nc+fOITU1Nb9CVIsq4/pQRkYGkpKSYGVllQ+RfTxVx1ZYvkNUGVdh/P4gyoteQQdAmnflyhXUrVsXycnJMDMzw44dO1C+fPls696+fRvjxo3DiRMnoKf3+X8cVBnbvXv3cPLkSRgZGWHHjh2Ii4vD4MGD8ezZs8/uvBxVxlWvXj1s2LABX375JZKTk5GWlgYfHx8sXbr0E0edt82bN+PChQsIDw9Xqn5sbCzs7e0Vyuzt7ZGWloa4uDg4OjrmR5gqU3VcH1qwYAFevXqFrl27ajiyj6fq2ArLd4iq4ypM3x9EyuKRvv8gNzc3XLp0CWfOnME333yDPn364Pr161nqpaeno0ePHpg2bRrKlStXAJGqTtmxAe+OpkiShA0bNqBWrVpo3bo1Fi5ciODg4M/uX+uqjOv69esYNmwYpkyZgvPnz2P//v2IjIzE119//Ymjzl10dDSGDx+O9evXw8jISOntJElSeC3+/6FBH5YXFHXHlWnTpk0ICAjAli1bYGdnlw8Rqk/VsRWW7xB19llh+v4gUlpBzy9T/mvatKkYOHBglvLnz58LAEJXV1deJEmSy44cOVIA0aomp7EJIUTv3r1F6dKlFcquX78uAIhbt259ivDUltu4evbsKTp37qxQduLECQFAPH78+FOEp5QdO3Zk+XwBEJIkCV1dXZGWlpZlm4YNG4phw4YplIWEhAg9PT3x9u3bTxV6rtQZV6bNmzcLY2Nj8eeff37CiJWn6tgKy3eIOvusMH9/EOXk8z0WTxojhEBKSkqWcgsLC1y5ckWhbPny5Th69Ci2bduGkiVLfqoQ1ZbT2ACgfv362Lp1K16+fAkzMzMAwK1bt6Cjo4NixYp9yjBVltu4Xr9+nWUaTVdXV97uc9G0adMsn69+/frB3d0dY8eOlWN+X926dfHHH38olB08eBA1atSAvr5+vsarLHXGBbw7wufn54dNmzZ9tudfqjq2wvIdos4+K8zfH0Q5KtickzRt/Pjx4vjx4yIyMlL8888/YsKECUJHR0ccPHhQCCHEuHHjRK9evXLc/nO+8k7VsSUlJYlixYqJzp07i2vXroljx46JsmXLiv79+xfUELKl6riCgoKEnp6eWL58ubh79644efKkqFGjhqhVq1ZBDUFpH14x+eHY7t27J0xMTMTIkSPF9evXRWBgoNDX1xfbtm0rgGiVl9e4Nm7cKPT09MSyZctETEyMvLx48aIAolVNXmP70Of8HfK+vMZVWL4/iFTBI33/Mf/++y969eqFmJgYWFpaonLlyti/fz+aN28OAIiJicGDBw8KOEr1qDo2MzMzHDp0CN9++y1q1KgBa2trdO3aFTNnziyoIWRL1XH17dsXSUlJ+Pnnn/Hdd9+hSJEiaNKkCX744YeCGoLaPhxbyZIlsXfvXowcORLLli2Dk5MTlixZgk6dOhVglKr7cFwrV65EWloahgwZgiFDhsjlffr0QXBwcAFEqL7C/B2Sm8L6/UGkCkmIz2g+iIiIiIjyBa/eJSIiItICTPqIiIiItACTPiIiIiItwKSPiIiISAsw6SMiIiLSAkz6iIiIiLQAkz4iIiIiLcCkj4iIiEgLMOkjIspDWFgYJEnCixcvCjoUIiK1MekjrSJJUq5L3759CzpEjfPy8sKIESMKOgwAwMCBA6Grq4vNmzcXdCifXN++fdG+fXuFsm3btsHIyAjz5s0rmKCISKvw2bukVWJiYuS/t2zZgilTpiAiIkIuMzY2Loiw1JKamgp9ff1C09/r16+xZcsWfP/99wgMDES3bt00GF3h89tvv2HIkCFYtmwZ+vfvX9DhEJEW4JE+0ioODg7yYmlpCUmSFMqOHz+O6tWrw8jICKVKlcK0adOQlpYmby9JElauXIm2bdvCxMQEHh4eOH36NO7cuQMvLy+Ympqibt26uHv3rrxNQEAAPD09sXLlShQvXhwmJibo0qVLlqnCoKAgeHh4wMjICO7u7li+fLm8LioqCpIk4ffff4eXlxeMjIywfv16xMfHo3v37ihWrBhMTExQqVIlbNq0Sd6ub9++OHbsGBYvXiwfzYyKikJwcDCKFCmi0P/OnTshSVKWuFevXo1SpUrB0NAQQggkJCRg4MCBsLOzg4WFBZo0aYLLly/n+d5v3boV5cuXx/jx4/HXX38hKipKYX3mkbD58+fD0dER1tbWGDJkCFJTU+U6z58/R+/evVG0aFGYmJigVatWuH37trw+c1x//vkn3NzcYGJigs6dO+PVq1dYs2YNXF1dUbRoUXz77bdIT0+Xt1u/fj1q1KgBc3NzODg4oEePHnjy5Em243j16hUsLCywbds2hfI//vgDpqamSEpKyvO9mDdvHoYOHYqNGzcy4SOiT0cQaamgoCBhaWkpv96/f7+wsLAQwcHB4u7du+LgwYPC1dVVBAQEyHUACGdnZ7FlyxYREREh2rdvL1xdXUWTJk3E/v37xfXr10WdOnVEy5Yt5W2mTp0qTE1NRZMmTcTFixfFsWPHRJkyZUSPHj3kOqtWrRKOjo5i+/bt4t69e2L79u3CyspKBAcHCyGEiIyMFACEq6urXOfRo0fi4cOH4scffxQXL14Ud+/eFUuWLBG6urrizJkzQgghXrx4IerWrSsGDBggYmJiRExMjEhLS8sydiGE2LFjh3j/KyEzbm9vb3HhwgVx+fJlkZGRIerXry/atWsnwsPDxa1bt8R3330nrK2tRXx8fK7vd8OGDcXPP/8shBCiU6dOYsqUKQrr+/TpIywsLMTXX38tbty4If744w9hYmIiVq1aJdfx8fERHh4e4vjx4+LSpUvC29tblClTRrx9+1bep/r6+qJ58+biwoUL4tixY8La2lq0aNFCdO3aVVy7dk388ccfwsDAQGzevFluNzAwUOzdu1fcvXtXnD59WtSpU0e0atVKXh8aGioAiOfPnwshhBgwYIBo3bq1QvwdOnQQvXv3znH8ffr0Eb6+vmLs2LHCzMxMHDp0KNf3i4hI05j0kdb6MPFp2LChmD17tkKddevWCUdHR/k1ADFp0iT59enTpwUAERgYKJdt2rRJGBkZya+nTp0qdHV1RXR0tFy2b98+oaOjI2JiYoQQQhQvXlxs3LhRoe8ZM2aIunXrCiH+l/QtWrQoz3G1bt1afPfdd/LrRo0aieHDh+c6diGyT/r09fXFkydP5LIjR44ICwsLkZycrLBt6dKlxcqVK3OM6datW0JfX188ffpU7qt48eIiPT1drtOnTx/h4uIi0tLS5LIuXbqIL7/8Um4DgPjrr7/k9XFxccLY2Fj8/vvv8rgAiDt37sh1Bg0aJExMTERSUpJc5u3tLQYNGpRjvGfPnhUA5G0+TPr+/vtvoaurKx49eiSEEOLp06dCX19fhIWF5dhmnz59hIGBgQAgjhw5kmM9IqL8wuldov93/vx5TJ8+HWZmZvIyYMAAxMTE4PXr13K9ypUry3/b29sDACpVqqRQlpycjMTERLmsRIkSKFasmPy6bt26yMjIQEREBJ4+fYro6Gj4+/sr9D1z5kyFaWIAqFGjhsLr9PR0zJo1C5UrV4a1tTXMzMxw8OBBPHjwQCPviYuLC2xtbeXX58+fx8uXL+W+MpfIyMgssb4vMDAQ3t7esLGxAQC0bt0ar169wuHDhxXqVahQAbq6uvJrR0dHeZr1xo0b0NPTQ+3ateX11tbWcHNzw40bN+QyExMTlC5dWn5tb28PV1dXmJmZKZS9P3178eJF+Pr6wsXFBebm5vDy8gKAHN/HWrVqoUKFCli7di0AYN26dShRogS++OKLHN8D4N1nx9XVFVOmTFFqGpiISJN4IQfR/8vIyMC0adPQsWPHLOuMjIzkv9+/mCHzHLjsyjIyMnLsK7OOJElyvV9//VUhoQGgkAABgKmpqcLrBQsW4KeffsKiRYtQqVIlmJqaYsSIEXj79m3OAwWgo6MDIYRC2fvnzuXUX0ZGBhwdHREWFpal7ofnCGZKT0/H2rVrERsbCz09PYXywMBAtGjRQi778EKR99+fD+PNJIRQOBcxuzZya/fVq1do0aIFWrRogfXr18PW1hYPHjyAt7d3ru9j//798fPPP2PcuHEICgpCv379FOLIjrOzM7Zv347GjRujZcuW2L9/P8zNzXPdhohIU5j0Ef2/atWqISIiAmXKlNF42w8ePMDjx4/h5OQEADh9+jR0dHRQrlw52Nvbw9nZGffu3cNXX32lUrsnTpyAr68vevbsCeBdUnb79m14eHjIdQwMDBQuWgAAW1tbJCUl4dWrV3Jid+nSpTz7q1atmpy8ubq6KhXj3r17kZSUhIsXLyoksTdv3sRXX32F+Ph4WFtb59lO+fLlkZaWhr///hv16tUDAMTHx+PWrVsK41XVzZs3ERcXh7lz56J48eIAgHPnzuW5Xc+ePTFmzBgsWbIE165dQ58+fZTqr0SJEjh27BgaN26MFi1a4MCBA7CwsFA7fiIiZXF6l+j/TZkyBWvXrkVAQACuXbuGGzduYMuWLZg0adJHt21kZIQ+ffrg8uXLOHHiBIYNG4auXbvCwcEBwLsrZefMmYPFixfj1q1buHLlCoKCgrBw4cJc2y1TpgwOHTqEU6dO4caNGxg0aBBiY2MV6ri6uuLvv/9GVFQU4uLikJGRgdq1a8PExAQTJkzAnTt3sHHjRgQHB+c5jmbNmqFu3bpo3749Dhw4gKioKJw6dQqTJk3KMVEKDAxEmzZtUKVKFVSsWFFeOnXqBFtbW6xfv16p97Bs2bLw9fXFgAEDcPLkSVy+fBk9e/aEs7MzfH19lWojOyVKlICBgQGWLl2Ke/fuYffu3ZgxY0ae2xUtWhQdO3bE999/jxYtWihM3+elWLFiCAsLQ3x8PFq0aIGEhAS14yciUhaTPqL/5+3tjT///BOHDh1CzZo1UadOHSxcuBAuLi4f3XaZMmXQsWNHtG7dGi1atEDFihUVbsnSv39//PbbbwgODkalSpXQqFEjBAcHo2TJkrm2O3nyZFSrVg3e3t7w8vKCg4NDlhsAjx49Grq6uihfvrw8dWllZYX169dj79698m1eAgIC8hyHJEnYu3cvvvjiC/j5+aFcuXLo1q0boqKi5PMb3/fvv/9iz5496NSpU7ZtdezYEYGBgXn2mykoKAjVq1dH27ZtUbduXQghsHfv3o+6f6CtrS2Cg4PlW8rMnTsX8+fPV2pbf39/vH37Fn5+fir36+zsjGPHjuHFixdo3rw5n/ZBRPlOEjmdKENEGhEQEICdO3cqNX1KhcuGDRswfPhwPH78GAYGBgUdDhFRrnhOHxGRil6/fo3IyEjMmTMHgwYNYsJHRIUCp3eJiFQ0b948eHp6wt7eHuPHjy/ocIiIlMLpXSIiIiItwCN9RERERFqASR8RERGRFmDSR0RERKQFmPQRERERaQEmfURERERagEkfERERkRZg0kdERESkBZj0EREREWmB/wPi8DdccXsmYgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot deltaT-deltaP diagram\n",
    "\n",
    "# To use xarray scatter plot, create a new Dataset using both DataArrays\n",
    "ds_anomaly = xr.Dataset({\"TG\": full_tg_anomaly, \"PRCPTOT\": full_pr_anomaly})\n",
    "\n",
    "# Scatter Plot\n",
    "xr.plot.scatter(ds_anomaly, x=\"TG\", y=\"PRCPTOT\", hue=\"model_id\")\n",
    "plt.suptitle(\n",
    "    \"DeltaT-DeltaP Diagram  \"\n",
    "    + str(yearb)\n",
    "    + \"-\"\n",
    "    + str(yeare)\n",
    "    + \" vs \"\n",
    "    + str(yearrefb)\n",
    "    + \"-\"\n",
    "    + str(yearrefe)\n",
    ")\n",
    "plt.title(\n",
    "    \"Domain: lat=[\"\n",
    "    + str(minlat)\n",
    "    + \", \"\n",
    "    + str(maxlat)\n",
    "    + \"] lon=[\"\n",
    "    + str(minlon)\n",
    "    + \", \"\n",
    "    + str(maxlon)\n",
    "    + \"]\"\n",
    ")\n",
    "plt.ylabel(\"Precipitation Anomaly %\")\n",
    "plt.xlabel(\"Temperature Anomaly K\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}