13_tmin.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:32:00.782735Z",
     "start_time": "2020-05-16T15:31:55.302410Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "variable = 'tasmin'\n",
    "var_name = variable + \"/\"\n",
    "model_name_file = 'NorESM2-MM'\n",
    "import numpy as np\n",
    "import time\n",
    "import xarray as xr\n",
    "import os\n",
    "os.environ['KMP_DUPLICATE_LIB_OK']='True'\n",
    "import pickle as pkl\n",
    "import pickle\n",
    "import numpy as np\n",
    "from numpy import *\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from datetime import datetime\n",
    "import sys\n",
    "from rpy2.robjects.packages import importr\n",
    "import rpy2.robjects.numpy2ri\n",
    "rpy2.robjects.numpy2ri.activate()\n",
    "qmap = importr('qmap')\n",
    "from rpy2.robjects import pandas2ri\n",
    "from rpy2.robjects import r\n",
    "utils = importr(\"utils\")\n",
    "#extremes = importr(\"extRemes\")\n",
    "#ismev = importr(\"ismev\")\n",
    "#eva = importr(\"eva\")\n",
    "gnFit = importr(\"gnFit\")\n",
    "from datetime import datetime\n",
    "statis = importr('stats')\n",
    "import time\n",
    "from datetime import datetime\n",
    "from matplotlib.offsetbox import AnchoredText\n",
    "from netCDF4 import num2date, date2num\n",
    "statis = importr('stats')\n",
    "import rpy2.robjects.numpy2ri\n",
    "rpy2.robjects.numpy2ri.activate()\n",
    "import glob\n",
    "import time\n",
    "from tqdm import tqdm\n",
    "import sys\n",
    "sys.setrecursionlimit(10000)\n",
    "import sys\n",
    "import time\n",
    "import traceback\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import glob\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:32:00.789728Z",
     "start_time": "2020-05-16T15:32:00.785977Z"
    }
   },
   "outputs": [],
   "source": [
    "if sys.platform ==\"darwin\":\n",
    "\n",
    "    gcm_path = '/media/udit/cmip6/cmip6_raw_amar_regrided/'\n",
    "    obs_path = '/Volumes/cmip6/'\n",
    "    target_folder = '/Volumes/cmip6/'\n",
    "else:\n",
    "    gcm_path = '/media/udit/cmip6/cmip6_raw_amar_regrided/'\n",
    "    obs_path = '/media/udit/cmip6/'\n",
    "    target_folder = '/media/udit/cmip6/'\n",
    "    \n",
    "model_name = model_name_file+'/'\n",
    "\n",
    "scenarios = ['hist','ssp126','ssp245','ssp370','ssp585']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:32:00.804189Z",
     "start_time": "2020-05-16T15:32:00.791508Z"
    }
   },
   "outputs": [],
   "source": [
    "def bias_correction(x,y):\n",
    "    q_map = qmap.fitQmap(x, y,method=\"RQUANT\", qstep=0.01, wett_day=False)\n",
    "    qm1 = qmap.doQmap(y, q_map)\n",
    "    bias_corrected_output = {}\n",
    "    bias_corrected_output['params'] = q_map\n",
    "    bias_corrected_output['outputs'] = qm1\n",
    "    return bias_corrected_output\n",
    "\n",
    "def bias_correction_model(y,q_map):\n",
    "    qm1 = qmap.doQmap(y, q_map)\n",
    "    bias_corrected_output = {}\n",
    "    bias_corrected_output['outputs'] = qm1\n",
    "    return bias_corrected_output\n",
    "\n",
    "\n",
    "\n",
    "def save_as_pickled_object(obj, filepath):\n",
    "    \"\"\"\n",
    "    This is a defensive way to write pickle.write, allowing for very large files on all platforms\n",
    "    \"\"\"\n",
    "    max_bytes = 2**31 - 1\n",
    "    bytes_out = pickle.dumps(obj)\n",
    "    n_bytes = sys.getsizeof(bytes_out)\n",
    "    with open(filepath, 'wb') as f_out:\n",
    "        for idx in range(0, n_bytes, max_bytes):\n",
    "            f_out.write(bytes_out[idx:idx + max_bytes])\n",
    "\n",
    "\n",
    "def try_to_load_as_pickled_object_or_None(filepath):\n",
    "\n",
    "    max_bytes = 2**31 - 1\n",
    "    try:\n",
    "        input_size = os.path.getsize(filepath)\n",
    "        bytes_in = bytearray(0)\n",
    "        with open(filepath, 'rb') as f_in:\n",
    "            for _ in range(0, input_size, max_bytes):\n",
    "                bytes_in += f_in.read(max_bytes)\n",
    "        obj = pickle.loads(bytes_in)\n",
    "    except:\n",
    "        return None\n",
    "    return obj\n",
    "\n",
    "def read_netcdfs(files, dim, transform_func=None):\n",
    "    def process_one_path(path):\n",
    "        # use a context manager, to ensure the file gets closed after use\n",
    "        with xr.open_dataset(path) as ds:\n",
    "            # transform_func should do some sort of selection or\n",
    "            # aggregation\n",
    "            if transform_func is not None:\n",
    "                ds = transform_func(ds)\n",
    "            # load all data from the transformed dataset, to ensure we can\n",
    "            # use it after closing each original file\n",
    "            ds.load()\n",
    "            return ds\n",
    "\n",
    "    paths = sorted(glob.glob(files))\n",
    "    datasets = [process_one_path(p) for p in paths]\n",
    "    combined = xr.concat(datasets, dim)\n",
    "    return combined\n",
    "\n",
    "def preproc_hist(ds):\n",
    "    ds = ds.sel(time = slice('01-01-1951','2014-12-13'),lat = slice(5,38),lon = slice(60,100))\n",
    "    return ds\n",
    "def preproc_model(ds):\n",
    "    ds = ds.sel(lat = slice(5,38),lon = slice(60,100))\n",
    "    return ds\n",
    "\n",
    "\n",
    "\n",
    "def preproc_hist2(ds):\n",
    "    start = pd.to_datetime('1951-01-01', format='%Y%m%d', errors='ignore')\n",
    "    end = pd.to_datetime('2014-12-13', format='%Y%m%d', errors='ignore')\n",
    "    ds = ds.sel(time = slice(start,end),lat = slice(5,38),lon = slice(60,100))\n",
    "    return ds\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# here we suppose we only care about the combined mean of each file;\n",
    "# you might also use indexing operations like .sel to subset datasets\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:32:00.809716Z",
     "start_time": "2020-05-16T15:32:00.805775Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-MM/\n"
     ]
    }
   ],
   "source": [
    "model_path = gcm_path + model_name\n",
    "model_path_historical = model_path\n",
    "model_path_RCP126 = model_path \n",
    "model_path_RCP245 = model_path\n",
    "model_path_RCP370 = model_path\n",
    "model_path_RCP585 = model_path\n",
    "\n",
    "print(model_path_historical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:32:00.815013Z",
     "start_time": "2020-05-16T15:32:00.811245Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-MM/tasmin/tasmin_day_NorESM2-MM_historical_r1i1p1f1_gn_20150101-21001231.nc']\n"
     ]
    }
   ],
   "source": [
    "prcp_hist = model_path +var_name+ \"*historical*.nc\"\n",
    "file_lists = glob.glob(prcp_hist)\n",
    "print(file_lists)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:32:00.826099Z",
     "start_time": "2020-05-16T15:32:00.816859Z"
    }
   },
   "outputs": [],
   "source": [
    "if variable == \"pr\":\n",
    "    observed = (obs_path +\"final_kelvin.nc\")\n",
    "    prcp_hist = model_path +var_name+ \"*historical*.nc\"\n",
    "    file_lists = glob.glob(prcp_hist)\n",
    "    try:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist)\n",
    "    except:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist2)\n",
    "    \n",
    "    observed = xr.open_dataset(observed)\n",
    "    observed = observed['pr']\n",
    "    observed.isel(time = 1).plot()\n",
    "    model_hist.pr.values = model_hist.pr.values*86400\n",
    "    model_hist = model_hist['pr']\n",
    "    lats = observed.lat.values\n",
    "    lons = observed.lon.values\n",
    "    prcp_245 = model_path +var_name+ \"*245*.nc\"\n",
    "    prcp_126 = model_path +var_name+ \"*126*.nc\"\n",
    "    prcp_585 = model_path +var_name+ \"*585*.nc\"\n",
    "    prcp_370 = model_path +var_name+ \"*370*.nc\"\n",
    "#     print(\"reading 45\")\n",
    "    print(prcp_245)\n",
    "    print(prcp_126)\n",
    "    file_lists = glob.glob(prcp_245)\n",
    "    print(file_lists)\n",
    "    model_245 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "\n",
    "    file_lists = glob.glob(prcp_126)\n",
    "    model_126 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(prcp_370)\n",
    "    model_370 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(prcp_585)\n",
    "    model_585 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    \n",
    "    \n",
    "    model_245.pr.values = model_245.pr.values*86400\n",
    "    model_585.pr.values = model_585.pr.values*86400\n",
    "    model_126.pr.values = model_126.pr.values*86400\n",
    "    model_370.pr.values = model_370.pr.values*86400\n",
    "    \n",
    "    \n",
    "    model_245 = model_245['pr']\n",
    "    model_585 = model_585['pr']\n",
    "    model_126 = model_126['pr']\n",
    "    model_370 = model_370['pr']\n",
    "    model_245.close()\n",
    "    model_585.close()\n",
    "    model_126.close()\n",
    "    model_370.close()\n",
    "    model_hist.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:32:12.838405Z",
     "start_time": "2020-05-16T15:32:00.827472Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-MM/tasmin/*historical*.nc\n",
      "['/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-MM/tasmin/tasmin_day_NorESM2-MM_historical_r1i1p1f1_gn_20150101-21001231.nc']\n",
      "hello\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:14: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "  \n",
      "/home/udit/anaconda3/envs/master/lib/python3.6/site-packages/xarray/backends/api.py:931: FutureWarning: Also `open_mfdataset` will no longer accept a `concat_dim` argument.\n",
      "To get equivalent behaviour from now on please use the new\n",
      "`combine_nested` function instead (or the `combine='nested'` option to\n",
      "`open_mfdataset`).The datasets supplied have global dimension coordinates. You may want\n",
      "to use the new `combine_by_coords` function (or the\n",
      "`combine='by_coords'` option to `open_mfdataset`) to order the datasets\n",
      "before concatenation. Alternatively, to continue concatenating based\n",
      "on the order the datasets are supplied in future, please use the new\n",
      "`combine_nested` function (or the `combine='nested'` option to\n",
      "open_mfdataset).\n",
      "  from_openmfds=True,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-MM/tasmin/tasmin_day_NorESM2-MM_ssp245_r1i1p1f1_gn_20150101-21001231.nc']\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:27: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "/home/udit/anaconda3/envs/master/lib/python3.6/site-packages/xarray/backends/api.py:931: FutureWarning: Also `open_mfdataset` will no longer accept a `concat_dim` argument.\n",
      "To get equivalent behaviour from now on please use the new\n",
      "`combine_nested` function instead (or the `combine='nested'` option to\n",
      "`open_mfdataset`).The datasets supplied have global dimension coordinates. You may want\n",
      "to use the new `combine_by_coords` function (or the\n",
      "`combine='by_coords'` option to `open_mfdataset`) to order the datasets\n",
      "before concatenation. Alternatively, to continue concatenating based\n",
      "on the order the datasets are supplied in future, please use the new\n",
      "`combine_nested` function (or the `combine='nested'` option to\n",
      "open_mfdataset).\n",
      "  from_openmfds=True,\n",
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:30: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:33: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:36: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n"
     ]
    }
   ],
   "source": [
    "if variable == \"tasmax\" or variable == \"tasmin\":\n",
    "    observed = (obs_path +\"final_kelvin.nc\")\n",
    "    observed = xr.open_dataset(observed)\n",
    "    observed.close()\n",
    "    observed = observed[variable]\n",
    "    temp_hist = model_path +var_name+ \"*historical*.nc\"\n",
    "    print(temp_hist)\n",
    "    file_lists = glob.glob(temp_hist)\n",
    "    print(file_lists)\n",
    "    print(\"hello\")\n",
    "    try:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist)\n",
    "    except:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist2)\n",
    "    \n",
    "    model_hist[variable].values = model_hist[variable].values\n",
    "    model_hist = model_hist[variable]\n",
    "    lats = observed.lat.values\n",
    "    lons = observed.lon.values\n",
    "    temp_245 = model_path +var_name+ \"*245*.nc\"\n",
    "    temp_126 = model_path +var_name+ \"*126*.nc\"\n",
    "    temp_585 = model_path +var_name+ \"*585*.nc\"\n",
    "    temp_370 = model_path +var_name+ \"*370*.nc\"\n",
    "#     print(\"reading 45\")\n",
    "    file_lists = glob.glob(temp_245)\n",
    "    print(file_lists)\n",
    "    model_245 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "\n",
    "    file_lists = glob.glob(temp_126)\n",
    "    model_126 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(temp_370)\n",
    "    model_370 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(temp_585)\n",
    "    model_585 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    \n",
    "    \n",
    "    model_245[variable].values = model_245[variable].values\n",
    "    model_585[variable].values = model_585[variable].values\n",
    "    model_126[variable].values = model_126[variable].values\n",
    "    model_370[variable].values = model_370[variable].values\n",
    "    model_245 = model_245[variable]\n",
    "    model_585 = model_585[variable]\n",
    "    model_126 = model_126[variable]\n",
    "    model_370 = model_370[variable]\n",
    "    model_245.close()\n",
    "    model_585.close()\n",
    "    model_126.close()\n",
    "    model_370.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:33:36.468442Z",
     "start_time": "2020-05-16T15:32:12.841468Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading hist\n",
      "reading 45\n",
      "reading 26\n",
      "reading 85\n",
      "reading 70\n"
     ]
    }
   ],
   "source": [
    "print(\"reading hist\")\n",
    "model_hist = model_hist.interp(lat = lats, lon = lons)\n",
    "print(\"reading 45\")\n",
    "\n",
    "model_245 = model_245.interp(lat = lats, lon = lons)\n",
    "print(\"reading 26\")\n",
    "\n",
    "model_126 = model_126.interp(lat = lats, lon = lons)\n",
    "print(\"reading 85\")\n",
    "\n",
    "model_585 = model_585.interp(lat = lats, lon = lons)\n",
    "print(\"reading 70\")\n",
    "model_370 = model_370.interp(lat = lats, lon = lons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:33:38.735589Z",
     "start_time": "2020-05-16T15:33:36.470823Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9e3c016128>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9d5icV3mwfz/Ty/aq1WrVrC4XyZY7BmyMARNaDAn+gmOCiQOBjxpCdX4hCQnlC/wIXyjGJpSYamM7BlOMsbGNq7olS7Z62V1tb9Pb+f54y5Sd2Z3tq91zX9deeufMed95RzNznvN0UUqh0Wg0msWHY65vQKPRaDRzgxYAGo1Gs0jRAkCj0WgWKVoAaDQazSJFCwCNRqNZpLjm+gYmQkNDg1q5cuVc34ZGo9GcVezYsaNXKdVYOH5WCYCVK1eyffv2ub4NjUajOasQkRPFxsc1AYmIT0SeFZE9IrJfRD5rjq8SkWdE5JCI/EREPOb4O0WkR0R2m3/vLnHdi0TkeRE5LCL/ISIylTeo0Wg0molRjg8gDlyjlLoA2AK8VkQuA74AfEUptRYYAG7JOecnSqkt5t8dJa77DeBWYK3599rJvgmNRqPRTJxxBYAyCJkP3eafAq4B7jbHvwe8udwXFZEWoEop9ZQyUpG/P5HzNRqNRjN1yooCEhGniOwGuoGHgCPAoFIqZU45DbTmnHKDiOwVkbtFpK3IJVvNcywKz8997VtFZLuIbO/p6SnndjUajUZTBmUJAKVUWim1BVgGXAJsLDbN/PcBYKVS6nzgdxjaQSHF7P1FixIppW5XSm1TSm1rbBzlxNZoNBrNJJlQHoBSahB4FLgMqBERK4poGdBhzulTSsXN8W8DFxW51GnzHAv7fI1Go9HMDuVEATWKSI157AeuBQ4AjwBvNafdDNxvzmnJOf2N5tw8lFKdwIiIXGZG//yldb5Go9FoZodyNIAW4BER2Qs8BzyklPoF8HHgIyJyGKgH7jTnf8AMF90DfAB4p3Uh049g8V7gDuAwhk/hV1N8LxrNmPx2/xk6BqNzfRsazbxBzqZ+ANu2bVM6EUwzGVLpDOtv+zX/65Ll/PObz53r29FoZhUR2aGU2lY4rmsBaRYF/ZEE6YxiX8fQXN+KRjNv0AJAsyjoCyUAONA5TCqdmeO70WjmB1oAaBYFvSEjMC2WzHC0NzzHd6PRzA+0ANAsCiwBALCvXZuBNBrQAkCzSLBMQG6nsK99eI7vRqOZH5xV5aA1mokyHEtypDtETyiOx+ng3NYq7QjWaEy0BqBZ0Hzq58/zlq8/yQsdwzRUeDi3tZoXOobJZM6e8GeNZqbQAkCzoDncbRSyffxQL/UVXs5dWk0onuJEf2SO70yjmXu0ANAsaFzObN3BpkovDZUeAIaiybm6JY1m3qB9AJoFi1KK470R2ur8XL2+iRsvWW5HA+lcAI1GCwDNAqY3lCAUT/HR69bxV1euAuCpI30AJNPaB6DRaBOQZsFyvM9I+FrZELTHPC7DJJTUGoBGowWAZuFysNOI91/TWGGPuRzGVz6V0QJAo9ECQLNg2XVykIYKL8tq/faY5RROpLQJSKPRAkCzYNl1apALl9dg9Bwy8Di1BqDRWGgBoFmQDIQTHOsNs3V5bd64yxIA2gms0WgBoFmYtJudv85pDOaNuy0TkHYCazRl9QT2icizIrLHbPX4WXN8lYg8IyKHROQnIuIxxz8iIi+IyF4ReVhEVpS47qMi8qKI7Db/mqb3rWlmmvt3t3PRPz9EJJGa61sZhZXoVe135427tQag0diUowHEgWuUUhcAW4DXishlwBeAryil1gIDwC3m/F3ANqXU+cDdwBfHuPZfKKW2mH/dk34Xmjnhn39xgL5wgicO9c71rYzCFgCBfAHgcugwUI3GYlwBoAxC5kO3+aeAazAWeIDvAW825z+ilLIKrTwNLJvWO9bMGzymOeWhF7rm+E5GYwmAKl+BBuAyvvJaAGg0ZfoARMQpIruBbuAh4AgwqJSydP/TQGuRU28BfjXGpf/LNP/cJrmhGvmvfauIbBeR7T09PeXcrmYW6AvF6RiKAfDU0b45vpvRlDQB2XkA2gSk0ZQlAJRSaaXUFozd/CXAxmLTch+IyDuAbcCXSlz2L5RS5wFXmX83lXjt25VS25RS2xobG8u5Xc0Mk84o/vkXLwBwQVsNg5H5V1htKJrE5RACHmfeuJUHkExpDUCjmVAUkFJqEHgUuAyoERGrltAyoMOaJyLXAp8G3qiUihdex7xWu/nvCPBDDMGiOQv49b4z3Le7g/dfvYZXrGskFE+Rnmc76qFokmq/m0LF0vYBzLP71WjmgnKigBpFpMY89gPXAgeAR4C3mtNuBu4352wFvoWx+Bd17IqIS0QazGM38CfAvqm9Fc1scftjR1jVEOTDr15Hlc/YA4Ti8ysSaNgUAIWICG6naB+ARkN5GkAL8IiI7AWeAx5SSv0C+DjwERE5DNQDd5rzvwRUAD8z7fv/Y13I9CMAeIHfmNfcDbQD356ON6SZWSKJFHtOD/GWra04HUKlKQBGYvPLDDQUTVJVRACAEQqqy0FrNGWUg1ZK7QW2Fhk/ShGzjVLq2jGutcX8NwxcNKE71cwLTpqdtFabCVaVZpTNSGweagABT9HnXA7R5aA1GnQmsGaCHO81BMCKOkMAVM1TATBUwgQEhgagTUAajRYAmglywqyxv7w+ADAhE1A8lSaeSs/czeVgCIDiCq5hAtIagEajBYBmQpzoj1AbcNu766wAGF8D+PBPdvPhn+wed14h/eEEe04Nlj0/HE8xHEtR4y9hAtJOYI0G0AJAM0FO9IVZUZ8tsGb5AIbL0AAOd4c4eGZkwq/5Nz/Yzpv+849EE/naw2cf2M+7v/fcqPk/3X6KdEZx9Ybi5aU8TocOA9Vo0AJAM0E6BmO05jRYmYgG0B9O0DNcNC1kTPa1G529dhdoAQc6h+3ncvnZ9tNsaavhohW1o54DUwPQiWAajRYAmonRPRyjudJnP/a5nXicjnE1gExGMRBJMhJPTbh6aEuN8Xrbj/fnjQ9GkvSHEyiV3c0rpTjeF+bC5cUXfzDaQuqGMBqNFgCaCRCOpwgn0jRVefPGK32ucTWAkVg2W7h7glpAJG6Yfp47MZA3PhRNkkhnGMlJQhuIJIkk0nltIAtxuxw6DPQsJJnO8L9/tIuDZ0ZrfZrJoQWApmy6R4yFu6ly4gKgL5xd9LuGY2W/ZiqdoXvEmH9mKJr3nFWDqD+UsMdOmXkKYwoAh3YCn410DsZ4YE8HTx6ef8UHz1a0ANCUTbe5cDflmIAAaoMeDnWN8MNnTjIQThQ7lYFIdrxrpHwNoCcUx/LXhuNZJ3AsmSaaNB7nCpfTA4aQaKsLlLymyyk6DPQsxDIzhudZ2ZGzGS0ANGVjawAFJqC3XdTGwTMjfOre57lvd3vRc/tydundE9AAOgaNuW11fsI5vgOr3DNAb861Tw8YGkDrWBqA06FbQp6FDJufeWgedqA7W9ECQFM2pUxAb9u2jAuX1wAQSxZfWHM1gO4JaABnzJ4D5zRW2L4AIK8EdX84VwBEqfa7RzWCycXt1E7gsxFLA8j9HmimhhYAmpIopfIibLqHY3hcjqJ9du957xW4nVIyGqg/bIxX+lwMTaB/wGDUWNyX1fpJpDMkzPDNwRyB0hfKCpTOoShLa0rv/o371Sags5HhqLHz1yag6UMLAE1JbrrzWT5zX7ZK96mBCM1V3lE19sEos1zlc9tqeiH94Tg+t4P6oMe23ZeD9WO3/A5WMthgzuv05WgAPSPxURpKIS5tAjorsTYX8630+NmMFgCaokQSKZ462scjB42WDsl0hscP9XL56vqS51T53SWjgU71R2mp9uNzO4kkJiIAjLmN5qJu+QEsLcLndtiOXzAEQEPF2ALA7dAawNnIsPndmsj3RzM2WgBoirLn1BDpjKJjKMaZoRjPHe9nJJbiVRubS55T5XOVNAG92DXC+uZKAh4nsQloAJFEioDHSYXXyDi2NALLNHT9uS08fqiHcDyFUoreUMIWFqXQ1UDPTmwnsNYApg0tADRF2Xkym3S16+QATx/txyHwsjUNJc+pLGECiibSHO8Ls6GlEr/HOaFM4FA8TcDjIug1evuGzd3fQCSJ2ym8/ZLlxJIZfnegi+FoikQ6M64AcDl1ItjZiA4DnX60ANAU5fnTQ7TW+PE4Hew+Pcip/ggt1X6C3tI9hKr8LltNz+WlrhGUgg1LKvG7XURLRAoVI5JIEfQ6CXiM141YGkAkSbXfw7YVtbidwsEzI/SEjIih8QSAxyk6CugsxHICaxPQ9FFOT2CfiDwrIntEZL+IfNYcXyUiz4jIIRH5iYh4zHGv+fiw+fzKEtd9rYi8aM77xHS+Kc3UOdEfYV1zBec0VXCwc4ST/RGWj5FcBUZzmGJ9AV7sMiqArl9Shd/jJDoBDSAcTxP0uAiaAsDSAIaiCWoDbhwOoSbgYSCcsMNLGyqKl4G2cDkdecXg9rUPseWffmvnEGjmJ9oJPP2UowHEgWuUUhcAW4DXishlwBeAryil1gIDwC3m/FuAAaXUGuAr5rw8RMQJ/CfwOmATcKOIbJrqm9FMD0opTppln9c3V/BSV5kCwO+2d2m59JgLc0u1j4DbOaEoIEsDsExAlvloMJKkJmCEo9YG3AxEEvbrjB8FJHnloL/12FEGI0kePtBd9n1pZh/LvGj5ezRTZ1wBoAxC5kO3+aeAa4C7zfHvAW82j99kPsZ8/lUyOm7wEuCwUuqoUioB/Ng8TzMP6AsnCCfSLK8LsH5JFZ1DMXpG4rTVjR1fX+l1EU2m7Vh9i6FoEo/Lgc/tNH0AEwsDNXwAhgYQKjABAYYGEEnaAqCxwlf8YiaeAiewlUfQMRgtdYpmHmBFmKUySofxThNl+QBExCkiu4Fu4CHgCDColLK2e6eBVvO4FTgFYD4/BBTGDtpzipxf+Nq3ish2Edne09NTzu1qpojV+H1FfYANSyrt8bHq64ChAcDo9pBDkWx/Xv8Eo4DCiTQVXhcBj6kBmGGhg5GErQHUmSagnlAcj9NBVYlWkBYuhwOlIJ1RpDOKvaeHAHihU1eZnK8opRiIJHA7jb1kWGcDTwtlCQClVFoptQVYhrF731hsmvnv6Cyh7HNMYI712rcrpbYppbY1NjaWc7uaKXKyzxAAy+sCbGmrscdXN1SMeZ618BbmAuQ2aPe7nSTTquwwzEjcCAMN2D4AKww0SY15zdqgm4FIkt4RIwS0WKJaLl638bWPJtMc6QkRiqdwOYQDncW7lSml+KcHXuDrjx7W4aNzxHAsRSSRZlWD0Y1ORwJNDxOKAlJKDQKPApcBNSJibbWWAR3m8WmgDcB8vhrI7+SRM6fI+Zo55oQpANrqAtQGPez5h+v4/rsu4bxl1WOeZ9XfGSwIBc0VANZOvlw/QCieIuh14XQIfjOJLJ5KE0mkbQ2gJuBhMJKgeyQ2rgMYoD5ozOkLxTncbVg3b7hwGb2hOIe6RguBZ4/1850/HuOLv36RN3ztiQmVs9ZMD5Z5bm2ToZGGdUG4aaGcKKBGEakxj/3AtcAB4BHgrea0m4H7zeP/MR9jPv97Ndpj8xyw1owk8gBvN8/TzANO9kdYUuXD5zYW6+qAm5evG1/7WlFvmIisRdUiVwBY1yzs71sMpRSRRNp2AAe9TsLxlF0JtCZgLOR1AQ+pjOJYb3jcEFDIhon25giA91+zBpdDuHvH6VHz73ziGHVBD1+7cSsn+yN84Ee7tBNyluk0e0GsaTK0UK0BTA/laAAtwCMishdj4X5IKfUL4OPAR0TkMIaN/05z/p1AvTn+EeATACKyVEQeBNs38H7gNxjC5KdKqf3T97Y0U+Fkf5jl9WPb+4uxqqGCgMfJvvahvPGiGkARAaCUyvthJ9IZUhllm38qvEaegVUGIqsBGP+eHohOSAD0jBgCYFmtn7a6AFdvaOLuHafzfBTDsSSPvNjNn25t5Q0XLOXvrlvPM8f6Rwk5zczSbpYF39hiaABDJWpOaSZGOVFAe5VSW5VS5yulzlVK/ZM5flQpdYlSao1S6m1Kqbg5HjMfrzGfP2qOdyilrs+57oNKqXVKqXOUUp+bqTeomTgn+sYP+SyG0yFsaqlif0e+ABgu8AFAcRPQd/54nM3/32/sfgGWo88qA9FU6aN7OGabmGrMKKDaQNbs0zhOHSAYLQCsXeVfXbGSvnCCDbf9mt+90AXAwwe6SKYV15/fAsBVa41M6OcLhJxm+tl5coB33PEMvaE4nYNRXA5hbbMhAAbCWgBMBzoTWJNHNJGmeyTOikkIAIBzW6vZ3zFs9/9NZxQj8VReFBAUz+a0zC9dZs9gSxuwtIamKi/dI3G765idBxDMlqduKEMDqAt4EIEzwzGO9IRY02gIgMvPqeeiFUYz+e89dRyAnzx3itYaP1uWGc7w1Y2GlmNFDmlmhqFIkj/9+pM8cbiXJ4/00TEYpbnKR0PQ+Hxz+0toJo8WABoAjvWGuXfXafaeHgSYlAkI4OKVdUQSaZ47bvj9reSdQg2gWCioFWGTNu3rlpCwcgCaq3x0DcfsBjC1pjN3ZX3QvkY5GoDLaZSlfuJQL/FUhvPNSCcR4Yd/fSl/tm0Zu04O2jWQbr5iBQ6HEVnkdAibl1bZ/09AyQJ4msmz42Q2buRA5zAdQzFaa/xU+lw4JL8hkGbyaAGwwGkfjJbVgOV9d+3kwz/Zw6fN+v8rchbViXD1hkb8bicP7DGCuoYKBIBd06eIBmAlkIXMMNI95iK7pNpI7Gqu8hJJpDnaGwayi319hdcWLDWB8aOAABoqvOwxd/Fbc0JdvS4nV65pIBRP8blfHsDrcnDjJcvzzr38nAZ2nhzkdV99nI/8dDdbPvtbvvzQSyUdw7/Zf4ZvPHokr3GNJksqneHffnWAm+58xtbudp4YxOkQVtYHDAEwGGVpjS9b+kNrANOCFgALnCs//3uu+fdHx503EjcW6sPdIURgbdPYMf+lCHhcXLOhyS6rUCgA/J5sDH4h8ZQxFoonUUrxnSeOsWFJpb1AN1cZgmDv6UHqgh48ruzX946bt7G8LmA7CcfD8gM0VHhYVtA/+LLV9YjA9hMDvHxdI5UF7SX/9zVruPnyFRzoHObnO40eyP/x8CE++8ALeZ3KADIZxQd/vIsv/Pog33/qRFn3tth4+mg/3/rDUR4/1MtTR/sA2HVqgI0tlVy4opZ97cN0DcdoMTu91QTcWgOYJrQAWATkdswqRTKV3b2ubgiOWfVzPDYtreLMcCwvZLM6YAkA47rFCsJZJZpD8TTDsRQHz4zwpi2tdmKX1RVsX/vwqHo/V65p4LG/v7psDWC1mVB0+TkNoxLHmqt8XLfJ6HtwzYamUee6nQ4++6Zzefijr8DrcvCp6zdy02Ur+O6Tx7n+q4/zq+c77SinY31hu0/ywTM607gYj76YrcG0r93oQ7H75CBb22rZ1FJFbyhOMq3sVp+1WgOYNib/K9fMe1I5WauxZNqOwS9Gbm2Vc1vHTvgaD8smf7wvPNoEZN5DYSr/4e6Q7fQNxZL2cW0gu/turjIW/VA8RVPV2PV+xuPTr9/EzVesLFne4gs3nM/apmO8acvSktc4p7GCpz/5KmoCbkSEt1zYyru++xzvvWsn737ZKj7zJ5tsX8HqxiAvnimeabzYefSlHq5a20BfKMG+jmEOdY8QTqS5cEUNzZXZz3mpaQqsDbjtsFDN1NAawAImtzb/eItPMpVhXbNh9jl36dQEgJUQdqIvMkoABAs6exn3meTaL/+BuOUDiKfsqp+BHE2kOWfRby4j2mcsPC4HqxsrcDuL/wRqAh7+7jXrbZ9FKWqDHluDuHB5LY/9/dW8bE0DP3z2JIORBHtPD+F3O/mT81o40R+ZUDOcxUA6ozjSE2JrWw3ntlaxv32IHSeMZkRb22rZ2FJlz11qm4A8o0xtmsmhBcACJjdZ5sA4hc7i6QyXr67nrRcts2PeJ8tK07xyrHe0BuBxOfC6HHk13ftC+T/mkXjK1hCCnqzWEvS6WG/GgTdVTU0AzBRVPjefeN0GIok0v3y+k+dPD7F5aRWbllahFLzUpRPIchmOJlHKEKRXrmmgL5zg0/fuozbgZkW9UYqkxdz5L622TEBubQKaJrQAWMDkCoDcxumFKGUUZ6vyu/k/b7uA1pqxyz6PR4XXRUOFlxN9YYZzSkFbVPpcjOQIgMIfcziesjWEQl+EVY9ovJ35XLJ5aRVtdX5+u7+LfR1DnL+shq3La3EIfOPRw/xyb+dc3+K8wU7qC7h54wVLeecVKwF4xbpGW7Pa1FJF0OO0iw3WBDzEkpkJVZXVFGf+/oo0UyZXAIzlCE5lFEoZdfKni9WNQV7qCrFhSaW9+7eo8LrsUE9glDofiqXszl/BgoX+1ZuauXvH6XGbvswlIsKrNjTz3SePA3D+smqaq3xcs6GJ3+zv4jf7u3A6LuS1505N01oIWMK/JmCY0v7xjZv5xOs25H0Xb7lqFVetzTrrrczvgUiCluqpbVYWO1oDWMBYAsDpEPrDpWPQrQQst2v6vg5b2mp4oWOYnpH4aAHgc+WZgKy0/stX19Na48/zAViF4Cxes3kJ973vSm64cNm03etM8NaLsvdnaS0ffvU6Xrm+kQ1LKrnt/v1F6yElUhl+8PSJRZMzYAn/mpzviM/ttBPvAK44p4F3XrnKfmwFBuhyEFNHC4AFjCUAVjUE7ezZYlghoNOpAWxtqyGRzvDs8f5xNQBrF/jNd1xkC4BQCRMQGMIld4GYj5zbWs21G5vwOB2sMqOiNi+t5rt/dQn/9KZz6RmJ84Onj4867+uPHua2+/Zx9f95lE/f+/yCrzpqxfPXlhm+C9lkP+0InjpaACxgrDIMK+uDoxytucTTxk50OjWAC82aOiOx1CgBUOlz5/kABiNJHGL4BiztwOr8NZV8hLnmWzdtY8dt144SVpesqmPz0ir+8FJ+h7uhaJKvP3qELW01DMdS3PXMyTF9NwuBgYLKruVg1X4a0MlgU0YLgAXMYCSBz+2gpdo3pg/ASsDyTqMG0Fzls+O2RwkAr4tQPPvjHYgkqAl4cDjE1g4sDcA/Ru7CfMfpkFFZxBYbW6pGRQQ9crCbRCrDbX+yif95/5UAC77o3GAkgUi2mVA55PoANFNDC4AFjFWHvy7oYSiaLNnO0KrB43ZNr1llq6kFFPUBxPI1AGsHGPS6CMXTRBIp/G4nznlu6pks65sr6RmJ26a53lCce3Yazu2tbTVsWFKFx+nIKzq3EBk0+0VPxKRnfVe0CWjqaAGwgLEEgNUmsdSOyRIMHuf07ratGj5VxXwA8ZRt3x6IJOxdnVHnJcFILHVWm3/GY90SI5/hpa4R7nj8KNv+5Xc8fqiXv7zcqDzqcTnY2FLJrpMLWwDkfvbl4nU5CXic2gQ0DWgBsIAZjCSp8XuoM2uol/ID2BqAc3p32xeOoQEk08rO/B2IJO3IjqU1flIZxfG+8KgIoIWEldC2+9Qg9+xsp9Lr4u73XM77r1lrz3nl+iaePd4/bhLfRDnWG543Jaxztb+JoOsBTQ/l9ARuE5FHROSAiOwXkQ+a4xeIyFMi8ryIPCAiVeb4X4jI7py/jIhsKXLdfxSR9px51xfO0UyN/nCCuqCHOrNufqlIoMQMhIECnNdazTuvWDmqoFqlubO37PyDpg8AYJmZhHaoKzQqB2Ah0Vzl5ZJVdXz+Vwc50DnMB69dy7aVdXlz3nXlKiq9Lr792NGyrjkUTdqNeEoxHEvyhq89wQd+tGvS9z6d5LYLnQi6Iuj0UM4vPgV8VCm1EbgMeJ+IbALuAD6hlDoPuBf4GIBS6i6l1Bal1BbgJuC4Ump3iWt/xZqrlHpwyu9Gk8dAJEFdhSevBWIxLA1gOp3AYFTN/Mc3bmZVQ35vgQqfKQBMP8BILGU7Aa16L33hxILWAESEr/z5Fja1VOF3O3ndeaOTwqoDbl6+rpFnjvUXuUKWz/3yBa78/O+5+HO/4+bvPFs0v8Dip8+dIhRP8eiLPTx5pHfK72MslFL86vlOu8x3MWLJ9KQc/XVBrQFMB+X0BO5USu00j0cwmri3AuuBx8xpDwE3FDn9RuBH03Ori4upprlnMoqBSJK6gMduqNI5VLyC4kwkgo1FpddY7EdiKdIZRSieotIUCq05tfkXsg8AoLXGz4MfvIr9n31NyfIbW5fX0D4YpWu4+Gd3oHOYbz9+jPbBKFetaeCJw73ct7vdfn4ommTHiQH7M/7Z9tNcsKwal0P442FDAAxGEkW1w1P9kSk1X9/fMcx779rJ/bs7Ss6JpzJ4J/G9MwrCaQ1gqkzof15EVgJbgWeAfcAbzafeBrQVOeXPGVsAvF9E9orId0SktsRr3ioi20Vke09PT7EpC46njvSx4bZf8+w4O7+xGI4Z5oDaoIcKr4tKn4szQ8Vjyi0NYDoTwcbCWuyHoknbDGSNVXhd2cqhC9gElMtYETCWH2XXSaNCZiaj8sp8/2z7abwuB3v+4TruuHkbAY/Trvz66IvdXPavD3PDN57k7h2n6RiM8mLXCK8/v4WVDUE7DPU9/72Dv/7+9rzXPXhmmKu++Aifvvf5vPHD3SMlNclCrPt4oaO0DyOeSuN1TVwDqA966C3zPjSlKfsXLyIVwD3Ah5RSw8C7MMxBO4BKIFEw/1IgopTaV+KS3wDOAbYAncC/F5uklLpdKbVNKbWtsbGx3Ns9q7H66f7+YPc4M7Mk05m8ME9rR1dv2v9bqn3jawCzJACyGkmUEdMZaQkAyAqk+VrxczbZvNQwEf3uQDdKKd73w53c8M2n7OefOtrHtpW1VJs9Cc5prOBIT4jOoSjv/+EuuzLr3tODPPqisYG6en0T65orONQ1Qn84wTPH+tl5ciCv/MS//OIAAM8c6yeSSPHrfZ0opbj2y4+V1WEO4HCPIWDGEgCJVAave+Lfu6U1PkbiqXnjzD5bKet/XkTcGIv/XUqpnwMopQ4qpa5TSl2Escs/UnDa2xlj96+U6lJKpZVSGeDbwCWTeQMLkQrT9DEygS/3x+/ey/t/uNN+bNlHa20B4OdMCTOCFY3jmSUTUD4L93oAACAASURBVEu1HxGjX3FWA8g6Am+6fAWXrqrjQ69aNyv3M5/xupy8bdsy7t/dzv/9/WF+te8Me04Ncqo/wrcfO8qBzmEuW1Vvz1/TVMHh7hD//0OHSGUyfOsdF3Hlmnpe6Bjm0Re7aa3xs6apgrVNlZzoj/Db/WdQCpSC3x3o4n/2dBCKp+xNyFA0yR2PH+M9/72T377QBRimu7Hs+haHuw0BcKBzuGRJi3gqMynNs7XG6DnRvsAzpWeacqKABLgTOKCU+nLOeJP5rwP4DPDNnOccGGahH49x3Vyv11swTEoaslE5I7Hym4e82DXCyf7sj8EK+awLlKMBTH8toLHwuBw0VXo5PRC132OuBvCp6zfyk7+53G4judj566tW4xDh3x96yW62c9UXH+FzDxq79KvWZTXjNU0VdA7FuG93O6/ZvITl9QE2Lqni+fYh/vBSD69Yb5RZXtdciVLwiZ8/T0u1j6ZKLx+/53k+8KNd/O1dO4mnMrxmczOJVIbbzSikr/7ukP06zx0bGPe+j3SHcDqEkXiqZEmL+BQ0AMgXAD946jib/uHX/P3deyZ8vcVKOf/zV2JE81xTELJ5o4i8BBwEOoD/yjnn5cBppVRe/JqI3CEi28yHXzRDSPcCVwMfnuqbWShYjre+MSp4FtIbihPPcRxbGkCdmQS2pNpHbyhum1dysRPBZkkDAMMB2j6QNQFVLHCH71Roqwtw580Xc8U59Xz/XVlF+Ws3buXBD1zFFjPhDrL5BfFUhldtNPoab1paRUYZY1evN0JyX7G+kYtW1OJ2GtFId958Mee2Gt23HjNrFL3jshVANlz3hZx8hCdMB/I9O07ztYezgsEinVGc6I9w7lLjmsU2H6l0hnRGTcoHYAULtA9mBcBvX+gikkjz0+2n2d+xsEtoTBfj/uqUUk8ApbxUXy1xzqMYIaOF4+/OOb6pvFtcfFgCoFz1VilFXyiBozL7MfWbpXJzNQCloGs4NqoP7kwlgo1Fa22APacGczQAvdsfi5etbeBlaxsA+NJbz8fndvKGC0b3K37l+kbecMFSnjrSxytMzeDl6xq5ck09PSNxrjjHMBdVeF387G8uN4r1mZrWL/73VZwZivEfvz/E0mofF6+sY2V9gPoKL2/Z2spn7ttnZihXsdN0Sn/0Z8Zu+8q1DVy4PBvH0ReOk84oNrdWs+f0EL1Fyltbmu5kooAagl48LoctAJRSvNAxzGs3L+GPh3v5z0cO8/W/uGjC111s6G3XPMQWAINRMhk1bp2UoWiSVEYRzdEA+kJxfG4HfrOl4hKzccaZIgJgrjSAX+/rtN9rlU9/FcvlbduKBdwZuJwOvnbjVtIZZddRaqjwcte7R+3HcDhklJltSbWPf33LefbjRz92tX28aWkVAty/u4MfP3cyr7/xd544xoX/KysArEihjWbJi2L9DeLJyX/vHA6xtUiAruE4feEEl62uY01TBf/3kcNc+fnfc9e7L7Ud4ZrR6FIQ8xCrjHMyregvI9ml17T35yYAdQ7F7B6qgN1XtZgqHk/NbhQQQFudn2Ra2aGCWgOYXmaiiN6Fy2vZuryWC1fUEktmuG9XNr6/MNLH+k6ua65EBHqKlCGxvneTMQGBaUY0NQDL5LO5tZp3vWwVDRUe2gejfP+pE5O69mJBC4B5SG7yTTg+viPYUq/jqQwZsxRAx1CUFtNRBtnQy2K5ANlicLP3ddi81OiS9dTRPpwOwTcJR6BmbrjinHoCHiefMnME3rxlKcf6wkQTaTv82NIAllT7qA14imsAZiTRZExAkC8ATvRFADinsYK6oIftn3k1rz+vhXt3nS7q99IY6F/dPGQomrR/FOH4+OF2uUXerF1V52Asr19qlc9NhddVVANIpDK4HDKrXbY2LKnE5RCO9oSp9Lnsfq+a+U9DhZePvNoI0V1W6+c1m5egFPzd3Xu48J8f4lR/xBYADRVeGio8RX0AtgYwSeHfWuunZyROLJlmMJpEJL/w4Ks2NjEQSXKyPzyp6y8GtOF1HjIUTbK0xs+x3jDhRPkaAEA0mcblFLpGYnZdHYsl1T7OFBEAyXRmVs0/YPR9Xb+kkv0dw3khoJqzg1tetorrNi1hSbWPTlOr/OXeTsCIFuoNxQl4nAS9LuqD3qKVaKeagW6Vz+gcijEUSVDlc+eZvqwN0JmhOGuaKif1GgsdrQHMMzIZxXA0acc5h8owAfUVCICu4RhKYXfksiiVC5BIZWbVAWxhRY2UW1pAM38QEZbXB/C4HLTVBrjKjFACONoTpmckbhchbKj0ltAATBPQJLu+2aGgA1EGo6PLSltmz1J1lDRaA5h3hBIpMgrbgVuODyDXwRZNpO0cgFEaQJWPQ12jK0Am0mrWNQCAD127ln0dQ3bsuubsxOEQfnDLpfSF4rzuq49ztCdEbyhOQ4UpACo8tlM4FysKaCo+AID2wYjZ+6JAAFSZfi8tAEqiBcA8Y9CM37d2N5GyfADZ3VUsmabDdIy1FGgAbXUBukZijMSSeVE3iUlWZJwq9RVe7v3bK2f9dTUzQ32Fl3MaKzjcE6J7OM4FbYajv7XGTyieons4RlNV9juZjQKa3HdvSbUPh2Q1gOqCzmJ+j5Mqn4vuHAGglNL+phy0CWieYWX/Ljdj9csxAfWG4limz2gybZ9T2GhjS1sNShldqHIxfAD6R6GZOuc0Bdl1cpD2wSjXmpnIl5q1ip462pc3d6o1qNxOB81VPk4PRhmKJEZpAADNVT5bAzjcPcKqTz7IE4dmtg/C2YQWAPMMK4zOStYqxwTUF07YDq9YMk0sWTy+euvyGhwC24/n13GJp9Jz4gPQLDzeeEGrffy6c41yX5uWVlHlc/HUkUIBYIWBTr7xj5UMNlTEBwBm4MOwsal61qxf9N9P69wAC20Cmmf0mQJgSZUPj9NBqJwooJE457ZW0z4YJZpI5zjX8hf1Sp+b9Uuq2HEiXwCE42ldi0czLVyyqo7733cl4UTKzkJ3OoSLV9aN+t5N1QQEhql0+/GBkq0lm3P8XlYOTE8Rh/RiRW/75hmWBlAX9BD0Osf1AUQTacKJNMtqDY0hmqcBjP54Nyyp5Fhvflz0SCxJhc7E1UwTF7TVcMU5DXljVjHCXKaaBwDZZLCMGm3yBFjdGOTMcIyekThHeozv/Qsdw3m9MxYzWgDMM/pCcbwuhx1DPZ4JyPpRLavNmoDiyTRel6Oos6u5ykf3SCyvPvtITktGjWYmqA14GIom7Ux1yOlFPRUTUE4L0ZoCJzDAy9caBfHe98Od/PJ5I08hmkyXLE+92NACYJ7RF05QH/QgIlR4XeM6gS2TkSUADBNQ6aie5iqvUWMopwfsSCxFpTYBaWaQmoCbjMrvcTHVUhBAXi/lYk7gTS1GOWqrvao1f3gKvY4XEvpXP8/oDyfsGv4Bj3PcTGCrL6rlNI6lMsSSaXwlkmtyY6PrzTjtUExrAJqZpdbcnQ9EEnYFUrsa6BRyUJblaACbzN4DuTgcwm1/somTfWGiyTTrmiv5l18e0K0kTfSvfp7RH05QFzQW5qDXNW5XMCvEzdrZWBpAKQHQnJMduXlpNcl0hmgyTYVX+wA0M0dt0Ph+DUQSBIadNFZ67XaQU6lBtbwuyGWr63jnFatGJT5a3PKyVfbxwTNG1dLhaPnd9hYy5bSEbBORR0TkgIjsF5EPmuMXiMhTZlevB0SkyhxfKSLRnO5h3yxx3ToReUhEDpn/1habt9joCyXsRu4VZfgAdp4coD7oYVmtH6/LYYaBpkuq1ZYG0GWGxoXjo1syajTTjWWfP3hmhEv+9WE+c98+4qnS39Ny8bgc/PjWy3ntuUvKml9lBjtMpN/2Qqac//0U8FGl1EaMLl/vE5FNwB3AJ5RS5wH3Ah/LOeeIUmqL+feeEtf9BPCwUmot8LD5eEY41R8ZFYI2XxmMJOx45rGcwMl0hrff/hQ/39nOZavrERH8HifR5NgaQGOlFxHsonCWhlGhBYBmBrFMQHtPG0mIdz1zkhN9kVnPP6ky/QTaBGQw7v++UqpTKbXTPB4BDgCtwHrgMXPaQ8ANE3ztNwHfM4+/B7x5gueXzTf/cIRbv799pi4/bWQyinAibTtk64MeesOJvMgJizNDMZ4+aji2tq00lKeA20k4PrYG4HY6qA967QJZ1g9Bd+TSzCS15qbGagAEsK99aNZLkAQ9ThximICUUnzpNwd5/vTi7R88of99EVkJbAWeAfYBbzSfehuQ26dulYjsEpE/iMhVJS7XrJTqBEPIAE0TuZeJUBNwMxhN5oU+zkciZkvHoCkAltX6SaQyRSspWmNBj5PXn29kXFYHPAxFE2NqAGD8f1g7/5ClAWgfgGYGqfS5EYFDXSF7rHskbu/IZwsRodLnZjiWpGckzn8+coQ3/N8nyiq5shApWwCISAVwD/AhpdQw8C4Mc9AOoBKw4go7geVKqa3AR4AfWv6BySAit4rIdhHZ3tPTM6lr1Pg9pDNq3n/IlrknKwCMyJ7Tg6Njlq3qij+69TKaKg27fl3QTX84MaYGANi+AiCnKbvWADQzh9MhVPvdjBT8Bgv7U88GVX4Xw9FkXkLkIwe78+bsODFAahEki5UlAETEjbH436WU+jmAUuqgUuo6pdRFwI+AI+Z4XCnVZx7vMMfXFblsl4i0mNdvAbqLzEEpdbtSaptSaltjY+PE3p2JFXY2GJnfdr+sADB271aIW7GkFasCqBXKCYajbTCSHDMMFIxmLDEzBtsSitoHoJlpLD9AS7XPLhW9Yi4EgM/QgK02kgAn+7PH9+1q54ZvPMnPd7XP+r3NNuVEAQlwJ3BAKfXlnPEm818H8Bngm+bjRhFxmsergbXA0SKX/h/gZvP4ZuD+yb+NsbFSxIfmefKH1f4x6DEW41ZbAERGzbUSwKyIIYC6gIeBSGLMRDAAn9thl4uwoiG0BqCZadY0VQDQVOm1Gx4tr58bATAcS3KsL4zLIdQFPZzMEQZ3PGEsV+2LIFu4HA3gSuAm4Jqc0M7rgRtF5CXgINAB/Jc5/+XAXhHZA9wNvEcp1Q8gIneIyDZz3ueBV4vIIeDV5uMZoWYGBMDh7pFpjyyyd+OmCSjgcVEf9BTVAHpDcSq8rrydfq3p64gm0mN2WfK5nFkTkBUGqn0Amhnm0lV1ACiy4cjL58wElOJEX5jldQFWNQQ5YfYNDsdT7Gs3cgXai5heFxrjbvuUUk8ApTI1vlpk/j0Y5qJi13p3znEf8KrybnNqWDHI02kCevvtT9MbSrDnH66zTUxTJZLI9wGAYQY63B0aNbc3lKChIr/2SW3Qg1KGdjC2BuC0C3ENRZO4nYJvCgW5NJpyuMQUAMd6w2xtqwHmRgBYTuBjvRFW1AeoCXjsUhG57SNztYKFyqL41Vtx9YPR0W3pJksybUQU/fcz01dbPFTgAwC4ZkMzzx7r52hPVgh87eFDPLCnI8/+D1kbKzCmD8DrzjqBB8IJ6szaQxrNTGLV5XnPK87h3NZqGiq8dqDDbFLlczMUTXKiL8zKhiBtdQE6hqLEU+m8zHpLK1jILAoBMBM+gLY6wz6/r336YohtH0COBnDjpW24HMLdO07bY//+0EtANpnLojbHHzB2FJDT9gH0h5N5gkOjmSlcTgfH/u163nf1Gt560TKe/uQ1c9KIqK3OTySRJpJIs7I+yNqmCpSCrz9yhG4zQ/6SVXV0DcftjdJCZVEIAJ/bidflYGgaTUDWYm2VtJ0OipmAmip9LK3x5/kBVpiOM0ultqjNMUWNHQXkIG5+sfvDceqCWgBoZgdL0xQRXFMoAjcV1jVX2scrG4K87twlvP78Fr768CFe7DIS1S7NMVeVIpHK8PlfHczrOXy2sSgEAJjJYNMoAKz4+cQ0xgpbJqBAweJdF/TklW/2OB1csrKOf/vT8/Lm5ZuAxvYBWGGgA5GkFgCaRcXa5gr7eGV9AJfTwU2XrQDgwec7qfC6uHCFkV1vFY8rxm/2n+GbfzjCV353aEr3s+fUID957uSUrjFZFo0AqPa7p9UEFIob14pPowYQjqfwuR2jdkb1BQIglkrTWusftcuvr8g1AY0dBZRMK9IZZVYf1QJAs3horPBS7XfjcohdRXdLWw1up3CiL0JTlZfVDUE8LgcHOkdKXuf5aTL/vuk//8jH73m+aMmXmWbRCICaQP4iOhVS6YxtQ59OE1CoRG/e2kIBkMwU3eEHPC78plAYWwMwngsnUgxFtQagWVyICOuaK2irC9ibLZ/byfnLjMikoMeFy+lgfXMlL3QU1wDu29XO7Y8Z+QI9I5M3AeUu+r3h2e9VvGgEwLJaP6eKJFRNhnBOn97p9gEEiwiA+qCH/kjCrmVklHoovsO/oK0aMFLvS2E5iLtMJ7IWAJrFxidet5F/etPmgrEN1Ac9XLXW6Ge8qaWKFzqHR9UQC8VTfPyevfbj3Czif33wAE8e7i37Po7kRPfNRZvKRSMAVtYH6RyKTYtXfySeNSVZbe2mg3A8RcAzWgDUBT0kUhnCCeO14snSxd4uWVUPGOGdpbDOtRJddBSQZrFx0YparlqbX1rm4pV17Ljt1fz9azcAhq+gP5wY5Tv81fOdxFMZ7nnv5dzyslWc7I+glCKWTHP7Y0f5X3c8M+7rf+3hQ/zxcC87T2aTSeci83jR5P9bkTOn+iOszYkCmAyWs9bndky7E7jCO3pht3bo/aEEAbeTRLp0qYe/feU5ZDKKt25rK/o8ZAVAp6kB1GsNQKMZxcr6IADH+8J2iLVSiv9++gSrG4JcuLyW/R3DxJIZukfiEyo2aYVyf+CaNfaY1gBmkBX2h1ncDHS4O8RHfrq7LJOOVUK5PuidVhNQNJHGX0IDAOg36/xA6TBPn9vJ371mfVFfQnaO8bF3mhpAjdYANJpRrGzICgCLp472sef0ELdctQoRsTOZT/ZH8nbwYzl0c5/rGo7TWOmlJuCmfXD2M48XjQBYaWoAJ/qKx/Xe+oPt/HxnOy91lfb6W1iS3jLNTBexZAZfkZ29LQDC2cSUqZRusOoEWW0hrX6tGo0mS1udH4fA8d7swnz39tNU+93ccOEyIFvK4kRfJK92UMdQ/m4+t7R0LMds3DEUpbnKS2tBrs9ssWgEQE3AQ23Aza6Tg0WfP9pjCIaeIs1XCpkxAZAqXsa53mwS3xdK2F+esRK9xsNnOpCt92r1SdVoNFm8LidLa/z2pjGZzvC7A11cu7HZ/v0tqw3gEDjZF87TAA7l1O863D3Chf/8EP/24AGUUnlBJLtODtJc6aOl2jcqs382WDQCAODPLm7jV/s6RxVXy+2721XGh2CZgOqCnmnNAzDq+I/+SBorDQHQPRK3w0+npgE4zOvFcDqEgGfywkSjWcisrA/a2cDPHutnOJbiNZub7ec9Lgct1X7DBDQYtX1zuYXk/uH+/QzHUnzrsaPsOjWYt96E4imaqnw0VvqKdv6baRaVALj1qtW4HA5++Ex+1t2Bzmys75kx0roHIwne/J9/ZK+ZAFIf9JDKqGlL4CjVytHvcdploa2oI98YiV7jYZ3bPRyn0ufSheA0mhJsXmqEgkYSKZ473o8IXLGmIW/O8roA9+3u4N5d7Zy/rBq3U+wAi4FwgieP9Nkmo0NdI4QT+c7ipkovjZVe+sKJWe9CtqgEQH2Fl1eub+SBvR2kCxwxFmOpYfvah9l9apDf7u8CssXXrEighw908YOnJ18ddKxOXq21ftoHozkawBQEgK0BxLX5R6MZgyvXNJBMK547PsDOk4Osb64cFWBhadRLqnz87dVraK7yccb0Aew+bZic37K1FbdTONYbIZLIDx1vrvLRVOm1S7nPJotKAAC8cctSekbiec1crEy+1hr/mBpAp/mh9obiNFV67UXYMgPd8r3t3HbfvkndlxFHXNwJbN1b+0DEdgKPVe1zPHKFR5V/0UQCazQT5uKVdXicDp441MPukwNsXV47as6lZu7NT//mcq5e30RLtc/WAHadHMQhsHV5DW11AY73hm0T0CozysjtFNvM2zMyu2agRffrtz6s59uHONQ9Qjqj6B6J43IIG5ZUjtkFqDNHO1jXXGmXsi10BI/Xk7cYlhAp1cmrtcbPIy92E7UEwJQ0gBwBoDUAjaYkfo+Tza1VPLCnk+FYiq3La0bN+eurVvFn25bZ/Tlaqv3sMXf+u08Nsq65kqDXxar6IMf7wrYG8Lm3nMvPtp/muk1LONpr+CW7R2JA9ey8OcrrCdwmIo+IyAER2S8iHzTHLxCRp0TkeRF5QESqzPFXi8gOc3yHiFxT4rr/KCLtBW0mZ5xG09626+QAn753H/9w/356RuI0VHhZUu0rSwMAI0vQa9YRKcwG7pyENz9umnZK7exba/3Ekhk6B41rT8kJnPMaWgBoNGOzprHCXhespja5uJyOvOZMlgaglOKFjmHOazUW9JUNhgCwgkjaagN85c+3UB1wz5kGUM4qkgI+qpTaCFwGvE9ENgF3AJ9QSp0H3At8zJzfC7zBHL8Z+MEY1/6KUmqL+ffgpN/FBNnUUsUv9nbajzuGojRVeanwuUbZ53IpVwMo1sR9PMYL77SqFloRTFPxAfjdTjym8NImII1mbM4xm9k7JNvYfiyWVPtIpDIc7g7RG4qzwRQay+sCxJIZuyZZbt2vhop5KgCUUp1KqZ3m8QhwAGgF1gOPmdMeAm4w5+xSSnWY4/sBn4jk9y6cY9Y153+IO04M0Fjhxet0kExnRhV/srB23wBrmyrsnbTlBLYetw9ECcdTfPuxo2VHCMXHce5uNL9Ezx3vH3NeOTgcwupGw/6oNQCNZmzWNBrrxcqGYFm/u5Zqo+H9wwe7AdiwxCg9Y5Vrt4rH5bZ+9bmdVPlc808A5CIiK4GtwDPAPuCN5lNvA4oVn7kB2KWUKvWu3i8ie0XkOyIy2rsyQ/zphct41YYm/uPGrYCRgdtU5cXtdKAUpEos2h1DUa4/bwm3vnw1F7TVjNIArN7D7YNRPv+rg3zuwQP87kBXWfeU1QCKfyTLav20VPvsGuSlnMXlYjmgKrUA0GjGxNIArIV8PDYsMTZrP9t+Ku88K6P/VH8El0NsLdyi0ucmFJ/dFpRlryIiUgHcA3xIKTUMvAvDHLQDqAQSBfM3A18A/qbEJb8BnANsATqBfy/xureKyHYR2d7T01Pu7Y7JxpYq7nznxbwipxpgY4UXt7moJovE4saSaUZiKTYvreZT12/E7XSMEgBWo/iT/RHbZpgpoU0Uuz6Uju8XkbwWkFNxAoOhpgK4nDoHQKMZi7ZaP02VXi5ZWTf+ZIzCk3VBD0d6wjRVem3/gJXRf7I/SsDjHJV/4/c4iSbLLyg3HZQlAETEjbH436WU+jmAUuqgUuo6pdRFwI+AIznzl2H4Bf5SKXWk2DWVUl1KqbRSKgN8G7ikxLzblVLblFLbGhsbi02ZNNUBt91H98o1DbhNiZxMjV60h2NGSdgqX9Zu57GdwIYAiJr+g72nh+wFvdym11Z8v3cM5+5lq+vt46lqAM1VhgCYbZVToznbcDkdPP7xq/nLy1eWNV9EbDPQW7a22uOWCag3FC/a9yPgcY7pg5wJyokCEuBO4IBS6ss5403mvw7gM8A3zcc1wC+BTyql/jjGdVtyHr4Fw6Q069zz3it48hPXcOnq+uyOvogGYPUAzjWZ5GoASimiyTRel4NjvWG7sFO5tYKyRd5K7+xfuT4rAKfaUPs1m5cA8Lpzl0zpOhrNYsDrcuIYo8lSIa9YZ/xW32H2Ggaj74a16S9WfsXvnocCALgSuAm4piBk80YReQk4CHQA/2XOfz+wBrgtZ74lLO4QkW3mvC+aoaJ7gauBD0/j+yqb1Y0VLDUjbDymOaSYAAjZAiArua2uXPFUxtYCXm5+8Fb9kGiZDWjGMwGBEV88XaxqCHL886/n0hytQqPRTA8funYdj//91bSZ1ULB6NJX4zc2kKU0gOgsC4BxYwCVUk8ApUTfV4vM/xfgX0pc6905xzeVeY+zRtYENEENIJ2xP7hLVtbxh5d67J1/NFGeBpCt8z+2TP7M6zfmdRHSaDTzD4/Lkbf4W9RXeBmIJIv26wh4XEQSs9sTQAeB52ALgKImIMMHkPvBeXNMQNZOv8rvYkVdwC4HG0mU59QpxwQE8O6rVpd1PY1GM3/ZVsSh7HPPvgaw6GoBjYUlAIr6AOKjTUCWBhBPpW0B4HM77RBLoOwexLHU2JnAGo3m7MdK5LxuU/Oo5wIeZ9km4+lCrzY5eO0w0NFRQJYJKDdxyooCSqSyJiB/gQAo9wONT0ONH41GM7+55WWrAKPMdCHzMgpoMeF25sf152KbgHKdwO7sfGun7/c47Rh7oOwPdDpaPWo0mvnNbX+yiaP/en3RHhx+j5N4KpNXqn6m0atNDm4zCqi4DyBFwOPEmRMKlpsHYO30/W4ntTlN1ss2ASUziDAqO1Cj0SwsSoWTWqGhs2kG0qtNDu4x8gBCsVSe/R+MWPwKr4uBSMI2AfncTq7b3Mz15y1BpHwNIJ5K43ONzg7UaDSLA79p/i03cGQ60AIgB89YYaDxZNG6OQ0VHnpDiawG4HES8Lj4+l9cxPrmyrK9+rFkRpt/NJpFjN9jbDBnMxJIrzg5jJcJXKgBgFHGtXcknvUB5Dhx/RPw6o/EkgQ8OipXo1msaBPQHDN2HkCqaPJGQ4WX3lA8LwrIwj+BuN7jfRFW1I9OHNFoNIsDv8cyAWkBMCfYTuAixeBGYsmitfMbKj2GADCLufk9BQKgDGmulOJoTygvfFSj0SwuAubmUZuA5gjPGIlgQ9EUVf5iPgAjtbt9MILP7chL5CrXBDQQSTIcS2kBoNEsYrQGMMd4SvQDUEoxFE3YDV9ysVq5/WZ/F1vaavKieHJNQH/69T9y1zMnAKMD2f272+15x8yG0FaXLo1Gs/gIeGY/Ckh7HXMolQgWTaZJphXVJTQAn6gKBgAAER5JREFUMOrqv/3i/KZoVmp3PJVm58lBdp4c5C8uXcEN33gSgDdtMWqFv3jGEACrGsbvN6rRaBYmxaKA/u5ne9h+vB+/x8U/vWkzF5fZlKZctAaQQykn8GDEyAKuKSIAGiuzSV+FBZ58Zmr3kHl+IUop+kJxvvzQi6xuDNJWO33lnjUazdlFwD3aBPS7A12ICK01/rwAk+lCawA5uO1+APlOYFsAFDEBbVhSxas2NNFY6eXSVfkCwO92kkhl6A1lu2XmduCKJtP88UgfvaEE33zHRVNu8qLRaM5eAt7RJqBwPMWNlyzn46/dMCOvqQVADiJGo+ZRGkDUWMCr/Z5R5wS9Lu5858VFr1dvNoE+0hOyx544nO1rPBhJcqrfqP+9qUhxKI1Gs3jwupy4nULY1ADiKcP0XCz8fLrQW84C3E4ZlQlsmXCK+QDGwuq7++KZEXvs8Zd67eOBSIKTfREaKrw6CUyj0RDwuAibpecj8bQ5NnMVgsvpCdwmIo+IyAER2S8iHzTHLxCRp8y2jg+ISFXOOZ8UkcMi8qKIvKbEdVeJyDMickhEfiIio7fXc4Db5RgVBjoYLW0CGgtLABzMEQDPHOu3j4ciSU72R1hep23/Go3GaDgVNhf+kCkIirWPnC7K0QBSwEeVUhuBy4D3icgm4A7gE0qp84B7gY8BmM+9HdgMvBb4uogUE2FfAL6ilFoLDAC3TPXNTAfuIiagoUkKAKss9MEzwwC01vhpH4zazw9Gk5waiLC8SOs4jUaz+Ah4nLYGEDZ9AXNqAlJKdSqldprHI8ABoBVYDzxmTnsIuME8fhPwY6VUXCl1DDgMXJJ7TTGC5a8B7jaHvge8eWpvZXrwOB0kUqOdwB6nY8Je+IYKLw6B0wNRPC6HHedvVYPtDcXpGIxqAaDRaABjt28t/OH5YALKRURWAluBZ4B9wBvNp94GWEHwrcCpnNNOm2O51AODSqnUGHOs17xVRLaLyPaenp5iU6YVj6uYBpCgyu+ecKlmp0NorDTyBGr8bpaZYZ5XrmkA4EDnCBkFrTr8U6PRAEFvtiuYpQnMCyewiFQA9wAfUkoNA+/CMAftACoBK9ax2CpZWFynnDnGoFK3K6W2KaW2NTY2lnu7k8btlFGJYIORJNX+yX0IS0w/QG3AYzeK2by0Gq/LYUcHWclkGo1mcRPMcQKH54kPABFxYyz+dymlfg6glDqolLpOKXUR8CPgiDn9NFltAGAZ0FFwyV6gRkRcY8yZE4r5AELxVNFeAOVg+QFcTrFbvTVVeqkJuDnaEwagLjgv/N8ajWaOyTMBmZpAcAYjBMuJAhLgTuCAUurLOeNN5r8O4DPAN82n/gd4u4h4RWQVsBZ4NveaSikFPAK81Ry6Gbh/am9lenA7R0cBRRPpSWfhvecV59BS7eNlaxv465ev5s+2LePPLm6jNmBUEQUtADQajUHQ67Rt/1kNYOZ8AOWIliuBm4DnRWS3OfYpYK2IvM98/HPgvwCUUvtF5KfACxgRRO9TSqUBRORB4N1KqQ7g48CPReRfgF0YQmbOKeYDiCbTE84BsNi6vJanPvkq+/EX33oBAI2VXjs8VAsAjUYD+Sag2QgDHffKSqknKG6zB/hqiXM+B3yuyPj1OcdHKYgOmg94nI5RJZyjiTS+afbEWw5hj9lXWKPRaIJeF/FUhlQ6QySRwumQvBLz043OBC7A43IQTxUIgGTaLtQ0XbTWGALA79GN4DUajYEV8hlOpAnH0wRneH3QAqAAr8tBPJlvAook0nmdvqYDK/QzNov9PzUazfzGsgaE4ylC8eJtaKcTLQAK8LmdxIpoANMuAGqM5K94anT3MY1GszgJmAt+JJEikkjZj2cKLQAK8LkdxHI0gHRGkUhlpr0Wt07+0mg0hVSYET+heJqRWGpGHcCgBcAovC4n8RyzjOUQnu507GYzQ7ihQkcAaTQaAyvmPxRLcaw3PONlYnT4SQFet4NYjlnGas823RqAy+ngWzddxMYlug+ARqMxaDIrBxzrC3N6IMqfb2sb54ypoQVAAT6X0cUrk1E4HJIVADOQjfeazUum/ZoajebsxSod88Qho+7Z2ubKGX09bQIqwOs2G8ObyWCWCWgm+nFqNBpNLn6Pk9qAm8fMxlFrmytm9PW0ACjA5zIWeis80+rPOZMlWTUajcaipdpPNJnG43SwYoZ9AFoAFOBzWwIgXwPwaQ1Ao9HMAi1mAcnzllXjcs7sEq0FQAFW2rWlAWR9AFoAaDSamceqO3bZ6roZfy0tAAqwdvpWgtZMhYFqNBpNMXrMKsHnL6uZ8dfSUUAF+Nz5GkBkhsJANRqNphh//5oNOB0v8vK1M98ASwuAArwFTmDrX20C0mg0s8F5y6r57l/NTqFkbQIqwNIAbBOQ1gA0Gs0CRQuAArJRQNoEpNFoFjZaABRgRwHlOIG9LgcOh67Zr9FoFhbl9ARuE5FHROSAiOwXkQ+a41tE5GkR2S0i20XkEnP8Y+bYbhHZJyJpERkVzyQi3xWRYzlzt0z/25s4dhSQqQHMRk1ujUajmQvKWdlSwEeVUjtFpBLYISIPAV8EPquU+pWIXG8+fqVS6kvAlwBE5A3Ah5VS/SWu/TGl1N1TfxvTh1UKwtIAIvGZL8mq0Wg0c0E5PYE7gU7zeEREDgCtgAKsUpbVQEeR028EfjQ9tzo7WFFAWQ0grXMANBrNgmRCPgARWQlsBZ4BPgR8SUROAf8H+GTB3ADwWuCeMS75ORHZKyJfERFvide81TQxbe/p6ZnI7U6KwiigsDYBaTSaBUrZAkBEKjAW8w8ppYaB92KYd9qADwN3FpzyBuCPY5h/PglsAC4G6oCPF5uklLpdKbVNKbWtsXHmEyM8Tgci+cXgtAlIo9EsRMoSACLixlj871JK/dwcvhmwjn8GFGYuvJ0xzD9KqU5lEAf+q8j5c4KI4HU5bAEQiqcIerUJSKPRLDzKiQISjN39AaXUl3Oe6gBeYR5fAxzKOafafO7+Ma7bknP9NwP7JnrzM4XP7cwxAaXtNm0ajUazkChnZbsSuAl4XkR2m2OfAv4a+KqIuIAYcGvOOW8BfquUCudeSEQeBN6tlOoA7hKRRkCA3cB7pvROphGfy2lnAIe1CUij0SxQyokCegJjkS7GRSXO+S7w3SLj1+ccX1PWHc4BVX4XQ9EkSinC2gSk0WgWKDoTuAh1QQ/94QSxZIaMQmsAGo1mQaIFQBHqK7z0hxOEzXaQOgxUo9EsRLQAKEJ90ENvKE44bvUD1gJAo9EsPLQAKEJ90MtwLMVgJAlAhfYBaDSa/9fe3cVIddZxHP/+wIVFoMu7AgUKBiPE6IorQWvQVKMtN2jCBTfKRZUENbGNJtI2Ia3RizZqo8bY1FjrW6Ra33pRo6St0ZtSaeU1lHaVqggWawstaa1F/148z8CwnLOdhWXOmZnfJ5lw5pmzs7/9szP/Pc+ZnKcLuQEUmDVtEgBHnnsJ8DkAM+tObgAFZk9NDeCvz74IeArIzLqTG0CBkQ3AJ4HNrBu5ARSYnaeADj9zCoDp/W4AZtZ93AAKzJ6aLkz6+D9eAGDOtMILlZqZdTQ3gAIDU/qY9JoJnHjxFWa+Nm2bmXUbv7MVmDBBLBjoB2De9P6K05iZXRpuACUWzpwCwLzLPP1jZt3JDaDEgoHUAOZOdwMws+7kBlBiwQw3ADPrbm4AJc5MAfkcgJl1KTeAEgtnNBqAjwDMrDu5AZR4+5KZbF67jLVvvPQL0ZuZVaGVNYEXSXpI0kFJByR9Oo8PSnpY0m5JuyStzuPvlXQyj++WtK3keZdK2inpSUn3SJo0vj/axenvm8iN61YwMKWv6ihmZpdEK0cAp4HPRMQKYA3wSUkrgduAWyJiENiW7zf8PiIG8+3zJc97K3B7RCwHngOuveCfwszMxuxVG0BEHIuIx/L2C8BBYCEQwGV5twHgaKvfVJKAq4B789B3gQ+1HtvMzC7WmK5yJukK4G3ATuA64NeSvkRqJO9q2vWdkvaQmsJnI+LAiKeaDZyIiNP5/hFSUyn6npuBzQCLFy8eS1wzMxtFyyeBJU0DfgpcFxHPA1uA6yNiEXA98O2862PAkoh4K/B14BdFT1cwFkXfNyLujIihiBiaO9cnZM3MxktLDUBSH+nN/4cR8bM8vAlobP8EWA0QEc9HxKm8fT/QJ2nOiKd8BpghqXEEcjljmEIyM7OL18qngET66/5gRHyl6aGjwHvy9lXAk3n/1+evIX8yaALwr+bnjIgAHgI25KFNwC8v/McwM7OxauUcwJXAR4B9knbnsRuBjwNfzX/F/5s8T096U98i6TTwErAxv+Ej6X7gYxFxFPgcsF3SF4A/cnYKyczM2kD5vbkjDA0Nxa5du6qOYWbWUSQ9GhFD5413UgOQ9E/gLxf45XNI5x7qzjnHVyfk7ISM4JzjrZ05l0TEeZ+i6agGcDEk7SrqgHXjnOOrE3J2QkZwzvFWh5y+FpCZWY9yAzAz61G91ADurDpAi5xzfHVCzk7ICM453irP2TPnAMzM7Fy9dARgZmZN3ADMzHpUTzQASVdLOiRpWNLWqvM0SHpK0r7Gojp5bJakHXmhnB2SZlaQ6y5JxyXtbxorzKXka7m2eyWtqjjnzZL+3rQg0bqmx27IOQ9J+mAbc5YtqlSbmo6SsVb1lNQv6RFJe3LOW/J44QJTkibn+8P58Ssqznm3pMNN9RzM49W8jiKiq2/AROBPwDJgErAHWFl1rpztKWDOiLHbgK15eytwawW51gKrgP2vlgtYB/yKdIXXNcDOinPeTLoE+ch9V+b/+8nA0vw7MbFNOecDq/L2dOCJnKc2NR0lY63qmWsyLW/3kS5Nvwb4MemyMwB3AFvy9ieAO/L2RuCeNv2fl+W8G9hQsH8lr6NeOAJYDQxHxJ8j4j/AdmB9xZlGs560QA5UtFBORPwOeHbEcFmu9cD3InmYdJXX+RXmLLMe2B4RL0fEYWCYfAXbSy3KF1WqTU1HyVimknrmmpzKd/vyLShfYKq5xvcC72tcrLKinGUqeR31QgNYCPyt6X7p4jMVCOA3kh5VWvgG4HURcQzSixKYV1m6c5XlqmN9P5UPo+9qmkKrRU6du6hSLWs6IiPUrJ6SJipdmPI4sIN09FG2wNSZnPnxk6QFqdqeMyIa9fxiruftkiaPzJm1pZ690ABaXnymAldGxCrgGtJay2urDnQB6lbfbwJvAAaBY8CX83jlOXX+okqluxaMtSVrQcba1TMi/htpLfLLSUcdK0bJUpuckt4M3AC8CXgHMIt0VeTKcvZCAzgCLGq6X5vFZyJdFpuIOA78nPTL/HTj0C//e7y6hOcoy1Wr+kbE0/mF9z/gW5ydlqg0p4oXVapVTYsy1rWeOdsJ4LekOfOyBabO5MyPD9D6tOF457w6T7VFRLwMfIeK69kLDeAPwPL8KYFJpBNB91WcCUlTJU1vbAMfAPaTsm3Ku9VpoZyyXPcBH82fYlgDnGxMa1RhxLzph0k1hZRzY/5UyFJgOfBImzKVLapUm5qWZaxbPSXNlTQjb08B3k86X1G2wFRzjTcAD0Y+61pBzsebGr5I5yma69n+11E7zjRXfSOdYX+CNFd4U9V5cqZlpE9R7AEONHKR5icfIK2w9gAwq4JsPyId7r9C+svk2rJcpEPXb+Ta7gOGKs75/ZxjL+lFNb9p/5tyzkPANW3M+W7S4fxeYHe+ratTTUfJWKt6Am8hLSC1l/TmuS2PLyM1oGHSErWT83h/vj+cH19Wcc4Hcz33Az/g7CeFKnkd+VIQZmY9qhemgMzMrIAbgJlZj3IDMDPrUW4AZmY9yg3AzKxHuQGYmfUoNwAzsx71fzvqyxYXb+e9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "obs_test = observed.isel(lat = 56, lon  = 56)\n",
    "hist_test = model_hist.isel(lat = 56,lon = 56)\n",
    "test_245 = model_245.isel(lat = 56,lon = 56)\n",
    "\n",
    "#test_245 = bias_correction_model(y_245,q_map)\n",
    "temp = bias_correction(obs_test.values,hist_test.values)\n",
    "q_map = temp['params']\n",
    "#test_245_corr = bias_correction_model(test_245.values,q_map)\n",
    "\n",
    "hist_corr = temp['outputs']\n",
    "\n",
    "# plt.plot(obs_test[0:365])\n",
    "# plt.plot(hist_test[0:365])\n",
    "# plt.plot(hist_corr[0:365])\n",
    "plt.plot(test_245[0:365])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:33:48.627635Z",
     "start_time": "2020-05-16T15:33:38.737080Z"
    }
   },
   "outputs": [],
   "source": [
    "bias_corrected_results_hist = np.zeros([len(model_hist.time.values),\n",
    "                                           len(model_hist.lat.values),\n",
    "                                           len(model_hist.lon.values)])\n",
    "bias_corrected_results_hist[:] = np.nan\n",
    "\n",
    "\n",
    "bias_corrected_results_585 = np.zeros([len(model_585.time.values),\n",
    "                                           len(model_585.lat.values),\n",
    "                                           len(model_585.lon.values)])\n",
    "bias_corrected_results_585[:] = np.nan\n",
    "\n",
    "\n",
    "\n",
    "bias_corrected_results_126 = np.zeros([len(model_126.time.values),\n",
    "                                           len(model_126.lat.values),\n",
    "                                           len(model_126.lon.values)])\n",
    "bias_corrected_results_126[:] = np.nan\n",
    "\n",
    "bias_corrected_results_245 = np.zeros([len(model_245.time.values),\n",
    "                                           len(model_245.lat.values),\n",
    "                                           len(model_245.lon.values)])\n",
    "bias_corrected_results_245[:] = np.nan\n",
    "\n",
    "bias_corrected_results_370 = np.zeros([len(model_370.time.values),\n",
    "                                           len(model_370.lat.values),\n",
    "                                           len(model_370.lon.values)])\n",
    "bias_corrected_results_370[:] = np.nan\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:34:07.284145Z",
     "start_time": "2020-05-16T15:33:48.630100Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(124, 147, 24106)\n"
     ]
    }
   ],
   "source": [
    "model_hist_values = model_hist.values\n",
    "hist_dict = {}\n",
    "hist_dict['time'] = model_hist.time.values\n",
    "hist_dict['lon'] =  model_hist.lon.values\n",
    "hist_dict['lat'] =  model_hist.lat.values\n",
    "\n",
    "model245_values = model_245.values\n",
    "rcp245_dict = {}\n",
    "rcp245_dict['time'] = model_245.time.values\n",
    "rcp245_dict['lon'] =  model_245.lon.values\n",
    "rcp245_dict['lat'] =  model_245.lat.values\n",
    "\n",
    "model585_values = model_585.values\n",
    "\n",
    "rcp585_dict = {}\n",
    "rcp585_dict['time'] = model_585.time.values\n",
    "rcp585_dict['lon'] =  model_585.lon.values\n",
    "rcp585_dict['lat'] =  model_585.lat.values\n",
    "\n",
    "\n",
    "\n",
    "model126_values = model_126.values\n",
    "\n",
    "rcp126_dict = {}\n",
    "rcp126_dict['time'] = model_126.time.values\n",
    "rcp126_dict['lon'] =  model_126.lon.values\n",
    "rcp126_dict['lat'] =  model_126.lat.values\n",
    "\n",
    "model370_values = model_370.values\n",
    "\n",
    "rcp370_dict = {}\n",
    "rcp370_dict['time'] = model_370.time.values\n",
    "rcp370_dict['lon'] =  model_370.lon.values\n",
    "rcp370_dict['lat'] =  model_370.lat.values\n",
    "\n",
    "observation_attr_values = observed.values\n",
    "print(observation_attr_values.shape)\n",
    "#observation_attr_values = observation_attr_values.reshape((-1,len(lats),len(lons)))\n",
    "#print(observation_attr_values.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T15:34:07.289300Z",
     "start_time": "2020-05-16T15:34:07.285866Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(23342, 124, 147)\n",
      "(124, 147, 24106)\n",
      "(31390, 124, 147)\n"
     ]
    }
   ],
   "source": [
    "print(model_hist_values.shape)\n",
    "print(observation_attr_values.shape)\n",
    "print(model126_values.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T19:21:03.467571Z",
     "start_time": "2020-05-16T15:34:07.290775Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.125 81.125\n",
      "7.375 79.875\n",
      "7.375 81.125\n",
      "8.625 77.375\n",
      "8.625 81.125\n",
      "9.875 77.375\n",
      "9.875 78.625\n",
      "11.125 76.125\n",
      "11.125 77.375\n",
      "11.125 78.625\n",
      "12.375 76.125\n",
      "12.375 77.375\n",
      "12.375 78.625\n",
      "12.375 79.875\n",
      "13.625 74.875\n",
      "13.625 76.125\n",
      "13.625 77.375\n",
      "13.625 78.625\n",
      "13.625 79.875\n",
      "14.875 74.875\n",
      "14.875 76.125\n",
      "14.875 77.375\n",
      "14.875 78.625\n",
      "14.875 79.875\n",
      "16.125 73.625\n",
      "16.125 74.875\n",
      "16.125 76.125\n",
      "16.125 77.375\n",
      "16.125 78.625\n",
      "16.125 79.875\n",
      "16.125 81.125\n",
      "17.375 73.625\n",
      "17.375 74.875\n",
      "17.375 76.125\n",
      "17.375 77.375\n",
      "17.375 78.625\n",
      "17.375 79.875\n",
      "17.375 81.125\n",
      "17.375 82.375\n",
      "18.625 73.625\n",
      "18.625 74.875\n",
      "18.625 76.125\n",
      "18.625 77.375\n",
      "18.625 78.625\n",
      "18.625 79.875\n",
      "18.625 81.125\n",
      "18.625 82.375\n",
      "18.625 83.625\n",
      "19.875 73.625\n",
      "19.875 74.875\n",
      "19.875 76.125\n",
      "19.875 77.375\n",
      "19.875 78.625\n",
      "19.875 79.875\n",
      "19.875 81.125\n",
      "19.875 82.375\n",
      "19.875 83.625\n",
      "19.875 84.875\n",
      "19.875 86.125\n",
      "21.125 71.125\n",
      "21.125 73.625\n",
      "21.125 74.875\n",
      "21.125 76.125\n",
      "21.125 77.375\n",
      "21.125 78.625\n",
      "21.125 79.875\n",
      "21.125 81.125\n",
      "21.125 82.375\n",
      "21.125 83.625\n",
      "21.125 84.875\n",
      "21.125 86.125\n",
      "22.375 69.875\n",
      "22.375 71.125\n",
      "22.375 72.375\n",
      "22.375 73.625\n",
      "22.375 74.875\n",
      "22.375 76.125\n",
      "22.375 77.375\n",
      "22.375 78.625\n",
      "22.375 79.875\n",
      "22.375 81.125\n",
      "22.375 82.375\n",
      "22.375 83.625\n",
      "22.375 84.875\n",
      "22.375 86.125\n",
      "22.375 87.375\n",
      "22.375 88.625\n",
      "22.375 89.875\n",
      "22.375 91.125\n",
      "22.375 92.375\n",
      "23.625 68.625\n",
      "23.625 69.875\n",
      "23.625 71.125\n",
      "23.625 72.375\n",
      "23.625 73.625\n",
      "23.625 74.875\n",
      "23.625 76.125\n",
      "23.625 77.375\n",
      "23.625 78.625\n",
      "23.625 79.875\n",
      "23.625 81.125\n",
      "23.625 82.375\n",
      "23.625 83.625\n",
      "23.625 84.875\n",
      "23.625 86.125\n",
      "23.625 87.375\n",
      "23.625 88.625\n",
      "23.625 89.875\n",
      "23.625 91.125\n",
      "23.625 92.375\n",
      "24.875 67.375\n",
      "24.875 68.625\n",
      "24.875 69.875\n",
      "24.875 71.125\n",
      "24.875 72.375\n",
      "24.875 73.625\n",
      "24.875 74.875\n",
      "24.875 76.125\n",
      "24.875 77.375\n",
      "24.875 78.625\n",
      "24.875 79.875\n",
      "24.875 81.125\n",
      "24.875 82.375\n",
      "24.875 83.625\n",
      "24.875 84.875\n",
      "24.875 86.125\n",
      "24.875 87.375\n",
      "24.875 88.625\n",
      "24.875 89.875\n",
      "24.875 91.125\n",
      "24.875 92.375\n",
      "24.875 93.625\n",
      "26.125 62.375\n",
      "26.125 63.625\n",
      "26.125 64.875\n",
      "26.125 66.125\n",
      "26.125 67.375\n",
      "26.125 68.625\n",
      "26.125 69.875\n",
      "26.125 71.125\n",
      "26.125 72.375\n",
      "26.125 73.625\n",
      "26.125 74.875\n",
      "26.125 76.125\n",
      "26.125 77.375\n",
      "26.125 78.625\n",
      "26.125 79.875\n",
      "26.125 81.125\n",
      "26.125 82.375\n",
      "26.125 83.625\n",
      "26.125 84.875\n",
      "26.125 86.125\n",
      "26.125 87.375\n",
      "26.125 88.625\n",
      "26.125 89.875\n",
      "26.125 91.125\n",
      "26.125 92.375\n",
      "26.125 93.625\n",
      "26.125 94.875\n",
      "27.375 63.625\n",
      "27.375 64.875\n",
      "27.375 66.125\n",
      "27.375 67.375\n",
      "27.375 68.625\n",
      "27.375 69.875\n",
      "27.375 71.125\n",
      "27.375 72.375\n",
      "27.375 73.625\n",
      "27.375 74.875\n",
      "27.375 76.125\n",
      "27.375 77.375\n",
      "27.375 78.625\n",
      "27.375 79.875\n",
      "27.375 81.125\n",
      "27.375 82.375\n",
      "27.375 83.625\n",
      "27.375 84.875\n",
      "27.375 86.125\n",
      "27.375 87.375\n",
      "27.375 88.625\n",
      "27.375 89.875\n",
      "27.375 91.125\n",
      "27.375 92.375\n",
      "27.375 93.625\n",
      "27.375 94.875\n",
      "27.375 96.125\n",
      "28.625 62.375\n",
      "28.625 63.625\n",
      "28.625 64.875\n",
      "28.625 66.125\n",
      "28.625 67.375\n",
      "28.625 68.625\n",
      "28.625 69.875\n",
      "28.625 71.125\n",
      "28.625 72.375\n",
      "28.625 73.625\n",
      "28.625 74.875\n",
      "28.625 76.125\n",
      "28.625 77.375\n",
      "28.625 78.625\n",
      "28.625 79.875\n",
      "28.625 81.125\n",
      "28.625 82.375\n",
      "28.625 83.625\n",
      "28.625 84.875\n",
      "28.625 86.125\n",
      "28.625 87.375\n",
      "28.625 88.625\n",
      "28.625 89.875\n",
      "28.625 91.125\n",
      "28.625 92.375\n",
      "28.625 93.625\n",
      "28.625 94.875\n",
      "28.625 96.125\n",
      "28.625 97.375\n",
      "29.875 67.375\n",
      "29.875 68.625\n",
      "29.875 69.875\n",
      "29.875 71.125\n",
      "29.875 72.375\n",
      "29.875 73.625\n",
      "29.875 74.875\n",
      "29.875 76.125\n",
      "29.875 77.375\n",
      "29.875 78.625\n",
      "29.875 79.875\n",
      "29.875 81.125\n",
      "29.875 82.375\n",
      "29.875 83.625\n",
      "29.875 84.875\n",
      "29.875 86.125\n",
      "29.875 87.375\n",
      "29.875 88.625\n",
      "29.875 89.875\n",
      "29.875 91.125\n",
      "29.875 92.375\n",
      "29.875 93.625\n",
      "29.875 94.875\n",
      "29.875 96.125\n",
      "31.125 67.375\n",
      "31.125 68.625\n",
      "31.125 69.875\n",
      "31.125 71.125\n",
      "31.125 72.375\n",
      "31.125 73.625\n",
      "31.125 74.875\n",
      "31.125 76.125\n",
      "31.125 77.375\n",
      "31.125 78.625\n",
      "31.125 79.875\n",
      "31.125 81.125\n",
      "32.375 68.625\n",
      "32.375 69.875\n",
      "32.375 71.125\n",
      "32.375 72.375\n",
      "32.375 73.625\n",
      "32.375 74.875\n",
      "32.375 76.125\n",
      "32.375 77.375\n",
      "32.375 78.625\n",
      "32.375 79.875\n",
      "32.375 81.125\n",
      "33.625 69.875\n",
      "33.625 71.125\n",
      "33.625 72.375\n",
      "33.625 73.625\n",
      "33.625 74.875\n",
      "33.625 76.125\n",
      "33.625 77.375\n",
      "33.625 78.625\n",
      "33.625 79.875\n",
      "33.625 81.125\n",
      "34.875 68.625\n",
      "34.875 69.875\n",
      "34.875 71.125\n",
      "34.875 72.375\n",
      "34.875 73.625\n",
      "34.875 74.875\n",
      "34.875 76.125\n",
      "34.875 77.375\n",
      "34.875 78.625\n",
      "34.875 79.875\n",
      "36.125 72.375\n",
      "36.125 73.625\n",
      "36.125 74.875\n",
      "36.125 76.125\n"
     ]
    }
   ],
   "source": [
    "correct_params = []\n",
    "for i,lat in enumerate(lats):\n",
    "    for j,lon in enumerate(lons):\n",
    "        params_dict = {}\n",
    "\n",
    "        if np.isnan(model_hist_values[0,i,j]) or np.isnan(observation_attr_values[i,j,0]):\n",
    "            bias_corrected_results_hist[:,i,j] = np.nan\n",
    "            params_dict['lat'] = lat\n",
    "            params_dict['lon'] = lon\n",
    "            params_dict['params'] = np.nan\n",
    "        else:\n",
    "            try:\n",
    "\n",
    "                y = model_hist_values[:,i,j]\n",
    "                x = observation_attr_values[i,j,:]\n",
    "                \n",
    "                \n",
    "                y_245 = model245_values[:,i,j]\n",
    "                y_585 = model585_values[:,i,j]\n",
    "                y_126 = model126_values[:,i,j]\n",
    "                y_370 = model370_values[:,i,j]\n",
    "                \n",
    "                temp = bias_correction(x,y)\n",
    "                q_map = temp['params']\n",
    "                temp_245 = bias_correction_model(y_245,q_map)\n",
    "                temp_585 = bias_correction_model(y_585,q_map)\n",
    "                temp_126 = bias_correction_model(y_126,q_map)\n",
    "                temp_370 = bias_correction_model(y_370,q_map)\n",
    "                bias_corrected_results_245[:,i,j] = temp_245['outputs']\n",
    "                bias_corrected_results_585[:,i,j] = temp_585['outputs']\n",
    "                bias_corrected_results_126[:,i,j] = temp_126['outputs']\n",
    "                bias_corrected_results_370[:,i,j] = temp_370['outputs']\n",
    "                bias_corrected_results_hist[:,i,j] = temp['outputs']\n",
    "\n",
    "                if i%5==0 and j%5==0:\n",
    "                    print(lat,lon)\n",
    "                params_dict['lat'] = lat\n",
    "                params_dict['lon'] = lon\n",
    "                params_dict['params'] = temp['params']\n",
    "\n",
    "            except:\n",
    "                bias_corrected_results_hist[:,i,j] = np.nan\n",
    "                bias_corrected_results_245[:,i,j] = np.nan\n",
    "                bias_corrected_results_126[:,i,j] = np.nan\n",
    "                bias_corrected_results_585[:,i,j] = np.nan\n",
    "                bias_corrected_results_370[:,i,j] = np.nan\n",
    "                params_dict['lat'] = lat\n",
    "                params_dict['lon'] = lon\n",
    "                params_dict['params'] = np.nan\n",
    "\n",
    "        correct_params.append(params_dict)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T19:21:03.536588Z",
     "start_time": "2020-05-16T19:21:03.469873Z"
    }
   },
   "outputs": [],
   "source": [
    "ds_hist = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_hist)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':hist_dict['time'] })\n",
    "\n",
    "\n",
    "ds_ssp585 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_585)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp585_dict['time'] })\n",
    "\n",
    "ds_ssp126 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_126)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp126_dict['time'] })\n",
    "\n",
    "ds_ssp245 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_245)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp245_dict['time'] })\n",
    "\n",
    "ds_ssp370 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_370)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp370_dict['time'] })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T19:22:18.526703Z",
     "start_time": "2020-05-16T19:21:03.537979Z"
    }
   },
   "outputs": [],
   "source": [
    "for scenario in scenarios:\n",
    "    name = \"ds_\"+ str(scenario)\n",
    "    ds = globals()[name]\n",
    "    name_netcdf =  model_name_file + \"_\" + variable + \"_\" + str(scenario) + '.nc'\n",
    "    name_pickle =  model_name_file + \"_\" + variable + \"_\" + str(scenario) + '.pkl'\n",
    "        \n",
    "    try:\n",
    "        os.makedirs(target_folder + \"/Bias_corrected_kelvin/\")\n",
    "        netcdf_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_netcdf\n",
    "        pickle_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_pickle\n",
    "        ds.to_netcdf(netcdf_corr)\n",
    "        if scenario == \"hist\":\n",
    "            save_as_pickled_object(correct_params, pickle_corr)\n",
    "\n",
    "\n",
    "    except FileExistsError:\n",
    "        netcdf_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_netcdf\n",
    "        pickle_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_pickle\n",
    "        ds.to_netcdf(netcdf_corr)\n",
    "        if scenario == \"hist\":\n",
    "            save_as_pickled_object(correct_params, pickle_corr)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T19:22:20.144381Z",
     "start_time": "2020-05-16T19:22:18.529520Z"
    }
   },
   "outputs": [],
   "source": [
    "%reset -f"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (master)",
   "language": "python",
   "name": "master"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}