code/NAM_CMAM_variance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "76f4115e-d671-4c56-8596-30b37bcddcff",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-03-04T15:58:14.752421Z",
     "iopub.status.busy": "2022-03-04T15:58:14.751725Z",
     "iopub.status.idle": "2022-03-04T15:58:14.758687Z",
     "shell.execute_reply": "2022-03-04T15:58:14.757848Z",
     "shell.execute_reply.started": "2022-03-04T15:58:14.752341Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from aostools.climate import GlobalAvg, eof\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "baccca0f-c69c-47a4-afdd-22f5ea5d7cfc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-03-04T15:58:17.943465Z",
     "iopub.status.busy": "2022-03-04T15:58:17.942970Z",
     "iopub.status.idle": "2022-03-04T15:58:17.964821Z",
     "shell.execute_reply": "2022-03-04T15:58:17.964032Z",
     "shell.execute_reply.started": "2022-03-04T15:58:17.943425Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def ComputeAnnularMode(lat, pres, data, choice='z', hemi='infer', detrend='constant', eof_in=None, pc_in=None, eof_out=False, pc_out=False):\n",
    "\t\"\"\"Compute annular mode as in Gerber et al, GRL 2008.\n",
    "\t\tThis is basically the first PC, but normalized to unit variance and zero mean.\n",
    "\t\tTo conform to Gerber et al (2008), `data` should be anomalous height or zonal wind\n",
    "\t\t with respect to 30-day smoothed day of year climatology.\n",
    "\t\tINPUTS:\n",
    "\t\t\tlat    - latitude\n",
    "\t\t\tpres   - pressure\n",
    "\t\t\tdata   - variable to compute EOF from. This is typically\n",
    "\t\t\t\t\t\tgeopotential or zonal wind.\n",
    "\t\t\t\t\t\tSize time x pres x lat (ie zonal mean)\n",
    "\t\t\tchoice - not essential, but used for sign convention.\n",
    "\t\t\t\t\t\tIf 'z', the sign is determined based on 70-80N/S.\n",
    "\t\t\t\t\t\tOtherwise, 50-60N/S is used.\n",
    "\t\t\themi   - hemisphere to consider\n",
    "\t\t\t\t\t\t'infer' - if mean(lat)>=0 -> NH, else SH\n",
    "\t\t\t\t\t\t'SH' or 'NH'\n",
    "\t\t\tdetrend- detrend method for computing EOFs:\n",
    "\t\t\t\t\t\t'linear' -> remove linear trend\n",
    "\t\t\t\t\t\t'constant' -> remove total time mean\n",
    "\t\t\teof_in - if None, compute EOF1 as usual.\n",
    "\t\t\t\t\t if the EOF1 is already known, use this instead of\n",
    "\t\t\t\t    computing it again.\n",
    "\t\t\tpc_in  - if None, standardize PC1 to its own mean and std deviation\n",
    "\t\t\t\t else, use pc_in mean and std deviation to standardize.\n",
    "\t\t\teof_out- whether or not to pass the first EOF as output [False].\n",
    "\t\t\tpc_out - whether or not to pass the first PC as output [False].\n",
    "\t\tOUTPUT:\n",
    "\t\t\tAM     - The annular mode, size time x pres\n",
    "\t\t\tEOF    - The first EOF (if eof_out is True), size pres x lat\n",
    "\t\t\tPC     - The first PC (if pc_out is True). size time x pres\n",
    "\t\"\"\"\n",
    "\t#\n",
    "\tAM = np.full((data.shape[0],data.shape[1]),np.nan)\n",
    "\tif pc_out:\n",
    "\t\tpco = np.full(AM.shape,np.nan)\n",
    "\t# guess the hemisphere\n",
    "\tif hemi == 'infer':\n",
    "\t\tif np.mean(lat) >= 0:\n",
    "\t\t\tsgn = 1.\n",
    "\t\telse:\n",
    "\t\t\tsgn = -1.\n",
    "\telif hemi == 'SH':\n",
    "\t\tsgn = -1.\n",
    "\telif hemi == 'NH':\n",
    "\t\tsgn = 1.\n",
    "\tj_tmp = np.where(sgn*lat > 20)[0]\n",
    "\tif eof_out:\n",
    "\t\teofo = np.full((data.shape[1],len(j_tmp)),np.nan)\n",
    "\tcoslat = np.cos(np.deg2rad(lat))\n",
    "\tnegCos = (coslat < 0.)\n",
    "\tcoslat[negCos] = 0.\n",
    "\t# weighting as in Gerber et al GRL 2008\n",
    "\tsqrtcoslat = np.sqrt(coslat[j_tmp])\n",
    "\t# try to get the sign right\n",
    "\t# first possibility\n",
    "\tif choice == 'z':\n",
    "\t\tminj = min(sgn*70,sgn*80)\n",
    "\t\tmaxj = max(sgn*80,sgn*70)\n",
    "\t\tsig = -1\n",
    "\telse:\n",
    "\t\tminj = min(sgn*50,sgn*60)\n",
    "\t\tmaxj = max(sgn*60,sgn*50)\n",
    "\t\tsig = 1\n",
    "\tjj = (lat[j_tmp] > minj)*(lat[j_tmp] < maxj)\n",
    "\t# second possibility\n",
    "\t#jj = abs(lat[j_tmp]-80).argmin()\n",
    "\t#sig = -1\n",
    "\tif isinstance(pres,(int,float)):\n",
    "\t\tdata = np.reshape(data,(data.shape[0],1,data.shape[1]))\n",
    "\t\tpres = [pres]\n",
    "\tE_ls = []\n",
    "\tfor k, p in enumerate(pres):\n",
    "\t\t# remove global mean\n",
    "\t\tglobZ = GlobalAvg(lat,data[:,k,:],axis=-1,lim=lat[j_tmp[0]],mx=lat[j_tmp[-1]])\n",
    "\t\tvar = data[:,k,:] - globZ[:,np.newaxis]\n",
    "\t\t# area weighting: EOFs are ~variance, thus take sqrt(cos)\n",
    "\t\tvar = var[:,j_tmp]*sqrtcoslat[np.newaxis,:]\n",
    "\t\tvarNan = np.isnan(var)\n",
    "\t\tif np.sum(np.reshape(varNan,(np.size(varNan),)))==0:\n",
    "\t\t\tif eof_in is None:\n",
    "\t\t\t\teof1,pc1,E,u,s,v = eof(var,n=1,detrend=detrend)\n",
    "\t\t\t\tprint(p, E)\n",
    "\t\t\t\tE_ls.append(E)\n",
    "\t\t\telse:\n",
    "\t\t\t\tpc1 = eof(var,n=1,detrend=detrend,eof_in=eof_in[k,:])\n",
    "\t\t\t\teof1 = eof_in[k,:]\n",
    "\t\t\t# force the sign of PC\n",
    "\t\t\tpc1  = pc1*sig*np.sign(eof1[jj].mean())\n",
    "\t\t\tif eof_out:\n",
    "\t\t\t\teofo[k,:] = np.squeeze(eof1)\n",
    "\t\t\tif pc_out:\n",
    "\t\t\t\tpco[:,k] = pc1\n",
    "\t\t\t# force unit variance and zero mean\n",
    "\t\t\tif pc_in is None:\n",
    "\t\t\t\tAM[:,k] = (pc1-pc1.mean())/np.std(pc1)\n",
    "\t\t\telse:\n",
    "\t\t\t\tAM[:,k] = (pc1-pc_in.mean())/np.std(pc_in)\n",
    "\tif eof_out and pc_out:\n",
    "\t\treturn AM,eofo,pco\n",
    "\telif eof_out:\n",
    "\t\treturn AM,eofo\n",
    "\telif pc_out:\n",
    "\t\treturn AM,pco\n",
    "\telse:\n",
    "\t\treturn AM,E_ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4aa12e52-2529-4abc-884c-78d358a96ae4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-03-04T15:58:23.008656Z",
     "iopub.status.busy": "2022-03-04T15:58:23.008175Z",
     "iopub.status.idle": "2022-03-04T15:58:23.015499Z",
     "shell.execute_reply": "2022-03-04T15:58:23.014500Z",
     "shell.execute_reply.started": "2022-03-04T15:58:23.008620Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "root_path = '/mnt/data04/model/CMAM/0A.daily/'\n",
    "label_ls = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I','J','K','L']\n",
    "its = 10000\n",
    "plt.rcParams.update({'font.size': 22})\n",
    "sns.set_style('whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f05d63fa-5660-4963-ba70-462640561749",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-03-04T15:58:25.942164Z",
     "iopub.status.busy": "2022-03-04T15:58:25.941634Z",
     "iopub.status.idle": "2022-03-04T15:58:26.139146Z",
     "shell.execute_reply": "2022-03-04T15:58:26.138086Z",
     "shell.execute_reply.started": "2022-03-04T15:58:25.942122Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
       "Dimensions:  (lat: 48, plev: 63, time: 46720)\n",
       "Coordinates:\n",
       "  * plev     (plev) float64 1e+03 925.0 850.0 800.0 ... 0.0015 0.001 0.0007\n",
       "  * lat      (lat) float64 -87.16 -83.48 -79.78 -76.07 ... 79.78 83.48 87.16\n",
       "  * time     (time) datetime64[ns] 1997-01-01 ... 1990-12-31T18:00:00\n",
       "    month    (time) int64 1 1 1 1 1 1 1 1 1 1 ... 12 12 12 12 12 12 12 12 12 12\n",
       "Data variables:\n",
       "    ua       (time, plev, lat) float64 ...</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-b2c3109e-07a1-4719-9f64-b50fc0ea0122' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-b2c3109e-07a1-4719-9f64-b50fc0ea0122' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>lat</span>: 48</li><li><span class='xr-has-index'>plev</span>: 63</li><li><span class='xr-has-index'>time</span>: 46720</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-9df7892a-26c2-498c-9b16-5339b5f2cc03' class='xr-section-summary-in' type='checkbox'  checked><label for='section-9df7892a-26c2-498c-9b16-5339b5f2cc03' class='xr-section-summary' >Coordinates: <span>(4)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>plev</span></div><div class='xr-var-dims'>(plev)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>1e+03 925.0 850.0 ... 0.001 0.0007</div><input id='attrs-1782ffb4-0073-4032-9cdd-399b87412cfa' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-1782ffb4-0073-4032-9cdd-399b87412cfa' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-5c7f6032-951c-4bb2-b145-9c2f1f60d63f' class='xr-var-data-in' type='checkbox'><label for='data-5c7f6032-951c-4bb2-b145-9c2f1f60d63f' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>hPa</dd><dt><span>long_name :</span></dt><dd>pressure</dd></dl></div><div class='xr-var-data'><pre>array([1.00e+03, 9.25e+02, 8.50e+02, 8.00e+02, 7.80e+02, 7.50e+02, 7.00e+02,\n",
       "       6.50e+02, 6.00e+02, 5.00e+02, 4.50e+02, 4.00e+02, 3.50e+02, 3.00e+02,\n",
       "       2.85e+02, 2.50e+02, 2.00e+02, 1.70e+02, 1.50e+02, 1.30e+02, 1.15e+02,\n",
       "       1.00e+02, 9.00e+01, 8.00e+01, 7.00e+01, 6.00e+01, 5.00e+01, 4.00e+01,\n",
       "       3.50e+01, 3.00e+01, 2.50e+01, 2.00e+01, 1.50e+01, 1.00e+01, 7.00e+00,\n",
       "       5.00e+00, 4.00e+00, 3.00e+00, 2.00e+00, 1.50e+00, 1.00e+00, 7.00e-01,\n",
       "       5.00e-01, 4.00e-01, 3.00e-01, 2.00e-01, 1.50e-01, 1.00e-01, 7.00e-02,\n",
       "       5.00e-02, 4.00e-02, 3.00e-02, 2.00e-02, 1.50e-02, 1.00e-02, 7.00e-03,\n",
       "       5.00e-03, 4.00e-03, 3.00e-03, 2.00e-03, 1.50e-03, 1.00e-03, 7.00e-04])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lat</span></div><div class='xr-var-dims'>(lat)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>-87.16 -83.48 ... 83.48 87.16</div><input id='attrs-2af05570-63e4-43f2-aa72-924bf3668119' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-2af05570-63e4-43f2-aa72-924bf3668119' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-75a931d1-bc86-4ac2-833e-b3e4d099b743' class='xr-var-data-in' type='checkbox'><label for='data-75a931d1-bc86-4ac2-833e-b3e4d099b743' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>bounds :</span></dt><dd>lat_bnds</dd><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>axis :</span></dt><dd>Y</dd><dt><span>long_name :</span></dt><dd>latitude</dd><dt><span>standard_name :</span></dt><dd>latitude</dd></dl></div><div class='xr-var-data'><pre>array([-87.159094, -83.478937, -79.777046, -76.070244, -72.361581, -68.652017,\n",
       "       -64.941949, -61.231573, -57.520994, -53.810274, -50.099453, -46.388558,\n",
       "       -42.677606, -38.96661 , -35.25558 , -31.544523, -27.833444, -24.122348,\n",
       "       -20.411238, -16.700118, -12.988988,  -9.277853,  -5.566714,  -1.855571,\n",
       "         1.855571,   5.566714,   9.277853,  12.988988,  16.700118,  20.411238,\n",
       "        24.122348,  27.833444,  31.544523,  35.25558 ,  38.96661 ,  42.677606,\n",
       "        46.388558,  50.099453,  53.810274,  57.520994,  61.231573,  64.941949,\n",
       "        68.652017,  72.361581,  76.070244,  79.777046,  83.478937,  87.159094])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>1997-01-01 ... 1990-12-31T18:00:00</div><input id='attrs-3b5540d0-33a0-4c94-b56c-4b90cafe572c' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-3b5540d0-33a0-4c94-b56c-4b90cafe572c' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-bef853e1-bc21-44f3-926b-4a7645056c68' class='xr-var-data-in' type='checkbox'><label for='data-bef853e1-bc21-44f3-926b-4a7645056c68' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>bounds :</span></dt><dd>time_bnds</dd><dt><span>axis :</span></dt><dd>T</dd><dt><span>long_name :</span></dt><dd>time</dd><dt><span>standard_name :</span></dt><dd>time</dd></dl></div><div class='xr-var-data'><pre>array([&#x27;1997-01-01T00:00:00.000000000&#x27;, &#x27;1997-01-01T06:00:00.000000000&#x27;,\n",
       "       &#x27;1997-01-01T12:00:00.000000000&#x27;, ..., &#x27;1990-12-31T06:00:00.000000000&#x27;,\n",
       "       &#x27;1990-12-31T12:00:00.000000000&#x27;, &#x27;1990-12-31T18:00:00.000000000&#x27;],\n",
       "      dtype=&#x27;datetime64[ns]&#x27;)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>month</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-9164f6bc-0e13-4acd-b705-4de9459c6249' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-9164f6bc-0e13-4acd-b705-4de9459c6249' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-dce04113-50eb-4f65-a899-526a2a8ba424' class='xr-var-data-in' type='checkbox'><label for='data-dce04113-50eb-4f65-a899-526a2a8ba424' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([ 1,  1,  1, ..., 12, 12, 12])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-805b345d-e785-4f2e-be62-5db8d62e02f2' class='xr-section-summary-in' type='checkbox'  checked><label for='section-805b345d-e785-4f2e-be62-5db8d62e02f2' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>ua</span></div><div class='xr-var-dims'>(time, plev, lat)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-f01d92d9-1ef1-41c9-9e0b-5cd9759cbde3' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-f01d92d9-1ef1-41c9-9e0b-5cd9759cbde3' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-93135aae-9ee7-41e8-9048-9ce02c487721' class='xr-var-data-in' type='checkbox'><label for='data-93135aae-9ee7-41e8-9048-9ce02c487721' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>[141281280 values with dtype=float64]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-5dd39036-9c9d-47b3-9c18-e78cff8f743d' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-5dd39036-9c9d-47b3-9c18-e78cff8f743d' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lat: 48, plev: 63, time: 46720)\n",
       "Coordinates:\n",
       "  * plev     (plev) float64 1e+03 925.0 850.0 800.0 ... 0.0015 0.001 0.0007\n",
       "  * lat      (lat) float64 -87.16 -83.48 -79.78 -76.07 ... 79.78 83.48 87.16\n",
       "  * time     (time) datetime64[ns] 1997-01-01 ... 1990-12-31T18:00:00\n",
       "    month    (time) int64 ...\n",
       "Data variables:\n",
       "    ua       (time, plev, lat) float64 ..."
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = xr.open_dataset(f'{root_path}ua/ua_6hrPlev_CMAM_CMAM30-SD_r1i1p1_19790101-20101231_anomalies_zm.nc')#.dropna('time', how='all')\n",
    "#ds = ds.where(np.isfinite(ds), drop = True)\n",
    "ds['plev'] = ds.plev/100.\n",
    "ds['plev'].attrs['units'] = 'hPa'\n",
    "ds['plev'].attrs['long_name'] = 'pressure'\n",
    "ds_sel = ds.sel(lat = slice(0,None))\n",
    "\"\"\"mask =  (ds.time.dt.month == 12) | \\\n",
    "    (ds.time.dt.month == 1) | (ds.time.dta.month == 2) #| (ds.time.dt.month == 3)\n",
    "\n",
    "ds_sel = ds_sel.where(np.isfinite(ds_sel), drop = True)\n",
    "#ds['time'] = xr.cftime_range('1979-01-01T06', '2010-12-31T18', freq='6H', calendar='noleap')\n",
    "ds_sel\"\"\"\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9c9e87b2-9841-497c-bf6f-3f818d79e224",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-03-04T15:58:30.862903Z",
     "iopub.status.busy": "2022-03-04T15:58:30.862395Z",
     "iopub.status.idle": "2022-03-04T15:59:13.304409Z",
     "shell.execute_reply": "2022-03-04T15:59:13.303095Z",
     "shell.execute_reply.started": "2022-03-04T15:58:30.862865Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000.0 [0.28495444]\n",
      "925.0 [0.30510426]\n",
      "850.0 [0.32292304]\n",
      "800.0 [0.32680078]\n",
      "780.0 [0.32714974]\n",
      "750.0 [0.32704052]\n",
      "700.0 [0.32531014]\n",
      "650.0 [0.32204996]\n",
      "600.0 [0.31742462]\n",
      "500.0 [0.30552052]\n",
      "450.0 [0.29916022]\n",
      "400.0 [0.29359553]\n",
      "350.0 [0.29007271]\n",
      "300.0 [0.29048437]\n",
      "285.0 [0.29121662]\n",
      "250.0 [0.298128]\n",
      "200.0 [0.32699982]\n",
      "170.0 [0.36207709]\n",
      "150.0 [0.39198185]\n",
      "130.0 [0.42966897]\n",
      "115.0 [0.46788942]\n",
      "100.0 [0.5194843]\n",
      "90.0 [0.56116193]\n",
      "80.0 [0.60726204]\n",
      "70.0 [0.65351635]\n",
      "60.0 [0.6937017]\n",
      "50.0 [0.72235427]\n",
      "40.0 [0.73607915]\n",
      "35.0 [0.73605411]\n",
      "30.0 [0.73239803]\n",
      "25.0 [0.72295168]\n",
      "20.0 [0.71103766]\n",
      "15.0 [0.69501598]\n",
      "10.0 [0.68500942]\n",
      "7.0 [0.67964719]\n",
      "5.0 [0.68075764]\n",
      "4.0 [0.67947005]\n",
      "3.0 [0.67637964]\n",
      "2.0 [0.66802664]\n",
      "1.5 [0.66233281]\n",
      "1.0 [0.65697891]\n",
      "0.7 [0.64779029]\n",
      "0.5 [0.63269287]\n",
      "0.4 [0.62212502]\n",
      "0.3 [0.60574766]\n",
      "0.2 [0.58338721]\n",
      "0.15 [0.56965028]\n",
      "0.1 [0.55277931]\n",
      "0.07 [0.54189023]\n",
      "0.05 [0.53353867]\n",
      "0.04 [0.53167376]\n",
      "0.03 [0.52109934]\n",
      "0.02 [0.5092943]\n",
      "0.015 [0.5025293]\n",
      "0.01 [0.5083152]\n",
      "0.007000000000000001 [0.51548219]\n",
      "0.005 [0.51512212]\n",
      "0.004 [0.50234543]\n",
      "0.003 [0.4720415]\n",
      "0.002 [0.47798642]\n",
      "0.0015 [0.48055426]\n",
      "0.001 [0.48162979]\n",
      "0.0006999999999999999 [0.40635839]\n"
     ]
    }
   ],
   "source": [
    "am_arr, E_arr = ComputeAnnularMode(ds_sel.lat.values, ds_sel.plev.values, ds_sel.ua.values, \\\n",
    "                            choice = 'u')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a290ebc5-d5de-44be-9727-1749c78e5eab",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-03-04T15:59:31.022680Z",
     "iopub.status.busy": "2022-03-04T15:59:31.022255Z",
     "iopub.status.idle": "2022-03-04T15:59:31.028244Z",
     "shell.execute_reply": "2022-03-04T15:59:31.027590Z",
     "shell.execute_reply.started": "2022-03-04T15:59:31.022647Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(46720, 63)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "am_arr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3de0e0b7-b7b8-44b6-800e-4643ed6cac91",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-03-04T16:00:49.198804Z",
     "iopub.status.busy": "2022-03-04T16:00:49.198195Z",
     "iopub.status.idle": "2022-03-04T16:00:50.821879Z",
     "shell.execute_reply": "2022-03-04T16:00:50.821034Z",
     "shell.execute_reply.started": "2022-03-04T16:00:49.198731Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAElCAYAAABTQG2/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABV/ElEQVR4nO3deVxU1fsH8M8M+46yCSgiyyACIuIuLiHmirsWfQU1t9JKLTX72bevWGaLSZhWZrkjZZoLaoq7UihuICICsi+C7DDsMPf3B86FkcUBZpiF5/168WLm3HPnPnMGeDj3nnsOh2EYBoQQQogC4co6AEIIIaStKHkRQghROJS8CCGEKBxKXoQQQhQOJS9CCCEKh5IXIYQQhaMq6wAURWpqKn777TdERUUhISEBNjY2OHPmjKzDIoSQLomSl5gSEhJw/fp1uLq6QiAQgG6PI4QQ2eHQTcriEQgE4HLrz7Ju2LABjx49op4XIYTICF3zEpMwcRFCCJE9uTltmJSUhJs3byI6OhqPHj1CSkoKGIZBYGAgJk6c2Oq+ISEhCA4ORlxcHAQCAfr06YPZs2fDx8eHkg4hhCghuUlewcHBOHjwYJv38/f3x5EjR6ChoYHhw4dDVVUV4eHh2Lx5M8LDwxEYGAgVFRUpREwIIURW5CZ58Xg8LF68GM7OznB2dsbGjRsRERHR6j4XLlzAkSNHYGJigsOHD8Pa2hoAkJeXBz8/P1y8eBGHDx/GggULRPYrLS3F8+fPXxmThYUFtLS02v2eCCGESIfcJK+5c+e2eZ/du3cDANauXcsmLgAwNjbGpk2b4Ovriz179sDX11fk9OHFixfxySefvPL19+3bhxEjRrQ5LkIIIdIlN8mrrbKzsxETEwM1NbVmr4kNGTIEZmZmyMnJQWRkJAYOHMhumzVrFmbNmtWZ4RJCCJEghR3N8PjxYwCAvb09NDU1m63j4uICAIiNje20uAghhEifwva8MjIyANRfl2qJubm5SN2OqKiowPXr1wEAmZmZ4PP5OH/+PID6JGlpadnq/vfv3weHw+lwHK1pfC8aofZoTB7aoqxagKySGghvLFXlAhZ6atBU69y45KEt5Ik8tweHw8GAAQOa3aawyau8vBwAWh1QoaOjAwAoKyvr8PHy8/OxatUqkTLh861bt77yNKSWlhYcHR07HEdrYmNjpX4MRULt0UBe2uLYvQx8fPwh6gT1KUxTjYuAeQMwycW802KQl7aQF/LcHq2dNVPY5CWcGETavRmhnj17Ii4urlOORYiymuPeE+YGmnj38D2UVNaiskaAd4PuY90EB6wYa9tpv89E8clnX1EMwl6VsAfWHGGPS1iXECJ7I+2McWLlSFgbabNl316Iw0dHo1BVWyfDyIgiUdjkJbzGlJWV1WKd7OxskbqEEPlga6KLEytGYphNd7bsrweZeGvPbeTxq2QYGVEUCpu8+vXrB6B+tvfKyspm60RHRwOA3J7PJaQr66ajjoNvD8Wbg3uxZfdSCzFj1z+IzymVYWREEShs8jI3N4eTkxNqamrYUX+NRUREIDs7GyYmJnBzc5NBhISQV1FX5WLrLBd8OsURwstdGYUVmPXjv7ga9+pZcEjXpbDJCwCWLVsGANi2bRtSU1PZ8vz8fPj7+wMAli5dKrfDQAkh9YOuloyywR7fQdBRr5+HlF9Vi8X772DfP8m0dh5pltyMNoyJiWETDgA8ffoUABAQEIC9e/ey5UePHmUfT5w4ET4+PggODoa3tzdGjBjBTszL5/Ph5eWF+fPnd96bIIS0m1c/Mxx7dwSWHLiLzKIKCBjAP+Qxnj7n43/eTlBXpX9CSQO5SV58Ph9RUVFNylNSUlrdb9OmTXB3d0dQUBAiIiIgEAhgY2NDS6IQooAczfVxYuUILD90Dw/SigAAQbfT8CizGDt83NDbiEYOk3pyk7yGDh3a7vuovL294e3tLeGICCGyYKqnieClw7D+2EOcjqofTRyVUYwpO8KwZaYzpg+g0cNEwa95EUKUk6aaCgLfHIBPpzhCTaV+JAe/qharfo/E+mNRKK+ulXGERNYoeRFC5JJwIMfxd0egd6Mbmo/ezYD3D2F4nFUiw+iIrFHyIoTItf49DXHmfQ9MH9AwCXdibhlm/PgPDoan0GjELoqSFyFE7ulpquH7Nwbg2zn9oaVWP5y+ulaAz07FYPmheygqr5ZxhKSzUfIihCgEDoeDuYN64cwHHnA012fLQx/nYHLgTfzzNE+G0ZHORsmLEKJQ6udFHIGFI6zZsqziSvzn19v46GgUCsuoF9YVUPIihCgcTTUVbJrmhF983WGorcaWH7+fgXHbr+Pkg0y6FqbkKHkRQhTW6049ELpmNKb0b1jMsqCsGqv/iMSCfXeQXtDykklEsVHyIoQoNFM9Tex6ayB+9RsECwNNtvxGfC7GB1zHLzcSUVsnkGGERBooeRFClIJXPzOEfjgGC0dYszPUV9YI8OW5J5i+6x9EZxTLNkAiUZS8CCFKQ1dDFZumOeHEipHo20OPLY/JKsH0XWH44sxjVNZQL0wZUPIihCidAb0MEfK+B9ZPdIDGi9noBQzwa1gylp/KwDVaK0zhUfIihCglNRUuVoy1w4XVozHC1ogtf15Wi4X77mDV7w+Qx6+SYYSkIyh5EUKUmrWxDoKWDMW2ua4iw+pPRWZh3HfXcfRuOg2rV0CUvAghSo/D4WCOe09c+nAMXrPRZcuLK2qw/thD+Oy5hafPS2UYIWkrSl6EkC7DWFcD60eZYv+iwejZTYstv5VUgEmBN/HN+SeoqK6TYYREXJS8CCFdzlgHU4SuGY1lo22gwq0fV19Tx+DHa4kYH3AdV57kyDhC8iqUvAghXZK2uir+b7IjzrzvgYFWhmx5RmEF3t5/F+8cuodnxRWyC5C0ipIXIaRLczTXx7F3RuCrWS4iAzrOx2Rj3HfX8evNJJqhQw5R8iKEdHlcLgdvDrHC5Q/HYI57T7a8vLoOX5yNxdQfwnAvtVCGEZKXUfIihJAXjHQ1sG2uK44uHw6eWcOoxCfZpZj907/45K+HtPClnFBtacPJkycldpAZM2ZI7LUIIUTahvTpjrMfjMJvYckIvJSAipr6EYjBEem4HpeLHT5uGGTdXcZRdm0tJq8NGzaAI5zdsoMoeRFCFI2aChfvjLHF1P7m2HT6MS7F1o9AzCquxBu/3MKH43l4d4wtuFzJ/J0kbdNi8gIAKysrDBw4sN0vfu/ePaSnp7d7f0IIkbWe3bTx64JBOP/oGT4+Ho3iihrUCRh8eyEO4Yn52P6GK0z1NF/9QkSiWk1eAwcOxNatW9v94p988gklL0KIUpjobA6XnoZYFfwAd18M3gh7mofJgWEIeMMVo+xNZBxh1yLVARs0XxghRJlYGmrh92XDsPI1W3bNsDx+Ffz2RuDbC09oSH0nanXAhoGBQYdefNWqVVi4cGGHXoMQQuSJqgoX6yb0xXAbY6z+IxJ5/CowDLDraiLuJBdi51tuMNWn04jS1mLPq2/fvjA3N+/Qi5ubm6Nv374deg1CCJFHHvbGOLfKA6PsjdmyiJQCTPkhDHdSCmQYWddA93kRQkg7mepp4sCiIVg3wQHCQYe5pVXw+eUWfgtLpksnUkTJixBCOoDL5WDla3Y4tHgouuuoAwBqBQw+P/MY7wc/QFlVrYwjVE6tjjZsTk5ODi5fvoyUlBTw+fxm/7PgcDj48ssvJRIgIYQogpF2xjjzvgfeDbqPqPQiAMCZh88Ql12Kn33dYWui2/oLkDZpU/Lav38/vvvuO9TWNvwnIUxewhuaGYah5EUI6ZIsDLVwdPkwbA55jKDbaQCAhOd8TN/5D76d0x+TXDo2joA0EPu04c2bN/HVV19BQ0MDy5Ytw4ABAwAAmzdvxuLFi9GzZ/1kln5+fpS4CCFdloaqCrbMdMG2ua7QUK3/E8uvqsW7Qfex9VwsDaeXELGT16FDh8DhcLB3716sWbMG1tbWAIB58+Zh3bp1OHv2LObMmYNjx47Bzc1NWvESQohCmOPeE3+tGIFe3RtWbN59Iwnzf7uN3NIqGUamHMROXtHR0XB2dkb//v2b3a6uro7PPvsMurq62Llzp8QClCcpKSlYvHgx3NzcMGzYMHz++eeoqKDF6gghzXOyMMCZ90bBs68pW3YrqQAzdv2DfD4lsI4QO3mVlpaiV69e7HM1tfpF28rLy0XKBg4ciNu3b0swRPlQUlICPz8/lJWVITAwEBs2bMCZM2fwf//3f7IOjRAixwy01fCr3yB8NJ7HzsqRWVSBtX9G0VD6DhA7eXXr1g18Pp99Lpx9IzMzU6ReVVUVSkpKJBSe/Pj9999RUlKCH3/8EaNHj8aMGTPw6aef4ty5c0hISJB1eIQQOcblcvD+OHv8+FbDROdX43LxW1iyDKNSbGInL0tLS2RlZbHPHR0dwTAMzpw5w5bl5+cjIiICFhYWko1SDty4cQPDhg1D9+4Na/hMmDAB6urquHHjhgwjI4Qoikku5lg6qg/7/OvzT9hh9aRtxB4qP3z4cPz888/IysqChYUFxowZAwMDA/zyyy9ITU1Fjx49EBoaivLycnh5eUkzZlZSUhJu3ryJ6OhoPHr0CCkpKWAYBoGBgZg4cWKr+4aEhCA4OBhxcXEQCATo06cPZs+eDR8fH3C5TXN6YmIiZs+eLVKmrq4OKysrJCUlSfR9EUKU17oJfRGRXICojGLU1DF4P/gBzn7gAT1NNVmHplDETl5Tp05Fbm4um7x0dHSwZcsWrF27FufPn2fr9evXD++++65Ugn1ZcHAwDh482Ob9/P39ceTIEWhoaGD48OFQVVVFeHg4Nm/ejPDwcAQGBkJFRUVkn5KSEujr6zd5LX19fRQXF7f7PRBCuhZ1VS5+8BmIKTtuorSqFmkF5fi/E4+w480BElsAuCsQO3nZ2triiy++ECnz8vLChQsXcPXqVRQXF8PGxgaenp5N/vBLC4/Hw+LFi+Hs7AxnZ2ds3LgRERERre5z4cIFHDlyBCYmJjh8+DA75D8vLw9+fn64ePEiDh8+jAULFogVg/CmbEIIEZeVkTa+nOWC94MfAABCorLgYWeENwZbyTgyxfHK5CUQCBAdHY3MzEyoq6vD0dERlpaW7HYzMzO8+eabUg2yJXPnzm3zPrt37wYArF27lk1cAGBsbIxNmzbB19cXe/bsga+vr8jpQ319/WYHopSWlsLW1rbtwRNCujRvVwv8m5iH4Ij6BXv/dzoGblbdwDPTk3FkiqHV5HX37l1s2LChyYjCcePGYevWrdDTU6xGzs7ORkxMDNTU1Jq9JjZkyBCYmZkhJycHkZGRGDiwYWSQra0tEhMTRepXV1cjLS0Ns2bNknrshBDl89lUJ9xLLUR8Dh+VNQK8d+Q+Tr/nAU21zjl7pchaHG2Ynp6OpUuXIiMjAwzDQF9fH2pqamAYBpcvX8bq1as7MUzJePz4MQDA3t4emprNLxbn4uICAIiNjRUpHz16NG7duoXCwkK27OLFi6iursaYMWOkFDEhRJlpqatg51sDoalW/6c4PoePP++myzgqxdBiz+u3335DRUUFRowYAX9/f/Tq1QsCgQDXrl3Dp59+in///RcPHz5sccYNeZSRkQEArQ7lFy7AKawr9Oabb+Lw4cNYsWIFVqxYgfz8fHz11VeYPHky7OzsXnnsqqqqJglR0iorK6V+DEVC7dGA2qKBPLbFfFdD/Hq3fgHLg2FPMciw82bukcf2EEeLyevWrVswMDBAYGAge3qQy+XC09MT69evx4YNG3Dr1i2FSl7C2UC0tLRarKOjowMAKCsrEynX19fHgQMH8MUXX+D999+HhoYGpkyZgnXr1ol1bA0NDTg6OrYzcvHExsZK/RiKhNqjAbVFA3lsi/d6V+Pgg8uorhMgIb8K3G494dCjcy7LyGN7CLWWVFtMXtnZ2Rg8eHCz17VGjx4NAHj27JkEwus8Ly/f0lZ9+vTBb7/9JsmQCCEEhtrqGO9khrMP6/+mHruXjo1T+sk4KvnW4jWvyspKGBsbN7tNOMtEVZViTSwp7FU1no/xZcIel7AuIYR0hrnuPdnHJx5kooaWTmmV2NNDNUfRJpUUDvFvPM3Vy7Kzs0XqEkJIZxhlbwIzfQ0AQB6/GtfjcmUckXxrdah8eXl5q3/oW9suj/Mb9utX3w1PSEhAZWVlsyMOo6OjAUBuzwETQpSTCpeDWQN74qdr9bfk/HkvHV79zGQclfxqNXmFhoYiNDS02W0cDqfF7RwOhx2WLk/Mzc3h5OSEmJgYnD9/HjNmzBDZHhERgezsbJiYmNCCmoSQTjfHvSF5XY59jnx+FYx0NWQclXxq9bQhwzDt+hII5Pdc7bJlywAA27ZtQ2pqKluen58Pf39/AMDSpUubnZyXEEKkydZEFwOtDAEAtQIGpyJbPvPV1bXY83ry5ElnxtEuMTExbMIBgKdPnwIAAgICsHfvXrb86NGj7OOJEyfCx8cHwcHB8Pb2xogRI9iJefl8Pry8vDB//vzOexOEENLI3EG9cD+tCABw9G46Fo20pvlTmyH2xLzyiM/nIyoqqkl5SkpKq/tt2rQJ7u7uCAoKQkREBAQCAWxsbFpdEoUQQjrDlP7m8A+JQWWNAE+yS3E6KgvTB9AAspcpdPIaOnQo4uLi2rWvt7c3vL29JRwRIYR0jL6mGhaMsMbu6/XrBH799xNMcOpB8x2+hLoYhBAiZ957zQ5GOuoAgKziSvx6kxa8fVmbel4ZGRn45ZdfEB4ejufPn6O6urrZevI62pAQQhSBnqYaPnydh40nHgEAfryWiHmDe8FUr/kJxbsisXteCQkJmDlzJv7880+kp6ejqqpKIUcbEkKIInhjUC/wzHQBAOXVdfjuQryMI5IvYve8AgICUFpaijFjxmDlypWwsbGBrq6uNGMjhJAuS1WFi0+n9IPf3vrV4Y/eS4ffiN5wsjCQcWTyQeye1927d2FpaYmdO3eif//+lLgIIUTKRvNM8JqDCQCAYYDPzzxWuGn5pEXs5FVdXQ0XFxeoqalJMx5CCCGNbJziCBVu/X1et5IKcC46W8YRyQexk5e1tTVKS0ulGQshhJCX2JnqYcFwa/b5lrOPUV5dK7uA5ITYyWvu3Lm4e/dukxWGCSGESNfq8fYw1m0YOi+c/7ArEzt5/ec//8GkSZOwaNEiXL9+nUYUEkJIJ9HXVMP6CX3Z57tvJCEtv+V1CbuCFkcbjhs3rtnyzMxMvPPOO1BRUYGpqWmzc25xOBxcunRJclESQkgXN8e9J4Ii0hCVXoTqWgE+P/sYe/wGyTosmWkxeWVmZra4E8MwqK2tbXEtL5pEkhBCJIvL5cB/mhNm7PoHAHDxcQ6ux+diDM9ExpHJRovJ6/Lly50ZByGEkFcY0MsQ8wb1xNG79WMP/E/H4Pzq0VBX7Xoz/bWYvCwtaRZjQgiRN+sm9MXf0dkorapFUl4Z9v2TjOVjbGUdVqfreumaEEIUmImeBlaP57HPd1xOwPOSShlGJBuUvAghRMH4De8Ne9P6WY7Kqutw+FbqK/ZQPi0mr3feeQcHDx7s0IsfOHAA77zzTodegxBCiCg1FS5Wedmzz0Mf58gwGtlo8ZrXtWvX0K1btw69+JMnT3D9+vUOvQYhhJCmxjqYQl2Fi+q6+hWX0wvK0au7tqzD6jStzipfXl7e4nB4cZSVlbV7X0IIIS3T1VDFCDsjXIvLBQBcjs3BwpF9ZBxV52k1eYWGhiI0NLSzYiGEENIGXo5mbPK6SMmrnoWFRWfGQQghpI3G9zPDpyfrV1u+nVSA4ooaGGh1jZU/WkxeV65c6cw4CCGEtJGZviZcexogKqMYtQIG1+KeY/qArnGPLg2VJ4QQBTbO0Yx9fCupQIaRdC5KXoQQosBcexmyjxNyus6ai5S8CCFEgQlvVgaAhOd8MAwjw2g6DyUvQghRYOYGmtDVqB++UFxRg1x+lYwj6hyUvAghRIFxOBzYNep9Pc3hyzCazkPJixBCFJytSUPySsrrGpNDUPIihBAFVycQsI811VRkGEnnoeRFCCEKrvF1LhM9DRlG0nnanbyqq6vx/PlzFBUVSTAcQgghbZVb2ih56VLyatbJkycxe/ZsuLm5YcyYMfj666/ZbefPn8dHH32E9PR0iQZJCCGkZc8bJS9T/a6RvFqdmPdlGzZswKlTp8AwDLS1tVFeXi6y3dzcHGfPnoWjoyOWLFki0UAJIYQ0VVVbh6LyGgCACpeDbtrqMo6oc4jd8zpx4gROnjyJvn374vjx47h3716TOq6urjA1NcWNGzckGiQhhJDm5fOr2cdGOupQ4XJkGE3nEbvndfToUejo6ODnn3+GmZlZi/WsrKyQmZkpkeDkSWpqKn777TdERUUhISEBNjY2OHPmjKzDIoR0cSLXu7rIYA2gDT2v+Ph4DBgwoNXEBQCmpqbIy8vrcGDyJiEhAdevX0fv3r1ha2sr63AIIQRA/awaQpyu0ekC0IbkVVtbC23tVy8xXVRUBFXVNl1KUwienp64fv06duzYAScnJ1mHQwghAIC+5npQfXGq8FFmCWKyimUcUecQO3mZm5sjPj6+1Tp1dXVISEiAlZVVhwOTN1wu3RJHCJE/pnqamORizj4/8G+K7ILpRGL/Rfbw8EBaWhpOnTrVYp3ff/8dubm5GDNmjESCe5WkpCQcOHAAa9euxcSJE9G3b184ODjg/Pnzr9w3JCQEb731Ftzd3eHm5oZZs2YhKCgIgkZ3qhNCiCJYOMKafXwyMgsFZdUtV1YSYp/fW7JkCU6ePImNGzciMTEREyZMAFB/s3JiYiL+/vtv7N69G4aGhvD19ZVawI0FBwfj4MGDbd7P398fR44cgYaGBoYPHw5VVVWEh4dj8+bNCA8PR2BgIFRUusYUK4QQxTfQyhD9exrgYUYxqmsFCI5Iw8rX7GQdllSJ3fPq0aMHdu7cCW1tbezZswdz5swBh8PBuXPnMHXqVOzcuRMaGhrYsWMHjIyMpBkzi8fjYfHixQgICMDFixcxZMiQV+5z4cIFHDlyBCYmJjh9+jR2796NXbt2ITQ0FLa2trh48SIOHz7cCdETQohkcDgckd7X4VupqKlT7rNIbbqQM2zYMJw9exZvv/027OzsoKmpCTU1NVhZWWH+/Pk4c+aMWAlEUubOnYv169dj8uTJYl9n2717NwBg7dq1sLa2ZsuNjY2xadMmAMCePXvo9CEhRKFM6W8OY936G5SfFVciNCZHxhFJl9inDZ88eQIulwsej4d169Zh3bp10oxLKrKzsxETEwM1NTVMnDixyfYhQ4bAzMwMOTk5iIyMxMCBA2UQJSGEtJ2GqgreGtobOy4nAAD2/5uMKf3NX7GX4hK75zVjxgx8/vnn0oxF6h4/fgwAsLe3h6amZrN1XFxcAACxsbGdFhchhEjC/KFW7LD5OymF+DdR+e65FRK756Wvr//KG5TlXUZGBgDAwsKixTrm5uYidYUqKipw/fp1AEBmZib4fD47qtHFxQWWlpatHruqqkrqCbGyspKSbiPUHg2oLRooe1uMstbB1aT61ZRXHbmLH6f1hK56ywPQFLU9xE5effv2VfjZ4oUTCWtpabVYR0dHBwBQVia6Gml+fj5WrVolUiZ8vnXrVsyaNavVY2toaMDR0bHNMbdFbGys1I+hSKg9GlBbNFD2tthq0QcTA2+gqLwGuWV1OBxbg8A3nVusL8/t0VpSFTt5+fn54b333sONGzcwevRoiQTW2RiGAVA/Mqetevbsibi4OEmHRAghEtXDQBNbZ7rg3aD7AIBTkVnw7GuK6QNaPzukaMROXv369cP8+fOxYsUKzJ49G+PHj4eFhUWL145aOzUnK8Je1ctLuTQm7HEJ6xJCiKKZ5GKOue498ee9+ssfn558BPfe3dCz26un+FMUYievcePGAajvvRw9ehRHjx5tsS6Hw2EHR8gT4XWprKysFutkZ2eL1CWEEEX0v2lOuJ1cgLSCcpRW1uLDo1EIXjpMaZZMETt5CQcyKLJ+/foBqJ8hvrKystleY3R0NADI7TlgQggRh66GKgLeGIB5u8NRJ2AQkVyAX24k4d2xyrEqhtjJ68qVK9KMo1OYm5vDyckJMTExOH/+PGbMmCGyPSIiAtnZ2TAxMYGbm5tsgiSEEAlx790N771mh8AX935tvxiHUfbGcLY0kHFkHdflpkpftmwZAGDbtm1ITU1ly/Pz8+Hv7w8AWLp0Kc0iTwhRCu972mFAL0MAQE0dg1W/P0BFdZ1sg5IAhV54KyYmhk04APD06VMAQEBAAPbu3cuWN74+N3HiRPj4+CA4OBje3t4YMWIEOzEvn8+Hl5cX5s+f33lvghBCpEhVhYvv3xiAyTtuory6Dom5Zfjkr4fYOqs/tFq5/0veiZ28Whvk0JzOGG3I5/MRFRXVpDwlJaXV/TZt2gR3d3cEBQUhIiICAoEANjY2mD17Nnx8fKjXRQhRKtbGOtjk7YT1xx8CqF825V5aIbbMcIGJjGNrL7GTl6enp9j3R3XWaMOhQ4e2+94rb29veHt7SzgiQgiRT3MH9UR4Uj5OPMgEAKQXVMBvbwTG9tHFtz1tYKKnIeMI20bs5NVST4phGOTm5qK2trbVeoQQQmSHw+Hgu7muGGbTHV+ee4LiihoAwLVkPry2X8cnk/pi3qBe4CrIUHqJjDasra1FWFgYvvjiC7i7u+Prr7+WSHCEEEIkh8vl4I3BVhjnaIbPzzzGqcj6y0HFFTXY8Fc0/rqfiS9nOcPOVE/Gkb6aRC7uqKqqYuzYsfjpp59w7tw5BAUFSeJlCSGESIGxrgYC33TDgbeHoIduQx8mIqUAkwJvYvvFeFTWyPeIRImOTLC3t4eLiwv++OMPSb4sIYQQKRjDM8FP03vi3bG27FIqNXUMdlxOwOTAm3K9pIrEh9UZGxsjLS1N0i9LCCFECjRVufh4Yl+c+cADblaGbHlSXhne2nMb649FsdfH5IlEk1dNTQ1iYmKgoaFYo1YIIaSr69tDH8feGYHPpztBT6PhVOLRuxl4PeA6Lj3OkWF0TUkkeVVUVCAmJgarVq1CVlYWBg0aJImXJYQQ0olUuBz4DrfGpY/GYKJTD7Y8p6QKSw7exarfH6CgrFqGETYQe7ShOBPVMgwDXV1drFmzpkNBEUIIkR0zfU387OuOv6Of4b+nHiGPX5+wTkVmISwhD5unO2OyS492rY0oKWL3vBiGafFLRUUFFhYWmDNnDk6cOAE7OztpxkwIIaQTTHIxx8U1YzDLrWGJqPyyaqw8ch/vHr6P56WVMotN7J7XkydPpBkHIYQQOdRNRx3b3xiAqa7m+L+/HiG7pD5hnY/JRnhSPjZN64cZAyw7vRdGk/gRQgh5Jc++Zgj9cDR8hvRiy4orarDmjygsP3QPuaVVnRoPJS9CCCFi0ddUw9ZZ/RG0ZCh6dddiy0Mf5+D1gOs487BtE7h3hNjJKz09HefOnUNmZqZI+ZMnT+Dj4wM3Nzd4e3vj+vXrEg+SEEKI/BhpZ4zzq0Zj/jArtqywvAbvHXmAlUfud8qIRLGT1759+/DRRx+hrq5hyhA+n49FixbhwYMHqKioQEJCAt577z0kJiZKJVhCCCHyQUdDFV/McMHhxUNhYaDJlp99+KxT7gsTO3nduXMHtra2sLJqyLSnTp1CYWEhpkyZgosXL2LDhg2oqanBoUOHpBIsIYQQ+eJhb4zza0Zj3qCebFkevxpLD91FeGK+1I4rdvLKzc2FpaWlSFlYWBi4XC4++eQT9OrVCwsXLoS9vT3u3Lkj8UAJIYTIJ31NNXwzxxV7Fw6C6Yt1wRgGWHcsCvyqWqkcU+zkxefzoacnOk1+VFQUHBwcYGxszJbZ2toiJ0e+phEhhBAifZ59zRDyvgcMtNQAABmFFdhyNlYqxxI7eeno6OD58+fs86SkJBQUFMDNzU30BblcCAQCyUVICCFEYZjpa2LzdCf2eXBEGm7E50r8OGInLwcHBzx48ICdMf7PP/8Eh8PBkCFDROplZGTAxMREslESQghRGNNcLTDJuWFuxMO3UiV+DLGT1xtvvIGamhrMnDkTM2fOxP79+9G9e3e89tprbB0+n4/Y2FjY29tLPFBCCCGKgcPhwHd4b/a5NJZUETt5TZkyBStXrkRdXR1iY2NhYWGBwMBAkeVP/v77b9TU1DTpjRFCCOlauI2mi2IYyb++2HMbAsD777+P5cuXg8/no3v37k22jxw5EidPnkSvXr2a2ZsQQoiyYxgGp6Oy8PXf0p0Pt03JCwDU1dWbTVwAYGFhAQsLiw4HRQghRPHcTyvE52ce40FakUj5MJvmc0ZHtDl5Nefff//FkydPYGFhgfHjx0NFRUUSL0sIIUQBZBVV4OvzT3AqUnRuQyMddXz4Og8+g61a2LP9xE5eR48exf79+7F582aRlZI//fRTHD9+nH0+aNAg7N27F2pqapKNlBBCiFwpr67Fz9eT8MuNRFTWNNwipa7CxSIPa6x8zQ76mtLJBWInr9DQUOTl5cHV1ZUte/DgAY4dOwYdHR2MGzcODx48wN27d3HmzBnMnDlTKgETQgiRHYZhkFFYgRsJudhxOQE5JaJLoUx06oFPJvdFbyMdqcYhdvJKTEyEvb29SI/q7Nmz4HA42L59O8aMGYPCwkJ4enrir7/+ouRFCCFKoKK6Dg8zinA/rQj30wrxIK0Iefyma3c5Wejjv1P7YZiNUafEJXbyKiwsbDKbxt27d6Gvr48xY8YAALp164ZBgwYhPj5eslESQgiROmGvSpik7qUWIvZZCWoFLY91N9HTwLoJDpg9sCdUuJ23mrLYyUsgEKC6umGNFuESKKNHjxapZ2hoiMLCQslFSAghRCoqa+rwKKcSN54n4n5aIe6nFYm1IrKehioGWBlipJ0x5g/rDV0NiYz9axOxj9ijRw/ExjZMsPjPP/+grq4OAwcOFKlXWloKfX19yUVICCGkwwQCBplFFXiQXoT7qYV4kFaImCxhr6r1FZBtTXQw0Kob3Ht3w8De3WBnogtuJ/aymiN28vLw8EBwcDD8/f3h4eGBbdu2gcPhiEwPBYCdfYMQQkjnKq+uRXpBBVLzy5BWUI70gnKkFpQjraAcGQUVqK579aTpuhqqGNDLEAOtDOHWuxvcehnCUFu9E6JvG7GT1zvvvIPQ0FAEBwfj999/B8Mw8Pb2hp2dHVvn8ePHyMnJwcSJE6USLCGEdGUMwyC3tKo+IeXXJ6XGX+Kc8ntZT301DLM3w0CrbhjY2xD2pnqdeu2qvcROXqampjhx4gSOHj2K/Px89O/fH9OnTxepEx8fj3HjxuH111+XeKCEENIVVNbUIaPwRULKr+85pTdKUI3vp2qr7jrqcDTXq09UVt3gZmWIZ6mJcHR0lOA76BxtuspmYmKClStXtrh9xowZmDFjRkdjkkt///03QkJC8OjRI5SUlMDKygq+vr6YM2cOOBz5/y+FECJ7DMOgsLwGOSWVyCmpxPOSKjwrrkR6o2SVXVLZ7tdX5XLQs5sWenXXRm8jbVh1F37poFd3Leg1c8Pws468IRnq/CEiCmr//v2wtLTEhg0b0K1bN/z777/47LPP8OzZM3zwwQeyDo8QIkMMw6C4ogbPS6teJKaqF8npxePS+kSVW1ol1nWn1hhoqaG3kTZ6vUhMvV9879VdG+YGmlBVEXuxEIXW5uQVERGBw4cPIzIyEgUFBZg2bRq+/PJLAMDNmzdx584d+Pr6Kt2ClD/99JPIhMTDhw9HUVERDhw4gPfeew9cbtf4gSGkK2EYBqVVtQ1J6MX35y+SUU5JJXJK68uqayWzgrwKlwMLQ022xyTsPfU20kavbtow0Kap94A2Jq8ffvgBP/74I5hGi7M0fqyqqoo9e/bAzMwM//nPfyQXpRxobiZ9R0dHHD16FFVVVdDS0pJBVISQ9hAIGBSWVyOPX408fhXy+C/3mOoTVE5JFSpq6iR6bH1NVZjpa8JMXxOmehow1ddEz25a7Gk+C0MtqHWR3lNHiJ28rly5gl27dsHc3BwbNmzA4MGDMWLECJE6w4YNg6GhIa5evdopySspKQk3b95EdHQ0Hj16hJSUFDAMg8DAwFeOeAwJCUFwcDDi4uIgEAjQp08fzJ49Gz4+PmL3ou7duwdLS0tKXITIgZo6AQrLqpHLr6pPSqVVbGJqSFLVyC4qQ0lVMupamTWiPXQ1VGGqrwEzPU2Y6WvUJ6cXCao+WWnAVE8TWuq06oYkiJ28Dh06BHV1dfz666+wtbVttg6Hw4G1tTXS0tIkFmBrgoODcfDgwTbv5+/vjyNHjkBDQwPDhw+HqqoqwsPDsXnzZoSHhyMwMPCVy7rcvXsX586dw9q1a9sbPiHkFapq61pNRI3LC8slv9Q8AGipqaCHgSZMhElIT5iYNER6UDoymGWiKxO7tWNiYuDq6tpi4hLq0aMHnjyR7gqaQjweD4sXL4azszOcnZ2xceNGREREtLrPhQsXcOTIEZiYmODw4cOwtrYGAOTl5cHPzw8XL17E4cOHsWDBghZfIzs7G2vWrMHgwYOxcOFCCb4jQpQbwzAoq65D/ouEk1tajfyyKuSVNpy+y3+RnHL5VSitrJVaLAZaajDSVYexrgZMdDUaJaP63pPpi8e6Gqo0olgOiZ28KisrW1xBubGysrIOBdQWc+fObfM+u3fvBgCsXbuWTVwAYGxsjE2bNsHX1xd79uyBr69vs6cPS0pKsHTpUhgaGmLXrl208Cbp8gSC+pF2woSTz28uEVWzCasj9ym1hsMBumvXJyNjPXUY6Wiwj4UJylhXA4XZaRg2wAnqqnRdSZGJnbxMTEyQlJT0ynpPnz6FpaVlh4KSluzsbMTExEBNTa3Za2JDhgyBmZkZcnJyEBkZ2WTexsrKSixfvhylpaX4448/oKen11mhE9KpauoEKCirRm6paBISPhZeV8rnV6GgrLrVWcc7QoXLgZGOMCFpwFhXHSa6GmyPif3SU0d3bXWxhonHlqpS4lICYievoUOH4sSJEwgLC4OHh0ezdc6dO4esrCz4+flJLEBJevz4MQDA3t4empqazdZxcXFBTk4OYmNjRZJXbW0tVq9ejaSkJAQFBcHMzKxTYiZEUhiGQUllLXJL60fS5ZbW33eUy69CbsmL7y/KCsqrwUgnH0FDlfsi6TQkIDYZ6YmWG2qpyXwCWCKfxE5eixcvRkhICFatWoX169eLTAFVUVGBCxcu4IsvvoCWlhZ8fX2lEmxHZWRkAECrEwebm5uL1BXy9/fH1atXsWHDBvD5fERGRrLb7OzsoKur2+qxq6qqRGbll4bKykqpH0ORdJX2qKljUFRZh4LyWhRW1KGgsq7++4vnhS+2FVUmo7pOOhlJR40LQy0VGGqqsN+7NXpsqFX/3EBTBdpqnGauITEAKuu/qoG6AiCnAMiRQqxd5edCXIraHmInL1tbW3z11VfYsGEDNm3aBH9/f3A4HISEhODkyZMAABUVFXzzzTfo1auXtOLtkPLycgBodWi7jk790tUvX7v7559/AABfffVVk30OHjyIoUOHtnpsDQ0Nqc8fFhsbq5BzlEmLIrcHwzAoKq8R6Q293Ft6/qK3VCSFUXatXj/SabiOZKSrASMddWiqKc61X0X+uZAGeW6P1pJqm8Z2TpkyBXZ2dvjpp58QFhYGPp+P2tpaaGpqYvjw4Vi5ciWcnZ07HLC0CG+obs/IoStXrkg6HNJF1dQJkFNSiWfFL76KKvCsuBJZRRXIKalkk1ONFHpJWmoqMNXXgKmeBkz06gcxmOjV339koid8rIHuOuJdPyJEVsROXnw+HwDg4OCA77//vn6CycJCCAQCdOvWTSFG3Ql7VcIeWHOEPS5hXULaorZOgFx+FbKKKvGsuALZxZXs46ziSmQXV+B5aZVErydxOagfTdcoIZnqCxNTfVIqfZ6BYQP60b1IRGmI/ZM8aNAguLi44M8//wRQ33sRZ+i8PBGOgszKannV0OzsbJG6hAgJBAxy+VVsbylL2Gsqaeg9PS+tktjMDXoaqjDRqx/EwPaU2OSkyfaauuuov3L9pdjKHEpcRKmI/dOso6OD3r17SzMWqevXrx8AICEhAZWVlc2OOIyOjgYAuT0HTKRHuFyFcN2k9EYL/qUXliO7uFIiQ8I5L3pKFgaaMDfQgrmhJsyFjw002VN4NI0QIS1r04CNnBxpjP3pPObm5nByckJMTAzOnz/fZO2xiIgIZGdnw8TEBG5ubrIJkkhVVW0dMgsrGpIT+1WB9IJy8Ks6PqODsa46zA200MNAsz5BGWqJJCczfU26z4iQDhI7ec2dOxf/+9//8OjRI7kelPEqy5Ytw6pVq7Bt2za4ubmxvcn8/Hz4+/sDAJYuXUpLnCgohmGQX1aNtIJy3Eri41JWgkhP6llJZYeuN3XTVoO5gRYsDDUbEpShaGJSpJF3hCiqNiWvuLg4vP3221iyZAnGjx8PS0tLqKurSzO+VsXExLAJB6if3QMAAgICsHfvXrb86NGj7OOJEyfCx8cHwcHB8Pb2xogRI9iJefl8Pry8vDB//vzOexOkzeoEDLKKKpCSX4bkvDKk5peL9KTKqxsvYfG8Ta+to64CKyMdWHXXYtdR6vXiy9JQixITIXJC7OTV+BpQQEAAAgICWqzL4XDY2Sykic/nIyoqqkl5SkpKq/tt2rQJ7u7uCAoKQkREBAQCAWxsbNq8JAqRHoZh8Ly0Csl5ZUjJq09Swq/UgvJ2L/zH5QDmBg2JyarRirRW3bXRTVuNJmElRAGInbyYNpxraUvdjhg6dCji4uLata+3tze8vb0lHBFpq8KyaiQ1TlD5ZUjOLUNqfhnKqtu3CKCehiqsjLRhqFYH595mIsnJwlCLrjcRogTETl6dtcwJUT41dQKk5pchIYePp8/5SGrUiyquaN/sEMa6GuhjrI0+xjrobaTDrkJr1V0bBlr1vSd5njmAENIxdOMHkZjKmjok5tYnqKfP+fXJKpePlLyydg0x19dURR8TXfQx0kYfY11YG2vD5sV3PU01KbwDQoiioORF2qy0sqY+OT3nI/HF96fP+UgvLG/zSD4tNRVYG+vAxlgH1sb1SarPi+90/YkQ0pI2J6/q6mqEhoay90RxOByYmppiyJAhmDBhgkxHHxLJyudXsUnqaaOv7JLKNr+WpaEWbE11YW+qCxsTHfQx1oGNsS7M9DUoQRFC2qxNyev+/ftYu3Ytnj171mRQxrFjx7B9+3Z8++23GDRokESDJNJXUV2HiJQChCXkIiqjGE+f81FQVt2m1+BygN5GOrA10YW9mS7sXny3NdGlqYkIIRIl9l+UhIQELF68GBUVFejVqxemTJnCzv+XmZmJv//+G6mpqVi6dCmOHj0Ke3t7qQVNOk4gYPD4WQluJuQh7Gku7qQUij38XE2FAxtjXdiZNnzZm+nC2kiH7oMihHQKsZPXjh07UFFRgeXLl2PVqlVN7oX64IMPEBgYiN27d+OHH37Ajh07JB4s6ZisogqEJeTh5tM8/PM075U9Ky01Fdia6sDeVE8kUfXurk3LZRBCZErs5BUREYE+ffpgzZo1zW7ncrlYs2YNQkNDcfv2bYkFSNqPX1WLW4n5CHuah5sJuUjMLWu1vp2pLjzsjDHSzhh9e+jB0lCLlmAnhMglsZNXVVUVOyt7a/r160cLN8pInYDB/bRChCXkISwhD/fTClsdom6ko46RdsbwsDfGKHtjmBu0vMI0IYTIE7GTV58+fZCbm/vKerm5uQq/dIoiqaypw8kHmbgWl4uwhOfgVye3WFddlYsh1t0xyr4+YTn20KeeFSFEIYmdvN588034+/vj3r17cHd3b7bOvXv3cPfuXfz3v/+VWICkZfdSC7Duz4dIymv5dKCjuT5Gv0hWg62704AKQohSEDt5vfHGG0hKSsKSJUvw1ltvwdvbGz179gRQP9owJCQER44cgZ+fH3x8fKQWMKkf1v5daBx++ye5yU3BPfQ12dOAI2yNYaKnIZsgCSFEito1q/zevXtFlhxp7MCBAzhw4IBIWWfNMt8V3E0pwLpjD5HcqLelq6GKFa/ZwkajDBOG9aebfgkhSk8qs8pLcl/S4HFWCd7acxvVdQ33Y42yN8ZXs/vD0lALsbGxlLgIIV0CzSqvQH69mcQmLi01FWya1g/zBvWihEUI6XJozh4Fcju5gH380/yBGOtgKsNoCCFEdmiaBAWSX1bFPh5k3V2GkRBCiGxR8lIQFdV1qKypP2WopsKBjjoNeSeEdF2UvBREUUXDPIQaqiqoEnMSXUIIUUaUvBSEkY4GdF8sK8KvqoV/CN16QAjpuih5KQh1VS4+ntSXfR4ckYZTkZkyjIgQQmSHkpcCmT/UClP7m7PP1x97iN8j0ug+OkJIl0PJS4FwOBxsneWCPsY6AICqWgE2/BWND36PRGlljYyjI4SQzkPJS8Hoaaph38LB4JnpsmUhUVmY+kMYnuRWyjAyQgjpPJS8FJC1sQ5OrfSAz5BebFlqfjnWnMvC2j+j8LyUkhghRLlR8lJQWuoq2DqrP3b4uLGjEAHg2L0MvPbtNfx8PRFVtXUyjJAQQqSHkpeCm+ZqgXMfjIKXoxlbVlZdh6/+foLx22/gdFQWBK2spkwIIYqIkpcSsDLSxq8LBuELrx6wM224FpZWUI4Pgh/Ae2cYbsTn0qhEQojSoOSlRNwttfH3qlH4n3c/GGipseUxWSXw2xuB//x6G1HpRbILkBBCJISSl5JRU+Fi0cg+uLH+NawYawtNtYaP+N/EfEzf9Q9WBt1HUi5fhlESQkjHUPJSUgZaalg/sS+ur3sNPkOsoMJtWPPrbPQzjA+4gU9PRqOovLqVVyGEEPlEyUvJmelrYussF4SuGY3JLj3Y8joBg8O30uD53XUcvZNOgzoIIQqFklcXYWuiix//446TK0diuI0RW15QVo31xx9i7u5wZBSWyzBCQggRHyWvLmZAL0McWToUe/wGwdJQiy2/l1qIaTv/wb9P82QYHSGEiIeSl5hCQ0Ph4+ODoUOHwsXFBV5eXvj6669RWloq69DajMPhYHw/M1z6cAxWvmYL1RfXwwrKqjH/t9v4O/qZjCMkhJDWqb66CgGA4uJiDB48GIsWLYKBgQHi4+Oxc+dOxMXFYe/evbIOr1201FWwbkJfvOZgineD7iO3tAoCpn62emdLA/Tqri3rEAkhpFmUvMQ0d+5ckedDhw6FhoYG/vvf/yInJwdmZmYt7Cn/Bll3x5n3PTD353CkFZSjtKoW20LjEPimm6xDI4SQZtFpww4wNDQEANTW1so2EAkw1dPASDtj9nlMVokMoyGEkNYpdM8rKSkJN2/eRHR0NB49eoSUlBQwDIPAwEBMnDix1X1DQkIQHByMuLg4CAQC9OnTB7Nnz4aPjw+43JZzel1dHWpra5GQkIBdu3bhtddeg6WlpaTfWqcorqjB0+eliM/h43RkFsKT8tltPfQ1ZRgZIYS0TqGTV3BwMA4ePNjm/fz9/XHkyBFoaGhg+PDhUFVVRXh4ODZv3ozw8HAEBgZCRUWl2X2HDh3KDtIYNWoUtm/f3qH30BmKK2qQkFOKhOd8xOeU4umL7zklVc3W19VQxfqJDp0cJSGEiE+hkxePx8PixYvh7OwMZ2dnbNy4EREREa3uc+HCBRw5cgQmJiY4fPgwrK2tAQB5eXnw8/PDxYsXcfjwYSxYsKDZ/Q8dOoSKigokJCTgp59+wjvvvIN9+/a1mOw6S1VtHTKKq5Ed9xxp+eVIzitjk9Tz0uaT1MtUuBz4DOmFVeN4MNHTkHLEhBDSfgqdvF4eRCGO3bt3AwDWrl3LJi4AMDY2xqZNm+Dr64s9e/bA19e32dOHjo6OAICBAwfCyckJs2fPxsWLF195mlISisqrkVZQjtT8cqQVlCMtvxypBWVIyy/Hs5JK1E8anyH266mrcGFjogOemR54ZrqY7GIOGxPdV+9ICCEyptDJq62ys7MRExMDNTW1ZpPNkCFDYGZmhpycHERGRmLgwIGtvp6joyO4XC7S0tKkEm9xRQ12XE7A7eR8pOWXo6SyfQND1FW5sDXRBc9MF/amurA304O9qS6sumtDVYXG7BBCFE+XSl6PHz8GANjb20NTs/kBCS4uLsjJyUFsbOwrk9f9+/chEAjQs2dPiccKAL9HpOG3sGSx6nI4gIm2KmzNDNDbSBtWRtovEpYerLpri0zMSwghiq5LJa+MjPpTahYWFi3WMTc3F6krtHjxYgwbNgz29vZQV1fH48eP8dtvv8HBwQFeXl6vPHZVVRViY2PbFO9oU2D0Aps27dOgBkAhKnMLEZ/bzpdQAm1tc2VGbdGA2kKUvLZHVVXL1+u7VPIqL6+feFZLS6vFOjo6OgCAsrIykfL+/fvj9OnTbFLr2bMn3nrrLSxatAjq6uqvPPaAAQPaGTUhhJCXdankxdSPaACH0/ZTaKtWrcKqVaskHRIhhJB26FJX64W9KmEPrDnCHpewLiGEEPnTpZKXcCaMrKysFutkZ2eL1CWEECJ/ulTy6tevHwAgISEBlZWVzdaJjo4G0HA/FyGEEPnTpZKXubk5nJycUFNTg/PnzzfZHhERgezsbJiYmMDNjWZUJ4QQedWlBmwAwLJly7Bq1Sps27YNbm5u6N27NwAgPz8f/v7+AIClS5e2OjlvZzp06BDu3r2L+Ph4FBQUgM/nQ09PD3379sXMmTMxbdq0FgegtHfyYUWyfft2dtaU9evXY/Hixc3WU8a22LBhA06cONHi9j59+jT7TxqgnO0BAJWVlTh06BDOnz+P1NRU1NTUwMjICM7OzliwYAHc3d2b7KNsbXH79m34+fmJVffq1atNbh1SlPbgMMIheAooJiaGTTgA8PTpU5SVlcHa2hoGBgZs+dGjR0X227RpE4KDg6GhoYERI0awE/Py+Xx4eXlhx44dMp+rUGj06NEoKCiAvb09zMzMoKWlhaysLERFRYFhGIwbNw47d+5s8kPV0uTDZWVlGD9+fKuTDyuKhw8f4s0334RAIADDMC0mL2VtC2HyGjhwIPtPWGMmJib46KOPmpQra3ukp6dj8eLFSE1NhZGREVxdXaGuro7MzEw8efIEK1aswIoVK0T2Uca2SExMxJ49e1rc/vDhQyQmJsLKygqhoaEi//wqVHswCuzWrVsMj8d75VdzTp8+zbzxxhuMm5sb4+rqysycOZM5fPgwU1dX18nvonV37txhysrKmpTHx8czI0aMYHg8HnPs2DGRbefPn2d4PB4zcuRIJjk5mS3Pzc1lJk2axPB4PGb//v3SDl2qqqqqmMmTJzMeHh7MihUrGB6Px/z6669N6ilzW3z88ccMj8djjh8/LvY+ytoeZWVljJeXF8Pj8Zjvv/+eqa6uFtleUFDAJCUliZQpa1u8yuTJkxkej8f89NNPIuWK1h4Knby6up07dzI8Ho/58MMPRcpnzpzJ8Hg85sSJE032uX37NvsDKm+Jui2++eYbhsfjMZcvX2b/iDeXvJS5LdqTvJS1PbZt28bweDxm/fr1Yu+jrG3Rmvv37zM8Ho9xdHRksrOzRbYpWnvIzwlM0maqqvWXLBvP8CHu5MO5ubmIjIzsrFAlKioqCvv27cPUqVPh6enZYr2u0BZtoaztUV1dzV4aWLZsmVj7KGtbvMrx48cB1K9FaGZmxpYrYnt0uQEbyiI9PR2///47AIj8AZf05MPypqqqCh9//DEMDAywcePGVusqe1sI3b59G3FxcSgvL4eRkRHc3d0xcuTIJtdBlbU9YmJiUFRUBHNzc9ja2uL+/fu4du0aioqKYGxsjFGjRjUZPaysbdGaiooKnDt3DgAwZ84ckW2K2B6UvBTE8ePHcefOHdTU1CAnJwcPHjyAQCDA8uXLMX78eLZeRyYfVgQBAQFITk5GQEAAunfv3mpdZW8LoZMnTzYps7Ozw/bt2+Hg0LAitrK2R3x8PACgd+/ezY7A3LVrFyZMmIBvvvmG/cOsrG3RmvPnz6OsrAxGRkYYO3asyDZFbA9KXgri/v37Ir+UqqqqWLVqFRYtWiRSryOTD8u7+/fv48CBA/Dy8sLkyZNfWV+Z2wIA+vbti08//RTDhw+HhYUF+Hw+Hj9+jICAADx58gSLFi3CiRMn2NNDytoexcXFAIC7d++irq4Ob7/9Nnx8fGBoaIg7d+7A398fFy5cgI6ODrZu3QpAeduiNcJThtOnT4eamprINkVsD7rmpSC2bNmCuLg4REVF4ezZs/Dz88POnTsxb9485OTksPWYDkw+LM8qKyvxySefQFdXF//73//E2kdZ20Jo4cKF8PX1hZ2dHbS1tWFqaoqxY8fizz//xIABA5Cfn8/eAwcob3sIBAIAQG1tLebMmYOPP/4YVlZW0NfXx7hx47Br1y5wOBycPHkS6enpAJS3LVqSmpqKO3fuAGh6yhBQzPag5KVgNDU1YWdnh48//hgffvghnjx5gs8//5zdrqyTD2/fvh0pKSnYsGEDTE1NxdpHWdviVdTV1dmBC9evX2fLlbU9Gsc6b968JttdXFzg5OQEgUCA27dvi+yjbG3REmGvy83NDba2tk22K2J70GlDBTZr1ix8/fXXuHr1KmpqaqCmpqa0kw9funQJXC4XJ0+ebHKNJykpCQAQHByMa9euwcrKClu2bFHathCHjU39IqaNe+XK2h6NY21pVfOePXvi0aNHyMvLE9lH2dqiOXV1dezvzOzZs5uto4jtQclLgenr60NVVRW1tbUoLi6GsbFxk8mHmxs5pKiTDwsEAkRERLS4PT09Henp6SgpKQHQdCJmZWqLVykqKgIg+l+ysraHk5MT+7iwsLDZgTyFhYUAAG1tbQDK2xbNCQsLQ05ODrS1tVu8VqyI7UGnDRXYnTt3UFtbC319fXTr1g2A8k4+fOXKFcTFxTX7NXPmTAD1cxvGxcXh1KlTAJS3LcTx999/AwCcnZ3ZMmVtDzMzM7i6ugIAbt261WR7cXExOxRc2B7K2hbNOXbsGABg0qRJLZ7yU8T2oOQlx+7evYvTp0+jurq6ybZ79+6x9znNmTNHZL4x4fWObdu2ITU1lS2X18mHpUlZ2yI2NhZXr15FXV2dSHltbS327duHQ4cOAagf1NGYsrbHO++8A6B+WHxsbCxbXlVVhU2bNqG0tBROTk4if3iVtS0aKygowNWrVwE0P1CjMUVrDzptKMfS0tLwySef4PPPP0e/fv1gbGyMsrIypKen4+nTpwCAsWPHYtWqVSL7TZw4ET4+PggODoa3t3ezkw/Pnz9fFm+p0ylrW2RmZmLlypUwNDSEtbU1zMzMUFZWhvj4eDx//hxcLhdr167FqFGjRPZT1vbw9PTE22+/jb1792Lu3LlwdXWFoaEhHj58iOfPn8PMzAzbt28XGU2nrG3R2OnTp1FTUwMbG5tX3lisaO2h0LPKK7v09HT89ddfuHv3LtLS0lBYWAiGYWBiYgJnZ2dMmzYNXl5eLe4fEhKCoKAgxMfHQyAQwMbGRi6XNugo4Y2pr1oSRZnaIj09HQcPHkR0dDQyMzNRVFQEDoeDHj16wN3dHf/5z39EThm+TNnaQ+jixYs4dOgQYmNjUVFRAQsLC3h6emLZsmUt3tSurG0BAN7e3oiPj8e6deuwZMkSsfZRlPag5EUIIUThyE8aJYQQQsREyYsQQojCoeRFCCFE4VDyIoQQonAoeRFCCFE4lLwIIYQoHEpehBBCFA4lL8Ly9PSEg4PDK7+Ey0rIgjAGSbh9+zYcHBzg6+srkdeTlg0bNsDBwQF//fWXrENpN+HPlryswisrvr6+Ir9LU6dObVInOTkZS5cuxYABAzBo0CB89NFHyM/Pb/E1V6xYgcGDB7Mz5jdn+vTpIseV9595cdD0UKQJDw8PmJiYtLjd2Ni4E6MhRPkIf8eEq1wLVVRUYMGCBcjJycHIkSNRXl6OM2fO4OnTpzh27FiTFZAvXryIy5cv4/PPP2/199LT0xOOjo7Izc1FWFiYVN5TZ6PkRZpYtmwZhg4dKuswpK5///44d+5cq0ufE8nYv38/ampqmvyx7qpa+h37448/kJOTg9WrV+Pdd98FAHzyySf466+/cOnSJUyaNImty+fz8cUXX8Dd3R1z585t9XjC+U9v376tNMmLThuSLktLSwu2trawsLCQdShKz8rKCra2tk16DkRUTEwMANFFI4WJKTIyUqTu999/j/z8fGzevFlkwuGugpIXaTeGYbBkyRI4ODjgv//9b5PtAoEACxYsgIODAzZv3syWZ2RkwMHBAZ6enqitrcUvv/yCSZMmwcXFBSNGjMDHH3/c6oquzXn69CkCAwPx5ptvwsPDA87Ozhg2bBiWLl2KGzduNLtPS9e8GsfHMAyCgoIwffp0uLq6YvDgwXj33XcRHx/fYiyFhYUICAiAt7c33NzcMGDAAMycOZPtfTSnvLwcAQEB8PLygrOzM8aMGYNNmzaxiyi2xbx58+Dg4IBLly61WOfrr7+Gg4MDvv76a7asoKAABw4cwOLFi+Hp6QkXFxe4u7tj3rx5CAoKarL8CtD0s/ztt98wbdo09nqNUEvXvDIzM7F79274+vpizJgxcHZ2xpAhQ+Dr64uQkJBmY2/8udXU1OCnn37CxIkT4eLiguHDh2Pt2rWt/vw8e/YMW7duxeTJkzFgwAAMHDgQkyZNwqZNm5r9XNvzebaXcBFRAwMDtkz4uKqqii2Ljo5GUFAQli5dCjs7O4nGoCgoeZF243A4+Oabb2BmZoajR4/izJkzItt37dqFW7duoV+/ftiwYUOzr7FmzRrs2LEDFhYW8PLygpqaGk6ePIk5c+YgKSlJ7Fj27duHH3/8ESUlJejbty+8vLxgaWmJGzduYOnSpdi3b1+73uOGDRvw1VdfwcjICGPGjIGenh6uXLkCHx8fpKenN6kfFxeHadOm4eeff0ZJSQmGDBmCwYMHIysrC1u3bsXSpUubrM9WXl4OPz8//PzzzygsLMTo0aPZU5rz5s1DcXFxm2IWLs554sSJZrfX1dWxiWHWrFls+c2bN/Hll1/i6dOnsLS0hJeXF/r164fY2Fhs3rwZ77//Plqax5thGLz//vsICAiAkZERPD09xfqjeurUKWzfvh3Z2dmwsbHB+PHjYWtri3v37mHt2rX44osvWty3pqYGS5cuxS+//ILevXtj9OjR4HK5CAkJwVtvvcWuqN1YWFgYpk6div3796O0tBQeHh7w8PCApqYm/vjjD1y4cEGkfns+z46wtLQEAJGf/eTkZABAz549AdR/fp999hmsrKzYdcy6JIaQF1577TWGx+Mxt27datN+d+7cYRwdHRk3NzcmOTmZYRiGCQ8PZ/r27cu4ubkxKSkpIvXT09MZHo/H8Hg8Zvjw4UxCQgK7raqqilm7di3D4/GY2bNnNzmWcL+X3b59m0lPT29SHhkZyQwcOJBxcnJinj17JrLt1q1bDI/HY+bPn99ifJ6enkxqaqpIfEuXLmV4PB6zceNGkf0qKioYT09PhsfjMT///DNTU1PDbissLGQWLlzI8Hg8ZseOHSL7bd26leHxeMzUqVOZvLw8try4uJh544032FiOHz/e5P01p7i4mHFxcWGcnJyY/Pz8JtuvXbvG8Hg8ZtasWSLlT58+ZSIjI5vUz8nJYaZPn87weDzm7NmzItsat9XYsWObfNZCwp+tlz+jqKgoJj4+vkn95ORkZsyYMQyPx2sSk/BzE76Hxm1WUlLCzJw5k+HxeMyPP/4osl9mZibj5ubG8Hg8ZteuXSKfj3B7dHQ0+7y9n2dr5s+f3+rvmPCzWbBgAZOTk8OkpKQw3t7ejKOjI/t7snfvXobH4zHh4eFiH1eopZ95RUQ9L9KEn59fi8PkG58KEho0aBA++OADlJWVYfXq1cjKysLatWshEAjwxRdfoHfv3i0ea8WKFSL/oaurq+Ozzz6Dnp4eoqOjce/ePbFiHjJkCPufaWOurq6YP38+ampqcPnyZbFeq7FPP/0UVlZWIvG99957AIDw8HCRun/99RcyMjIwadIkLF++HKqqDeOhDA0N8dVXX0FNTQ1BQUFsD6ayshJ//PEHeywjIyN2H319fWzatKnN1zP09fXh5eWFmpqaJr1hYZyAaK8LAGxtbeHq6tqkvqmpKdatWwcAzS4RL/Thhx+2+lk3p3///rC3t29Sbm1tjRUrVrR6TA6Hgy+//FKkzfT09Nh1q17+fPbt24eysjJMnjwZK1asEPl8AMDCwkJkDbT2fJ4dNWbMGEyYMAHh4eEYNWoUXn/9dcTFxeHtt9+GnZ0dsrKysGPHDsyaNQvDhg1j96upqZH4KUx5R6MNSROtDZXX1NRstnz58uW4c+cOwsLCMG3aNJSWluKNN97A5MmTWz3WtGnTmpTp6elh7NixCAkJQUREBNzd3cWKm8/n4/r164iNjUVxcTH7y5ySkgKg4fSLuFRVVZusRAwANjY2AIDnz5+LlAuvrU2cOLHZ1zMzM0Pv3r3x9OlTpKSkoE+fPnj06BHKy8thZmbW7Oizvn37wsHBAU+ePGlT7DNnzsTZs2dx4sQJ+Pn5seUlJSW4cuUK1NTUMGXKlCb71dbW4tatW4iMjERubi6qq6vBMAzKysoANLRlc8aPH9+mGIWqqqpw8+ZNREdHo7CwkD0Nl5ub2+oxLSwsmr3nr6XP5+bNmwDwypF5Qu35PCXh+++/x/nz5xEREQFVVVWMGTOG/TncvHkzNDQ0sH79egD11762bNnCDuZwdXXFp59+ChcXF4nEIs8oeZEm2jNUnsPh4Ntvv8W4ceNQWloKOzs7bNy4sdV99PX1oa+v3+w24bn/7OxssY5/6dIlbNy4kb3g3RzhH2BxmZiYNPnvHAB0dXUBoMm1DuE1MOGw5NYUFBSgT58+yMnJAYBme41ClpaWbU5eI0eOhJmZGR4/foy4uDj2j/zZs2dRXV2NCRMmwNDQUGSf5ORkrFy5EomJiS2+Lp/Pb7bcyMioxX9sWvPgwQOsXr261c+5pWOam5s3W97S5yMcxCFMbq/Sns9TErhcLiZPntzkH78LFy7g6tWr+Prrr9GtWzdkZmZi4cKF0NfXx1dffQUul4uAgAAsXLgQISEhSj+KlpIXkZhLly6hvLwcQH3SycnJETnl1h7inDLLzs7GRx99hMrKSixbtgxTp06FpaUltLW1weVy8ccff+Czzz5r86mdti55LhyNN3bsWHTr1q3Vui8nDknjcrmYMWMGdu/ejRMnTrADZk6ePAmg6SlDAPjggw+QmJgIT09PLFmyBLa2ttDT04OKigqSk5Nb7IEALffIW1NRUYH33nsPeXl5mDNnDnx8fNC7d2/o6OiAy+UiLCwMixcvbvU9tkVbT7/K0+fJ5/OxZcsWDB8+HDNmzAAABAcHg8/n44cffsCIESMA1E8gsGjRIgQHB+Ojjz6SakyyRsmLSER8fDy2bNkCNTU1TJw4ESEhIVi9ejV+//13qKurN7tPSUkJSktLoaen12RbZmYmgPrrLa9y9epVVFZWYsKECc3+wqamprbx3bSPubk5kpOT4ePjg7Fjx4q1j/CmXeH7bU5r21ozc+ZM7N69GyEhIVi7di3S0tIQGRkJExOTJqdDExMTER8fDyMjI+zcuRMqKioi29PS0toVQ2vu3LmDvLw8ODk5YcuWLU22S/pzE34+ycnJ6NGjh9j12/J5SktAQACKiorg7+/Plgl74wMGDGDLhI/b2lNXRDRgg3RYeXk5Vq9ejcrKSqxduxbffPMNhg4dipiYGHzzzTet7nv69OkmZaWlpbh27RqA+oEYryIcSt7cH6Tq6mqEhoaK8S46bvTo0QBaH9TwMicnJ2hrayM7Oxt37txpsj0+Ph5xcXHtiqdPnz5wc3NDXl4ebt68yfa6vL29myQnYRuampo22QY0/zl1lPCYLZ3+a26wSUd4eHgAAP7880+x6rfn85SGhw8f4siRI1ixYoXIgBhtbW0A9YN+hCoqKgC0vZepiCh5kQ7bvHkze7pp4cKF4HK52LZtG4yMjHDo0KFWb5b98ccfRa6x1NTUYMuWLSgtLYWTk1OzoxtfJryGERoaKjI5aXV1NT7//PNm78eShnnz5sHc3BwnTpzAjh072D8kjcXHx+P48ePscy0tLcyZMwcAsGXLFhQUFLDbSktL4e/v36GRbMLTg8ePH2cTUHOnDK2trcHlcpGQkNAkiR4/fhxnz55tdwwtEX5ut27dEvkZEAgE2LlzJ+7fvy/R4y1atAja2to4e/Ysdu/e3eSm62fPnuHRo0fs8/Z8npImvKfL1ta2ySlUHo8HACLHFz52dHSUWkzygk4bkiZ++eWXFm9wBYCpU6ey/8WePHkSJ06cgLm5ObZu3crWMTU1xTfffIMlS5bg//7v/+Do6MgOwhCysLCAk5MTpk+fjmHDhkFPTw+RkZHIyspCt27dXtlrE/L09ES/fv3w+PFjvP766xgyZAg0NDRw//598Pl8+Pr64tChQ+1oibbR0dHB7t278c4772DXrl0ICgoCj8eDiYkJ8vLykJGRgczMTLi6uopM/7N69Wrcu3cPMTExGD9+PIYNGwYVFRXcvn0benp68PT0xJUrV9oV0+TJk7FlyxZcvHgRAODs7Nzs0PTu3bvjrbfewuHDh+Hn54fBgwfDxMQE8fHxiI+Px/Lly7F79+72NUwLnJyc8Nprr+Hq1auYMWMGhg4dyt4i8ezZMyxZsgS//vqrxI5naWmJ77//HqtXr8b27dsRFBQEV1dXcDgcZGRkIDY2FitWrGCHy7f385Sk/fv348mTJzhy5EiTqbXmz5+PQ4cOYdu2bfj3338B1N8eYGhoiLfeeksq8cgTSl6kiVdN3Nm3b194eHggMTER/v7+UFVVxXfffdfkorWHhwc7A8KHH36Iw4cPi/wCcjgcfP/99/jll19w6tQpZGVlQVdXF9OmTcOqVataHYHXmKqqKg4dOoSffvoJly5dwj///AMDAwMMGTIE7733XpM54aTJwcEBp0+fxpEjR3D58mU8fvwYVVVV6N69O8zNzTF9+nRMmDBBZB8dHR0cOnQIP//8M86ePYvr16+je/fumDBhAlavXi12Em+Orq4uxo8f3+yMGi/buHEjHBwcEBwcjOjoaKiqqsLJyQl79uyBjY2NxJMXAOzYsQP79+/H6dOnERERAW1tbQwYMADfffcdKisrJZq8gPr7qE6fPo19+/YhLCwM165dg7q6OszMzODj4yMy8S3Qvs9TUjIzM7Fz507MmzcPAwcObLLd0NAQBw4cwLfffov79++DYRiMHDkSGzZs6BITIHMYSd1dR4iYMjIyMG7cOFhaWra7R0GIIvL19UVERAQOHjwok5Ubbt++DT8/PwwZMqRTzkZIE/W8CCGkkwlPzZuZmWHNmjVSP15gYCCePXvG3vitDCh5EUJIJxOemre3t++U5HXlyhWlGz5PyYsQQjqJrE7VnTp1SibHlSa65kUIIUTh0H1ehBBCFA4lL0IIIQqHkhchhBCFQ8mLEEKIwqHkRQghROFQ8iKEEKJw/h/2ktPzY5YYGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.array(E_arr)*100, ds_sel.plev.values, lw = 3)\n",
    "plt.yscale('log')\n",
    "plt.ylabel('pressure [hPa]')\n",
    "plt.xlabel('Explained variance [%]')\n",
    "plt.ylim(1000,0.1)\n",
    "plt.savefig('NAM_explained_variance.pdf', bbox_inches='tight')"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}