Kraus_Turner_1D_experiment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How the atmosphere mixing the ocean surface boundary layer\n",
    "I. Giddy 2023\n",
    "Adapted from Ilker Fer.\n",
    "\n",
    "More details on the model:\n",
    "Kraus EB, Turner JS (1966) A one-dimensional model of the seasonal thermocline. Tellus 19:88–97\n",
    "\n",
    "If not running on Binder.org, this code is dependent on numpy, matplotlib, and pandas. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This assignment is about the upper ocean, particularly the mixed layer and the pycnocline below, and its response to different atmospheric forcing. In order to simplify the physics and avoid complexity, one-dimensional (1D, vertical) solutions will be obtained and discussed. \n",
    "\n",
    "We will be using a bulk mixed layer model, the Krauss and Turner bulk mixed layer model, which solves for the depth of the mixed layer and its bulk properties (mixed layer average density, temperature, salinity) based on the net surface buoyancy flux and wind-induced mixing. \n",
    "\n",
    "The intricacies of the model are not the focus of this practical, but if you are interested in a more detailed description see Kraus EB, Turner JS (1966) A one-dimensional model of the seasonal thermocline. Tellus 19:88–97. \n",
    "\n",
    "For the sake of this experiment, we will assume that no horizontal processes are at play. \n",
    "\n",
    "Follow the steps below and run the model for a 15 day duration, and discuss the evolution of mixed layer depth and mixed layer density for different cases.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions\n",
    "\n",
    "The below functions are used in the practical. Run them before starting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this makes a generic CTD profile\n",
    "def make_CTDprofile_continuous(NBV=0.015,T0=10,S0=34,H=100):\n",
    "    \"\"\"\n",
    "    %\n",
    "    % profile=make_CTDprofile_continuous(NBV,T0,S0,H)\n",
    "    %\n",
    "    % Generates a CTD profile for given constant buoyancy frequency, NBV (in rad/s),\n",
    "    % starting from top level temperature T0 and salinity S0. \n",
    "    % Output (profile) is consistent with the PWP engine:\n",
    "    %     z: depth +ve downward, (pressure-like)\n",
    "    %     t: temperature\n",
    "    %     s: salinity\n",
    "    %     d: density\n",
    "    % \n",
    "    % Default values:\n",
    "    % NBV=0.015; # rad/s, buoy freq\n",
    "    % T0=10; S0=34; # top level temperature and salinity\n",
    "    % H=100; # total depth\n",
    "    % assumes typical coeff of thermal expansion and haline contraction,\n",
    "    % and a density ratio of -2 \n",
    "    %\n",
    "    % by Ilker Fer, November 2018, for GEOF310 Assignment\n",
    "    % adapted by I. Giddy November 2020 for python\n",
    "    % --------------------\n",
    "    \"\"\"\n",
    "\n",
    "    import numpy as np\n",
    "    import gsw as gsw\n",
    "    import pandas as pd\n",
    "\n",
    "    M=102\n",
    "    dz=H/(M-2)\n",
    "    rho0=gsw.rho(S0,T0,0)\n",
    "    alp=1.5e-4\n",
    "    bet=7.5e-4\n",
    "    Rrho = -2 #density ratio (alpha dT / beta dS)\n",
    "\n",
    "    z=np.arange(0,H,dz)*-1\n",
    "    rho=rho0-(NBV**2/9.81*rho0*z)\n",
    "\n",
    "    T=T0+(-(rho-rho0)/rho0)/(alp*(1-(1/Rrho)))\n",
    "    S=S0 + ((alp*(T-T0)/Rrho)/bet)\n",
    "\n",
    "    d = gsw.rho(S,T,0)  #density referenced to the surface\n",
    "    \n",
    "    # make pandas dataframe \n",
    "    \n",
    "    df = pd.DataFrame({'z':-z,'t':T,'s':S,'d':d})\n",
    "\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Kraus-Turner Bulk Mixed Layer Model\n",
    "def KrausTurner_engine(STRAT,FORCE,duration,PAR):\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "     A function to solve the Krauss Turner bulk layer model based on the \n",
    "     formulations in Olbers, Willebrand and Eden (Ocean Dynamics text book)\n",
    "\n",
    "     Initial stratification, STRAT = \n",
    "                 h0: mixed layer depth\n",
    "             dtop: top layer density\n",
    "             dbot: bottom layer density\n",
    "         TotalDepth: total water depth\n",
    "\n",
    "     Forcing, FORCE = \n",
    "          tau: wind stress, N/m2\n",
    "           FH: net heat flux, W/m2, +ve vertical up, >0 convection; <0 stabilizing\n",
    "        Bpyc: buoyancy flux induced by turbulent mixing in the pycnocline. Turbulent\n",
    "               mixing induces <0, put a negative value or zero\n",
    "     Model PARAMETERS\n",
    "     PAR.r1 = 0.1\n",
    "     PAR.r2 = 0.9\n",
    "     PAR.c = 1\n",
    "\n",
    "     Output returns time evolution of:\n",
    "         time: in days\n",
    "         zML: mixed layer depth\n",
    "          d: mixed layer density\n",
    "\n",
    "     by Ilker Fer, November 2018, for GEOF310 Assignment\n",
    "     adapted by I. Giddy November 2020 for python\n",
    "    \"\"\"\n",
    "    \n",
    "    import numpy as np\n",
    "    \n",
    "    \n",
    "    # Just to convert to formulations in p358 in \n",
    "    # Olbers, Willebrand and Eden (Ocean Dynamics) text book \n",
    "    # Follow notation, \n",
    "    # Buoy flux, letter Q, is positive downward in the ocean\n",
    "    # with Q<0 convection, and Q>0 stable buoy flux (that is turbulence in\n",
    "    # stable stratification induces downward buoy flux\n",
    "    #\n",
    "    #\n",
    "\n",
    "    # prescribed parameters in the model\n",
    "   # r1 = PAR.r1 # fraction of tau dissipated\n",
    "   # r2 = PAR.r2 # fraction of convective Bo dissipated\n",
    "   # c = PAR.c # approx 1, for work done by the wind stress\n",
    "    \n",
    "    r1 = 0.1\n",
    "    r2 = 0.9\n",
    "    c = 1\n",
    "    \n",
    "    # PARAMETERS\n",
    "    rho0=1028 #reference density\n",
    "    rhotop = STRAT.dtop.values[0]\n",
    "    rhobot = STRAT.dbot.values[0]\n",
    "    # Initial mixed layer depth\n",
    "    h0=STRAT.h0.values[0]\n",
    "    D=STRAT.TotalDepth.values[0]\n",
    "\n",
    "    #wind stress in  N/m^2, converted to kinematic:\n",
    "    tau=np.array(FORCE.tau/rho0)\n",
    "    heatflux = np.array(-FORCE.FH) # W/m2, net radiation, converted to Eden's \n",
    "    # heatflux<0 convection; >0 stabilizing\n",
    "    BPYC=np.array(-FORCE.Bpyc) # convert to Eden: buoy flux at the base of ML, turbulent mixing induces downward\n",
    "    # buoy flux. Set to nil or a small +ve (downward) number representative of shear induced\n",
    "    # mixing\n",
    "    \n",
    "    # ------- START --------- #\n",
    "    # convert to buoyancy\n",
    "    b00 = -(9.8/rho0)*(rhotop-rho0)\n",
    "    bst = -(9.8/rho0)*(rhobot-rho0)\n",
    "\n",
    "\n",
    "    dt=900; # time step, s\n",
    "    ti = np.arange(0,duration,dt/(3600*24))\n",
    "    \n",
    "    # In case of constant forcing, generate forcing at all times\n",
    "    # If time variable, interpolate FORCE.time to model time\n",
    "    if 'time' in FORCE.keys():# & (len(tau)>1) :\n",
    "        heatflux = np.interp1(FORCE.time,heatflux,ti)\n",
    "        tau = np.interp1(FORCE.time,tau,ti)\n",
    "        BPYC = np.interp1(FORCE.time,BPYC,ti)\n",
    "    else:\n",
    "        heatflux = np.ones(len(ti))*heatflux\n",
    "        tau = np.ones(len(ti))*tau\n",
    "        BPYC = np.ones(len(ti))*BPYC\n",
    "    ###\n",
    "    \n",
    "    # define arrays to fill\n",
    "    # sign convention for buoy flux, +ve upward\n",
    "    # B>0 convection; B<0 stable\n",
    "\n",
    "    h = [h0]\n",
    "   # h = np.ndarray(len(D))\n",
    "    b0 = [b00]\n",
    "    #h[0]=h0 # thickness of ML   # CHECK THESE\n",
    "    #b0[0]=b00 # buoyancy anomaly in the ML\n",
    "    \n",
    "    \n",
    "\n",
    "    for I in range(1,len(ti)):   #change from 2 to 1 becuase of indexing convention in python\n",
    "\n",
    "        Q0=((9.8*1.5e-4)/(rho0*4000))*heatflux[I-1] # surface buoy flux\n",
    "        Qst = BPYC[I-1] #buoy flux at the base of the pycnocline, turbulent mixing induces <0\n",
    "\n",
    "        windwork = c*(abs(tau[I-1])**(1.5))\n",
    "        if Q0<0: #convection\n",
    "          #  Diss = r1*windwork + r2*(h[I-1]/4)*(abs(Q0)-Q0)    # Dissipation integrated over h\n",
    "            Diss = r1*windwork + r2*(h[I-1]/4)*(abs(Q0)-Q0)    # Dissipation integrated over h\n",
    "        else: #stress only\n",
    "            \n",
    "            Diss = r1*windwork  # Dissipation integrated over h\n",
    "#        \n",
    "    \n",
    "        Qh = (2/h[I-1])*((1-r1)*windwork - (h[I-1]/2)*(((1-(r2/2))*Q0)+(r2/2*abs(Q0))))\n",
    "#\n",
    "        if np.any(Qh) < 0:# stable, Q0>0\n",
    "            h.append(((2*(1-r1)*windwork) / (Qh + (((1-(r2/2))*Q0)+(r2/2*abs(Q0))))))  # work out indexing\n",
    "            dbdt = Q0/h[I-1]\n",
    "            b0.append((b0[I-1]+(dbdt*dt)))\n",
    "#    #         Qh = 0\n",
    "#    #        h(I) = 2*(1-r1)*windwork/Q0; \n",
    "           \n",
    "        else: # deepening\n",
    "            dbdt = (Q0 - Qh)/h[I-1]\n",
    "            b0.append((b0[I-1]+(dbdt*dt)))\n",
    "            dhdt = (Qh - Qst)/(b0[I] - bst)\n",
    "            h.append(((h[I-1]+(dhdt*dt))))\n",
    "        \n",
    "        if np.any(h[I]>D):\n",
    "            h.append(D)\n",
    "\n",
    "    # convert lists to arrays\n",
    "    h=np.array(h)\n",
    "    b0=np.array(b0)\n",
    "            \n",
    "    # convert buoyancy anomaly in the ML to density:\n",
    "    rho =rho0 -(rho0*b0/9.8)\n",
    "#     print(len())\n",
    "    # create a dataframe for outputing\n",
    "    df = pd.DataFrame({'time':ti,\n",
    "                      'zML':h,\n",
    "                      'd':rho})\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a CTD profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a CTD profile\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "df = make_CTDprofile_continuous(NBV=0.015,T0=10,S0=34,H=300)  \n",
    "\n",
    "#let's make a 2-layer profile, 25-m mixed layer\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "pd.options.mode.chained_assignment = None  # default='warn'\n",
    "\n",
    "hML = 25  # initial mixed layer depth\\\n",
    "\n",
    "ind_top = df.index[np.abs(df.z)<hML]  \n",
    "df.t[ind_top]=np.mean(df.t.iloc[ind_top]) \n",
    "df.s[ind_top]=np.mean(df.s.iloc[ind_top])\n",
    "df.d[ind_top]=np.mean(df.d.iloc[ind_top])\n",
    "\n",
    "ind_bot = df.index[np.abs(df.z)>=hML]\n",
    "\n",
    "df.t[ind_bot]=np.mean(df.t.iloc[ind_bot])\n",
    "df.s[ind_bot]=np.mean(df.s.iloc[ind_bot])\n",
    "df.d[ind_bot]=np.mean(df.d.iloc[ind_bot])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Density')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhwAAAHUCAYAAABxrfE4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxL0lEQVR4nO3de1jU1b7H8c8wwACKpKRcUoFdqJmX8nJM09BU1LTy2K7UbmQXO6Y7tpUns53oOUGamZXlSU+hXexmah0rk523CjPzmpejad4qiXYRohIIrPOH2zmMgALOcgDfr+eZZzPrt2b9votfv8fPXrOGcRhjjAAAACzy83UBAACg7iNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAB11Ny5c+VwOPTNN99U6/UOh0MpKSnu59u3b1dKSor27dtXpm9SUpJiY2OrdZ7KvjYpKUkOh8P9cLlcatmypSZOnKg//vijWueuyMaNG5WQkKCwsDA5HA7NmDFDK1eulMPh0MqVK939UlJS5HA4vHpuoK7y93UBAGqmNWvWqGnTpu7n27dv16RJk9SzZ88yAeFvf/ubHnzwQes1BQcHa/ny5ZKknJwcvfXWW5o8ebL+93//V++8847XzjNixAgdPXpUb7/9tho2bKjY2FiFhIRozZo1at26tdfOA5xPCBwAynXllVdWuu/FF19ssZL/5+fn51HXgAEDtG/fPr377ruaPn26LrroonJfl5+fr+Dg4EqfZ+vWrbr33ns1YMAAj/aq/E4AeOItFeA8kpSUpPr162v37t269tprVb9+fTVr1kwPPfSQCgoKPPqWfktl7ty5uummmyRJvXr1cr+tMXfuXPe4p656vPjii7r66qvVpEkT1atXT23bttXUqVN1/Phxr87pZAjYv3+/JCk2NlaDBg3SwoULdcUVVygoKEiTJk2SdCJI3HDDDWrYsKGCgoJ0+eWXa968ee6xTr4NVVRUpFmzZrnnKanct1Qq8s4776hr166qV6+e6tevr379+mnjxo1enTdQ2xA4gPPM8ePHdf3116t379764IMPNGLECD377LOaMmVKha8ZOHCgUlNTJZ0IEmvWrNGaNWs0cODACl+zZ88eDR8+XK+//rqWLFmiu+++W08//bRGjhzp1fns3r1bktS4cWN324YNG/TII4/oL3/5i5YuXaobb7xRO3fuVLdu3bRt2zY9//zzWrhwoVq3bq2kpCRNnTrVPc81a9ZIkv785z+751kVqampGjZsmFq3bq13331Xr7/+uvLy8tSjRw9t377dS7MGah/eUgHOM4WFhZo0aZJ7xaJ379765ptvNH/+fD3xxBPlvqZx48aKj4+XJLVu3bpSby1Mnz7d/XNJSYl69Oih8PBw3XXXXXrmmWfUsGHDatVfVFQkSfr99981f/58LV68WJ07d3bXJ0nZ2dnavn27WrRo4W4bNmyYCgsLtWLFCjVr1kySdO211+r333/XpEmTNHLkSDVu3NgdXCIiIqr8FsrBgwc1ceJEjR49Ws8//7y7vW/fvoqPj9ekSZO8utcEqE1Y4QDOMw6HQ9ddd51HW7t27dxvSXjLxo0bdf311ys8PFxOp1MBAQG64447VFxcrF27dlVrzKNHjyogIEABAQFq3LixkpOTNWDAAC1atMijX7t27TzChiQtX75cvXv3doeNk5KSknTs2LEqr2SU59NPP1VRUZHuuOMOFRUVuR9BQUFKSEio1NsxQF3FCgdwngkJCVFQUJBHm8vl8upHSw8cOKAePXqoZcuWeu655xQbG6ugoCB9/fXXeuCBB5Sfn1+tcYODg7V69Wp3zTExMWrQoEGZflFRUWXafv3113Lbo6Oj3cfP1s8//yxJ6ty5c7nH/fz4/3g4fxE4AHjd4sWLdfToUS1cuFAxMTHu9k2bNp3VuH5+furUqdMZ+5X3tzHCw8N16NChMu0//fSTJOnCCy88q9pKj7FgwQKPeQMgcACoJJfLJUmVWp04+Q/+yddIkjFGc+bMsVNcJfTu3VuLFi3STz/95F7VkKTXXntNISEhXvnIa79+/eTv7689e/boxhtvPOvxgLqEwAGgUtq0aSNJmj17tkJDQxUUFKS4uDiFh4eX6du3b18FBgZq2LBhGjdunP744w/NmjVLOTk557pst4kTJ2rJkiXq1auXnnjiCTVq1EhvvvmmPvroI02dOlVhYWFnfY7Y2FhNnjxZEyZM0Pfff6/+/furYcOG+vnnn/X111+rXr167o/oAucb3lAEUClxcXGaMWOGNm/erJ49e6pz5876n//5n3L7tmrVSu+//75ycnI0ZMgQjRkzRpdffrnHJzfOtZYtWyozM1MtW7bUAw88oMGDB2vr1q1KT0/XI4884rXzjB8/XgsWLNCuXbt05513ql+/fho3bpz279+vq6++2mvnAWobhzHG+LoIAABQt7HCAQAArCNwAAAA6wgcAADAujoTOF566SXFxcUpKChIHTt21Oeff+7rkgAAwD/VicDxzjvvKDk5WRMmTNDGjRvVo0cPDRgwQAcOHPB1aQAAQHXkUypdunRRhw4dNGvWLHfbpZdeqsGDBystLc2HlQEAAKkO/OGvwsJCrV+/Xo8++qhHe2JiojIzM8t9TUFBgQoKCtzPS0pK9Ntvvyk8PLzcP4kMAADKZ4xRXl6eoqOjT/t9QbU+cPzjH/9QcXGxIiIiPNojIiKUlZVV7mvS0tL4a38AAHjRwYMH1bRp0wqP1/rAcdKpKxPGmApXK8aPH6+xY8e6n+fm5qp58+Y6ePBgud88WVUfbflJ//7+t5KkRwe01G1Xxp71mAAA1ESHDx9Ws2bNFBoaetp+tT5wXHjhhXI6nWVWM7Kzs8usepzkcrk8vlTqpAYNGnglcITUz5OfK0SSFFwv1CtjAgBQk51pS0Kt/5RKYGCgOnbsqIyMDI/2jIwMdevWzUdVAQCA0mr9CockjR07Vrfffrs6deqkrl27avbs2Tpw4IDuv/9+X5cGAABURwLHLbfcol9//VWTJ0/WoUOH1KZNG3388ceKiYnxdWkAAEB1JHBI0qhRozRq1ChflwEAAMpR6/dwAACAmo/AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsq9GBIy0tTZ07d1ZoaKiaNGmiwYMHa+fOnR59jDFKSUlRdHS0goOD1bNnT23bts1HFQMAgPLU6MCxatUqPfDAA/rqq6+UkZGhoqIiJSYm6ujRo+4+U6dO1fTp0zVz5kytW7dOkZGR6tu3r/Ly8nxYOQAAKM3f1wWcztKlSz2ep6enq0mTJlq/fr2uvvpqGWM0Y8YMTZgwQUOGDJEkzZs3TxEREZo/f75Gjhzpi7IBAMApavQKx6lyc3MlSY0aNZIk7d27V1lZWUpMTHT3cblcSkhIUGZmZoXjFBQU6PDhwx4PAABgT60JHMYYjR07Vt27d1ebNm0kSVlZWZKkiIgIj74RERHuY+VJS0tTWFiY+9GsWTN7hQMAgNoTOEaPHq0tW7borbfeKnPM4XB4PDfGlGkrbfz48crNzXU/Dh486PV6AQDA/6vRezhOGjNmjD788EOtXr1aTZs2dbdHRkZKOrHSERUV5W7Pzs4us+pRmsvlksvlslcwAADwUKNXOIwxGj16tBYuXKjly5crLi7O43hcXJwiIyOVkZHhbissLNSqVavUrVu3c10uAACoQI1e4XjggQc0f/58ffDBBwoNDXXvywgLC1NwcLAcDoeSk5OVmpqq+Ph4xcfHKzU1VSEhIRo+fLiPqwcAACfV6MAxa9YsSVLPnj092tPT05WUlCRJGjdunPLz8zVq1Cjl5OSoS5cuWrZsmUJDQ89xtQAAoCI1OnAYY87Yx+FwKCUlRSkpKfYLAgAA1VKj93AAAIC6gcABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKyr0R+LrQsO/pavzQd/V6C/34mH00+ukz//87m/k9wHAKjbCByWvfrlXr365d7T9vFzyB0+Av2d/x9InJ7BJLBUUHGd5ljZUOMs07fcc5R67em+/A4AgKoicFhwceP6VepfYqQ/jpfoj+MlkorsFFVFAU5HuWEk0N9ZtcBT7hilj1cuDPn7OQhBAFCLETgsaHNRmF5N6qRv9uWosKhEhcUlJ/63qEQFpX4+9djJnwuKSlRYVOx+XnLmP7jqdceLjY4XF+toYfG5P3k5HA55BpIKgkygv7OKgaf84OQ6NQyd+lqnn/z8CEAAUFkEDkuuaRWha1pFeGWsouJTQksFQaXcY6WCS2UCT8EZwtDJn881Y6SCf84v75yfvXz+fo7yV3nK2adTOgxVPvBULQwFOFkFAlBzEThqAf9/biwNCfR1JScYY6oQeM60wlNc5cBTXuAq9sEyUFGJUVFhsY7VkFUgSWXf7qoguJQfeJyVCjxnOoer1DhOVoEA/BOBA1XmcDjk8nfK5e/0dSluxSWmVKgpLj+0VHmFp+w4BWcIQ6XP4Qsna1CBT05fhtPPUcnVnzOt8DgrDjxn3DTt+TOrQIBvEDhQJzj9HAoOdCo40CkpwNflyBij48WmnNWe4tOu4pxphccz8BSfdvXn+CmvP1587leBikuM8kuKlX+8Bq0CVRRWKvmWVvmBx3nGwBNQwd4iNkTjfEHgACxwOBwK9D+xx0MuX1dzQknJPwNQBas85a/eFJcbZir3dllF5/DthuiTv4OasgrEhmicLwgcwHnCz8+hID+nggJqzlthbIhmQzQbos8fBA4APsOG6DMHLjZEn1gFCnCyIbq2I3AAwD+xIfo0qz8+3BBtDBui68KGaAIHANRgtXFD9PFi4xFsyg0zp4adSq0eFZfbjw3RJ1RmQ3RQgFPXtYvSTZ2anfP6CBwAgEpjQ3T5q0C1aUP0l7v/oWtaNVF4/XN7AQkcAIBajQ3RZ94fVHpDdHGJUc6x4wQOAABqu5q4Ifqh9zZr4YYffVaDn8/ODAAAzgmHwyE/H28qJXAAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsK5WBY60tDQ5HA4lJye724wxSklJUXR0tIKDg9WzZ09t27bNd0UCAIAyak3gWLdunWbPnq127dp5tE+dOlXTp0/XzJkztW7dOkVGRqpv377Ky8vzUaUAAOBUtSJwHDlyRLfeeqvmzJmjhg0butuNMZoxY4YmTJigIUOGqE2bNpo3b56OHTum+fPn+7BiAABQWq0IHA888IAGDhyoPn36eLTv3btXWVlZSkxMdLe5XC4lJCQoMzOzwvEKCgp0+PBhjwcAALDH39cFnMnbb7+tDRs2aN26dWWOZWVlSZIiIiI82iMiIrR///4Kx0xLS9OkSZO8WygAAKhQjV7hOHjwoB588EG98cYbCgoKqrCfw+HweG6MKdNW2vjx45Wbm+t+HDx40Gs1AwCAsmr0Csf69euVnZ2tjh07utuKi4u1evVqzZw5Uzt37pR0YqUjKirK3Sc7O7vMqkdpLpdLLpfLXuEAAMBDjV7h6N27t7799ltt2rTJ/ejUqZNuvfVWbdq0SX/6058UGRmpjIwM92sKCwu1atUqdevWzYeVAwCA0mr0CkdoaKjatGnj0VavXj2Fh4e725OTk5Wamqr4+HjFx8crNTVVISEhGj58uC9KBgAA5ajRgaMyxo0bp/z8fI0aNUo5OTnq0qWLli1bptDQUF+XBgAA/qnWBY6VK1d6PHc4HEpJSVFKSopP6gEAAGdWo/dwAACAuoHAAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArKv219OXlJRo9+7dys7OVklJicexq6+++qwLAwAAdUe1AsdXX32l4cOHa//+/TLGeBxzOBwqLi72SnEAAKBuqFbguP/++9WpUyd99NFHioqKksPh8HZdAACgDqlW4Pjuu++0YMECXXLJJd6uBwAA1EHV2jTapUsX7d6929u1AACAOqrSKxxbtmxx/zxmzBg99NBDysrKUtu2bRUQEODRt127dt6rEAAA1HqVDhyXX365HA6HxybRESNGuH8+eYxNowAA4FSVDhx79+61WQcAAKjDKh04YmJi3D+vXr1a3bp1k7+/58uLioqUmZnp0RcAAKBam0Z79eql3377rUx7bm6uevXqddZFAQCAuqVagePkXo1T/frrr6pXr95ZFwUAAOqWKv0djiFDhkg6sUE0KSlJLpfLfay4uFhbtmxRt27dvFshAACo9aoUOMLCwiSdWOEIDQ1VcHCw+1hgYKCuvPJK3Xvvvd6tEAAA1HpVChzp6emSpNjYWD388MO8fQIAACqlWn/afOLEiZKk7Oxs7dy5Uw6HQy1atFCTJk28WhwAAKgbqrVp9PDhw7r99tt10UUXKSEhQVdffbUuuugi3XbbbcrNzfV2jQAAoJarVuC45557tHbtWi1ZskS///67cnNztWTJEn3zzTfs4QAAAGVU6y2Vjz76SJ9++qm6d+/ubuvXr5/mzJmj/v37e604AABQN1RrhSM8PNz9iZXSwsLC1LBhw7MuCgAA1C3VChyPP/64xo4dq0OHDrnbsrKy9Mgjj+hvf/ub14oDAAB1Q7XeUpk1a5Z2796tmJgYNW/eXJJ04MABuVwu/fLLL3r55ZfdfTds2OCdSgEAQK1VrcAxePBgL5cBAADqsrP6OxwAAACVUa09HJL0+++/67//+781fvx49zfHbtiwQT/++KPXigMAAHVDtVY4tmzZoj59+igsLEz79u3Tvffeq0aNGmnRokXav3+/XnvtNW/XCQAAarFqrXCMHTtWSUlJ+u677xQUFORuHzBggFavXu214gAAQN1QrcCxbt06jRw5skz7RRddpKysrLMuCgAA1C3VChxBQUE6fPhwmfadO3eqcePGZ10UAACoW6oVOG644QZNnjxZx48flyQ5HA4dOHBAjz76qG688UavFggAAGq/agWOadOm6ZdfflGTJk2Un5+vhIQEXXLJJQoNDdWTTz7p7RoBAEAtV61PqTRo0EBffPGFVqxYofXr16ukpEQdOnRQnz59vF0fAACoA6ocOEpKSjR37lwtXLhQ+/btk8PhUFxcnCIjI2WMkcPhsFEnAACoxar0looxRtdff73uuece/fjjj2rbtq0uu+wy7d+/X0lJSfrXf/1XW3UCAIBarEorHHPnztXq1av12WefqVevXh7Hli9frsGDB+u1117THXfc4dUiAQBA7ValFY633npLjz32WJmwIUnXXHONHn30Ub355pteKw4AANQNVQocW7ZsUf/+/Ss8PmDAAG3evPmsiwIAAHVLlQLHb7/9poiIiAqPR0REKCcn56yLAgAAdUuVAkdxcbH8/Sve9uF0OlVUVHTWRQEAgLqlSptGjTFKSkqSy+Uq93hBQYFXigIAAHVLlQLHnXfeecY+fEIFAACcqkqBIz093VYdAACgDqvWd6kAAABUBYEDAABYR+AAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANYROAAAgHUEDgAAYB2BAwAAWEfgAAAA1tX4wPHjjz/qtttuU3h4uEJCQnT55Zdr/fr17uPGGKWkpCg6OlrBwcHq2bOntm3b5sOKAQDAqWp04MjJydFVV12lgIAAffLJJ9q+fbueeeYZXXDBBe4+U6dO1fTp0zVz5kytW7dOkZGR6tu3r/Ly8nxXOAAA8FClb4s916ZMmaJmzZp5fEttbGys+2djjGbMmKEJEyZoyJAhkqR58+YpIiJC8+fP18iRI891yQAAoBw1eoXjww8/VKdOnXTTTTepSZMmuuKKKzRnzhz38b179yorK0uJiYnuNpfLpYSEBGVmZlY4bkFBgQ4fPuzxAAAA9tTowPH9999r1qxZio+P16effqr7779ff/nLX/Taa69JkrKysiRJERERHq+LiIhwHytPWlqawsLC3I9mzZrZmwQAAKjZgaOkpEQdOnRQamqqrrjiCo0cOVL33nuvZs2a5dHP4XB4PDfGlGkrbfz48crNzXU/Dh48aKV+AABwQo0OHFFRUWrdurVH26WXXqoDBw5IkiIjIyWpzGpGdnZ2mVWP0lwulxo0aODxAAAA9tTowHHVVVdp586dHm27du1STEyMJCkuLk6RkZHKyMhwHy8sLNSqVavUrVu3c1orAACoWI3+lMpf//pXdevWTampqbr55pv19ddfa/bs2Zo9e7akE2+lJCcnKzU1VfHx8YqPj1dqaqpCQkI0fPhwH1cPAABOqtGBo3Pnzlq0aJHGjx+vyZMnKy4uTjNmzNCtt97q7jNu3Djl5+dr1KhRysnJUZcuXbRs2TKFhob6sHIAAFCawxhjfF2Erx0+fFhhYWHKzc1lPwcAoE56+L3NWrD+B0nS38cm6JIm9b0ybmX/Da3RezgAAEDdQOAAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANYROAAAgHUEDgAAYB2BAwAAWEfgAAAA1hE4AACAdQQOAABgHYEDAABYR+AAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANYROAAAgHUEDgAAYB2BAwAAWEfgAAAA1hE4AACAdQQOAABgHYEDAABYR+AAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANYROAAAgHUEDgAAYB2BAwAAWEfgAAAA1hE4AACAdQQOAABgHYEDAABYR+AAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANYROAAAgHUEDgAAYB2BAwAAWEfgAAAA1hE4AACAdQQOAABgHYEDAABYR+AAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANYROAAAgHUEDgAAYB2BAwAAWEfgAAAA1hE4AACAdQQOAABgHYEDAABYR+AAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANbV6MBRVFSkxx9/XHFxcQoODtaf/vQnTZ48WSUlJe4+xhilpKQoOjpawcHB6tmzp7Zt2+bDqgEAwKlqdOCYMmWK/uu//kszZ87Ujh07NHXqVD399NN64YUX3H2mTp2q6dOna+bMmVq3bp0iIyPVt29f5eXl+bByAABQWo0OHGvWrNENN9yggQMHKjY2Vn/+85+VmJiob775RtKJ1Y0ZM2ZowoQJGjJkiNq0aaN58+bp2LFjmj9/vo+rBwAAJ9XowNG9e3d99tln2rVrlyRp8+bN+uKLL3TttddKkvbu3ausrCwlJia6X+NyuZSQkKDMzMwKxy0oKNDhw4c9HgAAwB5/XxdwOv/+7/+u3NxctWrVSk6nU8XFxXryySc1bNgwSVJWVpYkKSIiwuN1ERER2r9/f4XjpqWladKkSfYKBwAAHmr0Csc777yjN954Q/Pnz9eGDRs0b948TZs2TfPmzfPo53A4PJ4bY8q0lTZ+/Hjl5ua6HwcPHrRSPwAAOKFGr3A88sgjevTRRzV06FBJUtu2bbV//36lpaXpzjvvVGRkpKQTKx1RUVHu12VnZ5dZ9SjN5XLJ5XLZLR4AALjV6BWOY8eOyc/Ps0Sn0+n+WGxcXJwiIyOVkZHhPl5YWKhVq1apW7du57RWAABQsRq9wnHdddfpySefVPPmzXXZZZdp48aNmj59ukaMGCHpxFspycnJSk1NVXx8vOLj45WamqqQkBANHz7cx9UDAICTanTgeOGFF/S3v/1No0aNUnZ2tqKjozVy5Eg98cQT7j7jxo1Tfn6+Ro0apZycHHXp0kXLli1TaGioDysHAAClOYwxxtdF+Nrhw4cVFham3NxcNWjQwNflAADgdQ+/t1kL1v8gSfr72ARd0qS+V8at7L+hNXoPBwAAqBsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOgIHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwDoCBwAAsI7AAQAArCNwAAAA6wgcAADAOp8GjtWrV+u6665TdHS0HA6HFi9e7HHcGKOUlBRFR0crODhYPXv21LZt2zz6FBQUaMyYMbrwwgtVr149XX/99frhhx/O4SwAAMCZ+DRwHD16VO3bt9fMmTPLPT516lRNnz5dM2fO1Lp16xQZGam+ffsqLy/P3Sc5OVmLFi3S22+/rS+++EJHjhzRoEGDVFxcfK6mAQAAzsDflycfMGCABgwYUO4xY4xmzJihCRMmaMiQIZKkefPmKSIiQvPnz9fIkSOVm5urV155Ra+//rr69OkjSXrjjTfUrFkz/f3vf1e/fv3O2VwAAEDFfBo4Tmfv3r3KyspSYmKiu83lcikhIUGZmZkaOXKk1q9fr+PHj3v0iY6OVps2bZSZmVlh4CgoKFBBQYH7eW5uriTp8OHDlmYDAIBvFRw7opKCY5KkvLzDOhxU4pVxT/7baYw5bb8aGziysrIkSRERER7tERER2r9/v7tPYGCgGjZsWKbPydeXJy0tTZMmTSrT3qxZs7MtGwCAGq/DDO+PmZeXp7CwsAqP19jAcZLD4fB4bowp03aqM/UZP368xo4d635eUlKi3377TeHh4Wccu644fPiwmjVrpoMHD6pBgwa+LuecYu7n39zP13lLzJ2525+7MUZ5eXmKjo4+bb8aGzgiIyMlnVjFiIqKcrdnZ2e7Vz0iIyNVWFionJwcj1WO7OxsdevWrcKxXS6XXC6XR9sFF1zgxeprjwYNGpx3N+JJzP38m/v5Om+JuTN3u063snFSjf07HHFxcYqMjFRGRoa7rbCwUKtWrXKHiY4dOyogIMCjz6FDh7R169bTBg4AAHBu+XSF48iRI9q9e7f7+d69e7Vp0yY1atRIzZs3V3JyslJTUxUfH6/4+HilpqYqJCREw4cPl3QiUd1999166KGHFB4erkaNGunhhx9W27Zt3Z9aAQAAvufTwPHNN9+oV69e7ucn91Xceeedmjt3rsaNG6f8/HyNGjVKOTk56tKli5YtW6bQ0FD3a5599ln5+/vr5ptvVn5+vnr37q25c+fK6XSe8/nUJi6XSxMnTizz1tL5gLmff3M/X+ctMXfmXnPm7jBn+hwLAADAWaqxezgAAEDdQeAAAADWETgAAIB1BA4AAGAdgaOWSUtLU+fOnRUaGqomTZpo8ODB2rlzZ5l+O3bs0PXXX6+wsDCFhobqyiuv1IEDByocd+7cuXI4HGUef/zxh0e/l156SXFxcQoKClLHjh31+eefe32OFbE19549e5Y794EDB7r7pKSklDl+8o/TnQuVmXt5c3A4HHr66adPO/b777+v1q1by+VyqXXr1lq0aFGZPr667rbmPWfOHPXo0UMNGzZUw4YN1adPH3399dceferqNa8r93p15l5X7vUjR45o9OjRatq0qYKDg3XppZdq1qxZZxzb5/e6Qa3Sr18/k56ebrZu3Wo2bdpkBg4caJo3b26OHDni7rN7927TqFEj88gjj5gNGzaYPXv2mCVLlpiff/65wnHT09NNgwYNzKFDhzwepb399tsmICDAzJkzx2zfvt08+OCDpl69emb//v3W5luarbn/+uuvHnPeunWrcTqdJj093d1n4sSJ5rLLLvPol52dbXO6Hioz91Ov3auvvmocDofZs2dPheNmZmYap9NpUlNTzY4dO0xqaqrx9/c3X331lbuPL6+7rXkPHz7cvPjii2bjxo1mx44d5q677jJhYWHmhx9+cPepq9e8rtzr1Zl7XbnX77nnHnPxxRebFStWmL1795qXX37ZOJ1Os3jx4grHrQn3OoGjlsvOzjaSzKpVq9xtt9xyi7ntttuqNE56eroJCws7bZ9/+Zd/Mffff79HW6tWrcyjjz5apXN5i7fmfqpnn33WhIaGetzgEydONO3btz+rcb2pvLmf6oYbbjDXXHPNace5+eabTf/+/T3a+vXrZ4YOHep+XpOuu7fmfaqioiITGhpq5s2b526rq9e8rtzrp6rOda+t9/pll11mJk+e7NGvQ4cO5vHHH69wnJpwr/OWSi2Xm5srSWrUqJGkE19E99FHH6lFixbq16+fmjRpoi5dumjx4sVnHOvIkSOKiYlR06ZNNWjQIG3cuNF9rLCwUOvXr1diYqLHaxITE5WZmem9CVWBN+de2iuvvKKhQ4eqXr16Hu3fffedoqOjFRcXp6FDh+r777/3yjyq49S5n+rnn3/WRx99pLvvvvu046xZs6bMNe3Xr5/7mta06+6teZ/q2LFjOn78eJlx6+I1l2r/vX6q6l732nqvd+/eXR9++KF+/PFHGWO0YsUK7dq1S/369atwnBpxr3sltsAnSkpKzHXXXWe6d+/ubjt06JCRZEJCQsz06dPNxo0bTVpamnE4HGblypUVjrVmzRrz+uuvm02bNpnVq1ebG2+80QQHB5tdu3YZY4z58ccfjSTz5ZdferzuySefNC1atLAzwdPw5txLW7t2rZFk1q5d69H+8ccfmwULFpgtW7aYjIwMk5CQYCIiIsw//vEPr86rMsqb+6mmTJliGjZsaPLz8087VkBAgHnzzTc92t58800TGBhojKlZ192b8z7VqFGjzMUXX+zxurp6zevCvX6q6lz32nyvFxQUmDvuuMNIMv7+/iYwMNC89tprpx2rJtzrBI5abNSoUSYmJsYcPHjQ3XbyP5phw4Z59L3uuus8ls7OpLi42LRv396MGTPGY9zMzEyPfv/5n/9pWrZseRazqB5bc7/vvvtMmzZtztjvyJEjJiIiwjzzzDNVK9wLypv7qVq2bGlGjx59xrECAgLM/PnzPdreeOMN43K5jDE167p7c96lnfzHavPmzaftV1eu+alq471+qurMvTbf608//bRp0aKF+fDDD83mzZvNCy+8YOrXr28yMjIqHKsm3Os19uvpcXpjxozRhx9+qNWrV6tp06bu9gsvvFD+/v5q3bq1R/9LL71UX3zxRaXH9/PzU+fOnfXdd9+5x3U6ncrKyvLol52drYiIiLOYSdXZmvuxY8f09ttva/LkyWfsW69ePbVt29b9+zlXKpp7aZ9//rl27typd95554zjRUZGnvaa1pTr7u15nzRt2jSlpqbq73//u9q1a3favnXlmp+qNt7rpVVn7rX5Xs/Pz9djjz2mRYsWuT9d065dO23atEnTpk2r8ItLa8K9zh6OWsYYo9GjR2vhwoVavny54uLiPI4HBgaqc+fOZT5GtWvXLsXExFTpPJs2bVJUVJR73I4dOyojI8OjX0ZGhrp161bN2VSN7bm/++67Kigo0G233XbGvgUFBdqxY4f792PbmeZe2iuvvKKOHTuqffv2Zxy3a9euZa7psmXL3NfU19fd1rwl6emnn9Z//Md/aOnSperUqdMZ+9eVa17eeWrbvV5adeZem+/148eP6/jx4/Lz8/zn2+l0qqSkpMJxa8S97pV1Epwz//Zv/2bCwsLMypUrPT62dezYMXefhQsXmoCAADN79mzz3XffmRdeeME4nU7z+eefu/vcfvvtHjuPU1JSzNKlS82ePXvMxo0bzV133WX8/f093t88+ZGpV155xWzfvt0kJyebevXqmX379tXquZ/UvXt3c8stt5R77oceesisXLnSfP/99+arr74ygwYNMqGhoTVq7sYYk5uba0JCQsysWbPKHefUuX/55ZfG6XSap556yuzYscM89dRTFX5UzhfX3da8p0yZYgIDA82CBQs8xs3Ly3P3qavXvK7c68ZUfe4n1fZ7PSEhwVx22WVmxYoV5vvvvzfp6ekmKCjIvPTSS+4+NfFeJ3DUMpLKfZT+HLkxxrzyyivmkksuMUFBQaZ9+/ZlPp+dkJBg7rzzTvfz5ORk07x5cxMYGGgaN25sEhMTy7yXZ4wxL774oomJiTGBgYGmQ4cOp/2YmrfZmrsxxuzcudNIMsuWLSv33LfccouJiooyAQEBJjo62gwZMsRs27bNm9M7rcrO/eWXXzbBwcHm999/L3ec8ub+3nvvmZYtW5qAgADTqlUr8/7775d5na+uu615x8TElDvuxIkT3X3q6jWvS/d6df57rwv3+qFDh0xSUpKJjo42QUFBpmXLluaZZ54xJSUl7j418V7n6+kBAIB17OEAAADWETgAAIB1BA4AAGAdgQMAAFhH4AAAANYROAAAgHUEDgAAYB2BAwAAWEfgAFAnJCUlafDgwb4uA0AFCBwAzkpSUpIcDoccDocCAgIUERGhvn376tVXXz3tl0l523PPPae5c+e6n/fs2VPJycnn7PwATo/AAeCs9e/fX4cOHdK+ffv0ySefqFevXnrwwQc1aNAgFRUVnZMawsLCdMEFF5yTcwGoOgIHgLPmcrkUGRmpiy66SB06dNBjjz2mDz74QJ988ol71SE3N1f33XefmjRpogYNGuiaa67R5s2b3WOkpKTo8ssv1+uvv67Y2FiFhYVp6NChysvLc/dZsGCB2rZtq+DgYIWHh6tPnz46evSoJM+3VJKSkrRq1So999xz7tWXvXv36pJLLtG0adM8at+6dav8/Py0Z88eu78k4DxH4ABgxTXXXKP27dtr4cKFMsZo4MCBysrK0scff6z169erQ4cO6t27t3777Tf3a/bs2aPFixdryZIlWrJkiVatWqWnnnpKknTo0CENGzZMI0aM0I4dO7Ry5UoNGTJE5X3/5HPPPaeuXbvq3nvv1aFDh3To0CE1b95cI0aMUHp6ukffV199VT169NDFF19s9xcCnOcIHACsadWqlfbt26cVK1bo22+/1XvvvadOnTopPj5e06ZN0wUXXKAFCxa4+5eUlGju3Llq06aNevToodtvv12fffaZpBOBo6ioSEOGDFFsbKzatm2rUaNGqX79+mXOGxYWpsDAQIWEhCgyMlKRkZFyOp266667tHPnTn399deSpOPHj+uNN97QiBEjzs0vBDiPETgAWGOMkcPh0Pr163XkyBGFh4erfv367sfevXs93sqIjY1VaGio+3lUVJSys7MlSe3bt1fv3r3Vtm1b3XTTTZozZ45ycnKqVE9UVJQGDhyoV199VZK0ZMkS/fHHH7rpppu8MFsAp+Pv6wIA1F07duxQXFycSkpKFBUVpZUrV5bpU3qjZ0BAgMcxh8Ph/qSL0+lURkaGMjMztWzZMr3wwguaMGGC1q5dq7i4uErXdM899+j222/Xs88+q/T0dN1yyy0KCQmp1vwAVB6BA4AVy5cv17fffqu//vWvatq0qbKysuTv76/Y2Nhqj+lwOHTVVVfpqquu0hNPPKGYmBgtWrRIY8eOLdM3MDBQxcXFZdqvvfZa1atXT7NmzdInn3yi1atXV7seAJVH4ABw1goKCpSVlaXi4mL9/PPPWrp0qdLS0jRo0CDdcccd8vPzU9euXTV48GBNmTJFLVu21E8//aSPP/5YgwcPVqdOnc54jrVr1+qzzz5TYmKimjRporVr1+qXX37RpZdeWm7/2NhYrV27Vvv27VP9+vXVqFEj+fn5yel0KikpSePHj9cll1yirl27evvXAaAc7OEAcNaWLl2qqKgoxcbGqn///lqxYoWef/55ffDBB3I6nXI4HPr444919dVXa8SIEWrRooWGDh2qffv2KSIiolLnaNCggVavXq1rr71WLVq00OOPP65nnnlGAwYMKLf/ww8/LKfTqdatW6tx48Y6cOCA+9jdd9+twsJCNosC55DDlPeZMgCow7788kv17NlTP/zwQ6UDD4CzQ+AAcN4oKCjQwYMHdd999ykqKkpvvvmmr0sCzhu8pQLgvPHWW2+pZcuWys3N1dSpU31dDnBeYYUDAABYxwoHAACwjsABAACsI3AAAADrCBwAAMA6AgcAALCOwAEAAKwjcAAAAOsIHAAAwLr/A55YEdhpGttWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualise your initial profile\n",
    "\n",
    "plt.figure(figsize=[6,5])\n",
    "plt.title('Initial Profile')\n",
    "\n",
    "plt.plot(df.d-1000,df.z,lw='2')\n",
    "plt.ylim(100,0)\n",
    "plt.ylabel('Depth')\n",
    "plt.xlabel('Density')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 1: Neutral conditions: Buoyancy flux is 0 and wind stress is weak\n",
    "Set tau to 0.05 and FH to 0\n",
    "## Case 2: Stabilising conditions: Buoyancy flux is negative (e.g. ocean is warming) \n",
    "### 2a wind stress is weak\n",
    "Set tau to 0.05 and FH to -500\n",
    "### 2b wind stress is strong\n",
    "Set tau to 0.2 and FH to -500\n",
    "\n",
    "## Case 3: Convection: Buoyancy flux is positive (ocean is cooling) \n",
    "### 3a wind stress is weak\n",
    "Set tau to 0.05 and FH to 500\n",
    "### 3b wind stress is strong\n",
    "Set tau to 0.2 and FH to 500\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#######\n",
    "### Edit this block for each scenario\n",
    "######\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "FORCE = pd.DataFrame({'FH':500,    #<<< EDIT THIS heat flux - note positive out of the ocean (cooling)\n",
    "                     'tau':0.05,    #<<<< EDIT THIS wind stress\n",
    "                     'fcorio':1e-4},index=[0])\n",
    "\n",
    "# For Kraus and Turner\n",
    "\n",
    "FORCE['Bpyc'] = -1e-8 #  W/kg\n",
    "\n",
    "STRAT = df\n",
    "STRAT['h0'] = hML # initial mixed layer depth\n",
    "STRAT['dtop'] = df.d.iloc[0] # surface density\n",
    "STRAT['dbot'] = df.d.iloc[-1] # density at the bottom of the profile\n",
    "STRAT['TotalDepth'] = np.nanmax(np.abs(df.z))\n",
    "\n",
    "PAR = pd.DataFrame({'r1': 0.1,\n",
    "                    'r2': 0.9,\n",
    "                     'c':1},index=[0])\n",
    "\n",
    "duration = 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "kt=KrausTurner_engine(STRAT,FORCE,duration,PAR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'density (kg/m3)')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/MAAAH/CAYAAAAboY3xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACVs0lEQVR4nOzdd3QUZcPG4d+mF5JACIEAoZfQawggILwqYANE6b2KFPG1o/hZX1ERGyq9SUdEBVEUFZAeeqiB0CGEEALpfff7I7CKAiaQZLKb+zpnz3FmJ8OdmOzsvTPzPCaLxWJBRERERERERGyGg9EBRERERERERCR3VOZFREREREREbIzKvIiIiIiIiIiNUZkXERERERERsTEq8yIiIiIiIiI2RmVeRERERERExMaozIuIiIiIiIjYGCejAxQ2ZrOZyMhIvLy8MJlMRscRERH5VxaLhYSEBMqWLYuDgz6nv07HdBERsUU5Pa6rzP9NZGQkgYGBRscQERHJtbNnz1K+fHmjYxQaOqaLiIgt+7fjusr833h5eQHZPzhvb2+D04iIiPy7+Ph4AgMDrccwyaZjuoiI2KKcHtdV5v/m+mV43t7eOvCLiIhN0aXkN9IxXUREbNm/Hdd1Y52IiIiIiIiIjVGZFxEREREREbExKvMiIiIiIiIiNkZlXkRERERERMTGqMyLiIiIiIiI2BiVeREREREREREbozIvIiIiIiIiYmNU5kVERERERERsjMq8iIiIiIiIiI1RmRcRERERERGxMSrzIiIiIiIiIjZGZV5ERERERETExqjMi4iIiIiIiNgYlXkRERERERERG2MTZf7UqVMMGTKEypUr4+7uTtWqVXn99ddJT0+/YTuTyfSPx9SpUw1KLSIiIiIiIpI/nIwOkBNHjhzBbDYzbdo0qlWrxoEDBxg2bBhJSUl8+OGHN2w7Z84cOnbsaF328fEp6LhWKelZODuacHK0ic9MREREREREJJeyzBaS0zPxcnMu0H/XJsp8x44dbyjoVapUITw8nClTpvyjzBcvXpwyZcoUdMSbeuuHgxy6kMCkbvWp5u9ldBwRERERERHJQ+mZZv67dC+RcSksGBKCp2vBVWybPWUcFxeHr6/vP9aPHj0aPz8/goODmTp1Kmaz+bb7SUtLIz4+/oZHXvjj6CUWh55l39mrPPTZJqZtOE6W2ZIn+xYRERERERFjpWZk8eT8nazef4E9Z64ycuFuLJaC63w2WeaPHz/O5MmTGTFixA3r3377bb7++mt+/fVXevbsyXPPPce77757231NmDABHx8f6yMwMDBPMnq7O1PFzxPI/rRmwk9HeGLqFiKiE/Nk/yIiIiIiImKMxLRMBs4JZV34JQBcnRwY2LISJpOpwDKYLAX50cHfvPHGG7z55pu33WbHjh00bdrUuhwZGcm9997Lvffey8yZM2/7tZMmTeKtt94iLi7ultukpaWRlpZmXY6PjycwMJC4uDi8vb1z+J3cXGpGFpN+CWfmppNc/ym7ODnwQvuaDG5VGUeHgvsfLSIi9is+Ph4fH588OXbZE/1cREQkP1xNTmfAnB3sO3sVAE8XR2YOCKZF1ZJ5sv+cHr8MLfMxMTHExMTcdptKlSrh5uYGZBf5du3aERISwty5c3FwuP2FBZs3b6ZVq1ZERUVRunTpHGXKjwP/zlOxvLA8jJMxSdZ1TSqWYOIT9alSqlie/BsiIlJ0qbTenH4uIiKS1y4lpNFv1naORCUA4OPuzLzBzWgYWDzP/o2cHr8MHQDPz88PPz+/HG17/vx52rVrR5MmTZgzZ86/FnmAPXv24ObmRvHixe8y6d1pWsmXH59uzYe/hDN7c/ZZ+l2nr/Dgpxt5oUNNBt2js/QiIiIiIiKF2fmrKfSdud16ktavmCsLhjYjqIwxHxjbxGj2kZGRtG3blgoVKvDhhx9y6dIl63PXR65ftWoVUVFRtGjRAnd3d9atW8err77K8OHDcXV1NSq6lbuLI689UpuOdcvwwtf7OHU5mbRMM++sPsyaA1FM7NaAytfusRcREREREZHC41RMEn1mbuf81RQAyvq4sWBoiKFXWttEmf/ll1+IiIggIiKC8uXL3/Dc9bsEnJ2d+fLLL3n22Wcxm81UqVKFt956i1GjRhkR+ZaCK/ny09g2fPDzEeZuOYXFAjtPX+HBT//ghQ5BDGpZCQedpRcRERERESkUwqMS6DtrO5cSssdaq1TSg4XDmlOuuLuhuQy9Z74wKsj767afuMwLy8M4E5tsXdeski8fPFGfSjpLLyIiOaR7w29OPxcREblb+85eZcCcUK4mZwBQs7QX84c2w9/LLd/+zZwev2xyajp7EVKlJGueac3AlpWs60JPxdLx0z+Ys/kkZs1LLyIiIiIiYojQk7H0mbndWuQblPdh6ZPN87XI54bKvME8XJx4o1MdlgxvTgVfDwBSM8y8ueoQPWds4/TlpH/Zg4iIiIiIiOSljccu0X/2dhLTMgFoVtmXBUNDKO7hYnCyP6nMFxLNq5Tkp7Gt6d+ionVd6MlYHvx0Iwu2nUZ3Q4iIiIiIiOS/Xw9dZMjcnaRmmAFoU6MU8wY1w8vN2eBkN1KZL0Q8XZ14q3NdFg0LoXyJ7MEUktOzGP/dAQbO2UFUXKrBCUVEREREROzXqn2RjFiwi/Ss7CLfvnZpZvRvgruLo8HJ/kllvhBqWdWPn59pQ++QCtZ1G45eosMnf7ByX6SByUREREREROzT1zvPMnbJHjKvjV3WuWFZvujTGFenwlfkQWW+0PJ0deLdx+oxZ1Aw/l6uAMSlZPD04j2MWrSbK0npBicUERERERGxD/O3nuKF5WFcH4O8Z3AgH3VviLNj4a3MhTeZANCupj+//LcNnRqUta5bHXaB9p/8wboj0QYmExERERERsX3TNhznte8PWpcHtqzEu4/Vw9HBZGCqf6cybwOKe7jwWa9GTO7ViOIe2YMuXEpIY9DcHYxbEWYdYVFERERERERyxmKx8PHao0z46Yh13ci2VXn90do4FPIiDyrzNuXRBmX5+Zk2tK1ZyrpucehZHvz0D0JPxhqYTERERERExHZYLBYm/HSET387Zl33fPsavNgxCJOp8Bd5UJm3OaW93ZgzMJgJXevhcW1ExbOxKfSYvpV3fzxMakaWwQlFREREREQKL7PZwmvfH2D6Hyes6157pDaj/1PdwFS5pzJvg0wmE72aVeCnsa0JrlQCAIsFpv9xgk6fb+LA+TiDE4qISFEyYcIEgoOD8fLywt/fny5duhAeHn7DNiaT6aaPiRMn5ujfWLJkCSaTiS5duuTDdyAiIkVFZpaZF5aHsWDbGQBMJnj3sXoMaVXZ4GS5pzJvwyqW9GTJ8BaMezAIl2ujLB69mEiXLzbz+e/HyLw2N6KIiEh+2rBhA6NGjWLbtm2sXbuWzMxM2rdvT1JSknWbCxcu3PCYPXs2JpOJxx9//F/3f/r0aZ5//nlat26dn9+GiIjYufRMM2OX7OWb3ecAcHQw8VH3BjdMCW5LTBaLxWJ0iMIkPj4eHx8f4uLi8Pb2NjpOjh2JiufZpfs4dCHeuq5hYHE+7tGQyn6eBiYTEZH8VtiOXZcuXcLf358NGzbQpk2bm27TpUsXEhIS+O233267r6ysLO69914GDRrExo0buXr1Kt99912OchS2n4uIiBgnNSOLUQt389u1GcGcHU181rMRD9YLMDjZP+X0+KUz83YiqIw33426h9HtqnF94MW9Z6/y8GcbWRJ6Bn1mIyIiBSUuLvt2L19f35s+f/HiRVavXs2QIUP+dV9vvfUWpUqVytG2aWlpxMfH3/AQERFJTs9k6Lyd1iLv6uTA9H5NC2WRzw2VeTvi4uTA8x1qsvypllQq6QFAcnoWL6/Yz/D5u7icmGZwQhERsXcWi4Vnn32WVq1aUbdu3ZtuM2/ePLy8vOjatett97V582ZmzZrFjBkzcvRvT5gwAR8fH+sjMDAw1/lFRMS+xKdm0H9WKJsiYgDwcHFkzqBg2gX5G5zs7qnM26HGFUqw+unW9Gr255uYtYcu0uGTjawLjzYwmYiI2LvRo0cTFhbG4sWLb7nN7Nmz6dOnD25ubrfcJiEhgb59+zJjxgz8/Pxy9G+PGzeOuLg46+Ps2bO5zi8iIvbjanI6fWduZ+fpKwB4uTkxf0gzWlbN2XGlsHMyOoDkD09XJyZ0rU+7mv68vGI/sUnpxCSmMWjODga0qMi4h2rh5uxodEwREbEjY8aMYeXKlfzxxx+UL1/+ptts3LiR8PBwli5dett9HT9+nFOnTvHoo49a15nN2QO7Ojk5ER4eTtWqVW/4GldXV1xdXe/yuxAREXtwKSGNfrO2cyQqAYASHs7MHxJC3XI+BifLOyrzdq59nTI0rFCcF5eHsT78EgDztp5m8/HLfNKjoV39MouIiDEsFgtjxozh22+/Zf369VSufOvpfWbNmkWTJk1o0KDBbfcZFBTE/v37b1g3fvx4EhIS+PTTT3UJvYiI3NKFuBT6zNjOiZjsWVVKebmyYEgINct4GZwsb6nMFwH+Xm7MGRjMgm2neWf1YdIyzUREJ/LYl5t59oGaDG9TBcfro+aJiIjk0qhRo1i0aBHff/89Xl5eREVFAeDj44O7u7t1u/j4eL7++msmTZp00/3079+fcuXKMWHCBNzc3P5xz33x4sUBbnkvvoiIyJnLyfSeuY1zV1IAKOvjxsJhze1yhi/dM19EmEwm+rWoxOqnW1G3XPb0BhlZFt5fc4TeM7Zx/mqKwQlFRMRWTZkyhbi4ONq2bUtAQID18fdL6ZcsWYLFYqFXr1433c+ZM2e4cOFCQUQWERE7FBGdSPdpW61FvmJJD5aNaGGXRR40z/w/FIU5adMzzXz861GmbjjO9f/7Xm5OvNOlLp0bljM2nIiI5FpROHbdCf1cRESKjsMX4uk7czuXk9IBqOZfjIVDQyjtfevBVgsrzTMvt+Ti5MBLHYNYMqw55YpnX/6YkJrJ2CV7GbtkD3EpGQYnFBERERERyZl9Z6/Sc/o2a5GvHeDN0uHNbbLI54bKfBEWUqUkP45tTZeGZa3rvt8byYOf/MG2E5cNTCYiIiIiIvLvdp6Kpc/M7dYTkg0Di7N4WHNKFrP/2U1U5os4H3dnPunZiE97NsTLLXs8xMi4VHrN2MYHa46QkWU2OKGIiIiIiMg/bT1+mf6zQ0lMywQgpLIvC4aG4OPhbHCygqEyLwB0bliONc+0oXkVXwAsFvhy/XG6Td3KmcvJBqcTERERERH50x9HLzFwTijJ6VkAtK7ux9xBzSjmWnQmbFOZF6tyxd1ZOLQ5L3UMwunaVHV7z17loc828t2e8wanExERERERgd8OX2TovJ2kZWZfRfyfIH9m9G+Ku4ujwckKlsq83MDRwcRTbauy/KmWVCzpAUBiWibPLN3Ls8v2Wi9hERERERERKWhrDlxgxIJdpF+7HbhDndJM7dsEN+eiVeRBZV5uoWFgcVY/3Zqujf6cqm7F7vM88tlG9p29alwwEREREREpklbui2TUoj1kZGXPr/1og7J83rsxLk5Fs9YWze9acqSYqxMf9WjIJz0aWu89OXU5mcenbGHqhuOYzRaDE4qIiIiISFGwfNc5nlmyh6xrHeTxxuX5pEdDnB2LbqUtut+55FiXRuVY/XQrGgQWByDTbOG9n47Qf3Yo0fGpxoYTERERERG7tmj7GV5Yvo/r5xJ7NQtk4hP1cbw2zldRpTIvOVKxpCfLR7RgZNuqmK79zWyKiKHjpxv57fBFY8OJiIiIiIhdmrv5JK98ux/LtSI/sGUl3n2sHg5FvMiDyrzkgrOjAy92DGLhkBBKe7sCEJuUzpB5O3lj5UFSM7IMTigiIiIiIvZi+h/HeWPVIevy8DZVeP3R2phMKvKgMi93oGU1P34a24b7a5W2rpu75RRdvtjMsYsJBiYTERERERF7MPm3Y7z74xHr8pj/VGPcg0Eq8n+hMi93xNfThRn9m/B25zrW0SOPRCXw6OebWLj9NBaLBscTEREREZHcsVgsTPolnElrj1rXPfdADZ5rX1NF/m9U5uWOmUwm+rWoxMrR91CjdDEAUjPMvPrtAUYv3kN8aobBCUVERERExFZYLNkDbU/+PcK67pWHghhzX3UDUxVeKvNy14LKeLNydCv6Na9oXbc67AKPfLaJsHNXjQsmIiIiIiI2wWKx8OaqQ0z744R13Zud6jC8TVUDUxVuKvOSJ9ycHXm7S12m9m2Mt1v2nPRnYrPnpJ+16aQuuxcRERERkZsymy288u0B5m45BYDJBO8+Vo8BLSsZmquwU5mXPNWxbgCrn25Nw2tz0mdkWXj7h0MM+2oXV5PTjQ0nIiIiIiKFSpbZwgvLw1gcegYABxNMfKIBvUMqGJys8FOZlzwX6OvB1yNa8GSbKtZ1vx6+yEOfbmTnqVgDk4mIiIiISGGRkWXmmaV7+Wb3OQAcHUx80rMRTzQpb3Ay26AyL/nC2dGBcQ/VYs7AYEp4OAMQGZdKj+nb+HJ9BGazLrsXERERESmq0jPNjFm0h1X7IgFwdjTxRe/GdGpQ1uBktkNlXvJVuyB/fhzbmmaVfIHsy2g+WBPOgDmhxCSmGZxOREREREQKWmpGFiMW7GLNwSgAXJwcmNavCR3rljE4mW1RmZd8F+DjzqJhITz9n2pcnxpy47EYHvx0I1uOxxgbTkRERERECkxKehbDvtrJ70eiAXBzdmBm/6b8J6i0wclsj8q8FAgnRweebV+TBUNC8CvmCsClhDT6zNzOx2uPkqXL7kVERERE7FpSWiaD5oay8Vj2CT0PF0fmDGxGmxqlDE5mm1TmpUDdU82Pn8a2pnV1PwAsFvj0t2P0mbmNi/GpBqcTEREREZH8EJ+aQf/ZoWw7kT0gtperE/OHNKNF1ZIGJ7NdKvNS4Ep5uTJvUDNe6FATR4fs6+63nYjlwU83suHoJYPTiYiIiIhIXopLzqDfzO3sOn0FAG83JxYMDaFJRV+Dk9k2lXkxhIODiVHtqrFkeHMCfNwAiE1KZ+CcUD76JVyX3YuIiIiI2IErSen0nrmNfefiACjh4czi4c1pEFjc2GB2QGVeDBVcyZcfn27Nf4L8gezL7j/7PYL+s7drtHsRERERERt2OTGN3jO3czAyHgC/Yi4sGd6COmV9DE5mH1TmxXAlPF2Y2b8pL3UMsl52vzniMg99upHQk7EGpxMRERERkdyKSUyj94ztHL6QXeT9vVxZMrwFNct4GZzMfqjMS6Hg4GDiqbZVWTQ0BH+v7NHuoxPS6DVjG1M3HMesy+5FRERERGxCdEIqvaZvI/xiAgBlvN1YMrw51fyLGZzMvqjMS6ESUqUkq59uTctro1pmmS2899MRhs/fSVxyhsHpRERERETkdi7Gp9Jz+jaORScCUNbHjaVPNqdKKRX5vKYyL4VOKS9X5g8JYcx/qlnX/Xo4mocnbyTs3FXjgomIiIiIyC1diEuh5/RtnLiUBEC54u4sfbIFFUt6GpzMPqnMS6Hk6GDiufY1mTsomBIezgCcu5LCE1O2Mn/rKSwWXXYvIiIiIlJYnL+aQo9p2zgZk13kA33dWfpkcwJ9PQxOZr9U5qVQa1vTn9VPt6ZRheIApGeZee37gzy9ZC+JaZnGhhMREREREc7GJtNj2lbOxCYDULGkB0uHt6B8CRX5/KQyL4Ve2eLuLB3egsH3VLauW7Uvkk6fbyI8KsHAZCIiIiIiRduZy8n0nL6Nc1dSAKji58nS4S0oW9zd4GT2T2VebIKLkwP/92htpvRpjJerEwAnLiXR+YtNLN91zuB0IiIiIiJFz8mYJHpM38r5q9lFvmopT5YMb04ZHzeDkxUNKvNiUx6sF8CqMa2oFeANQGqGmee/3se4FWGkZmQZnE5EREREpGg4fimRntO3ciEuFYDq/sVYMrwF/t4q8gVFZV5sTiU/T74d2ZKewYHWdYtDz9J92p+fCoqIiIiISP6IiE6g5/RtXIxPAyCojBdLhjenlJerwcmKFpV5sUluzo6893h9PuzWAFen7F/jsHNxPPLZRjYdizE4nYiIiIiIfQqPyi7ylxKyi3ztAG8WDWtOyWIq8gVNZV5s2hNNyrNiZEsqXJvy4kpyBv1nb+eLdRGYzZq+TkREREQkrxy+EE+vGduISUwHoG45bxYNC8HX08XgZEWTyrzYvDplfVg1uhXtapYCwGyBiT+H8+SCXcSnZhicTkRERETE9h04H0evGduITcou8g3K+7BwSHOKe6jIG0VlXuyCj4czswYE89/7a2AyZa9be+ginT/frOnrRERERETuwv5zcfSZuZ2rydknyhpVKM78oSH4eDgbnKxoU5kXu+HgYGLs/dWZPTAYH/fsF5aTMUl0+WIz3+89b3A6ERERERHbs/fsVXrP3EZcSnaRb1qxBF8Nboa3m4q80VTmxe60q+nPqtGtqH1t+rqUjCzGLtnLm6sOkpFlNjidiIiIiIht2HX6Cv1mbichNROAZpV9mTe4GV4q8oWCyrzYpQolPVgxsiWPNy5vXTdn8yl6z9hGdHyqgclERERERAq/Hadi6T9rOwlp2UW+RZWSzB0UjKerk8HJ5DqVebFbbs6OfNitPv97rC7Ojtk30u84dYWHJ29ix6lYg9OJiIiIiBRO205cZsDsUJLSswBoXd2P2QOD8XBRkS9MVObFrplMJvqEVGTZky0I8HED4FJCGr2mb2P2ppNYLJq+TkRERETkui0RMQycE0rytSJ/b41SzOjfFHcXR4OTyd+pzEuR0KhCCVaNaUXLqiUByDRbeOuHQzyzdC8p116oRERERESKsj+OXmLQ3B2kZmSPM/WfIH+m9WuCm7OKfGGkMi9Fhl8xV74a3IwR91a1rvt+byRPTN3CuSvJBiYTERERETHW+vBohn61k7TM7CL/QO3STOnbWEW+EFOZlyLFydGBlx8MYmrfxnhcu1ToYGQ8nT7fzJbjMQanExEREREpeL8fucjwr3aRfq3Id6xThi96N8bVSUW+MFOZlyKpY90Avh15DxVLegAQm5ROv1mhuo9eRERERIqUXw9d5Mn5u0i/NoXzw/UCmNy7ES5OqoqFnf4PSZFVs4wXK0e14t4apQDIunYf/XNf7yM1Q/fRi4iIiIh9W3voIk8t3EVGVvbJrEcblOXTng1xdlRNtAX6vyRFmo+HM7MHBvNU2z/vo1+x+zzdp20l8mqKgclERERERPLPLwejGPmXIt+5YVk+7t4AJxV5m6H/U1LkOTqYeKljEF/0boz7tQE+ws7F0enzTWw/cdngdCIiIiIieeuXg1GMWrTbWuS7NCzLR90bqsjbGP3fErnm4foBrBjZkkBfdwBiEtPpM3M7X209pfvoRURERMQu/HwwipEL/yzyjzUqx6TuDXF0MBmcTHJLZV7kL2oFeLNqdCtaV/cDsuej/7/vD/LSN2G6j15EREREbNqaA1GMWribTHN2ke/aqBwfdmugIm+jVOZF/qa4hwtzBgYzvE0V67plO8/RY/o2ouJSDUwmIiIiInJn1hy4wOhFfynyjcsxUUXepqnMi9yEk6MDrzxUi097NsTNOfvPZN/ZqzwyeRM7T8UanE5EREREJOd+2n+B0Yv2WIv8443LM/EJFXlbpzIvchudG5bjm6daUq749fvo0+g1YxtLd5wxOJmIiIiIyL/7cf8FRi/+s8g/0aQ8HzxRX0XeDqjMi/yLOmV9WDWmFS2rlgQgI8vCS9/s581VB8nMMhucTkRERETk5laHXWDM4j1kXSvy3ZqU5/3HVeTthcq8SA74errw1eBmDGxZybpuzuZTDJq7g7jkDOOCiYiIiIjcxA9hkTy9REXenqnMi+SQk6MDb3Sqw3td6+HsmP0iuPFYDI99uZnjlxINTiciIiIikm3VvkjGLtlrLfLdm2YXeQcVebuiMi+SSz2bVWDh0Ob4eroAcCImiS5fbGbD0UsGJxMRERGRom7VvkieWfpnke/RNJD3uqrI2yOVeZE70KyyL9+PuoegMl4AJKRmMmhOKLM2ncRisRicTkRERESKopX7Ihn7l0vrewYHMqFrPRV5O6UyL3KHAn09+OaplrSvXRoAswXe/uEQL30TRlpmlsHpREQKzoQJEwgODsbLywt/f3+6dOlCeHj4DduYTKabPiZOnHjL/c6YMYPWrVtTokQJSpQowf33309oaGh+fzsiIjbp+73neWbJHq71eHo1C+Tdx1Tk7ZnKvMhd8HR1YmrfJoz5TzXrumU7z9FnxnZiEtMMTCYiUnA2bNjAqFGj2LZtG2vXriUzM5P27duTlJRk3ebChQs3PGbPno3JZOLxxx+/5X7Xr19Pr169WLduHVu3bqVChQq0b9+e8+fPF8S3JSJiM77fe57/Lt37lyJfgf91UZG3dyaLrgm+QXx8PD4+PsTFxeHt7W10HLEhq/ZF8vzX+0jLzJ6urlxxd6b3b0Kdsj4GJxMRe1fYjl2XLl3C39+fDRs20KZNm5tu06VLFxISEvjtt99yvN+srCxKlCjB559/Tv/+/f91+8L2cxERyQ/f7TnPs8v+LPK9QyrwTue6KvI2LKfHL52ZF8kjjzYoy/IRLSnj7QbA+aspPDFlK2sOXDA4mYhIwYqLiwPA19f3ps9fvHiR1atXM2TIkFztNzk5mYyMjFvuNy0tjfj4+BseIiL27Ns9524o8n1U5IsUlXmRPFSvvA8rR99Dg8DiAKRkZDFiwW4+++2YBsYTkSLBYrHw7LPP0qpVK+rWrXvTbebNm4eXlxddu3bN1b5ffvllypUrx/3333/T5ydMmICPj4/1ERgYmOv8IiK24ts953hu2T5rke/bvAJvq8gXKTZR5k+dOsWQIUOoXLky7u7uVK1alddff5309PR/bDt37lzq16+Pm5sbZcqUYfTo0QYklqLM39uNpcOb81ijctZ1H609ypjFe0jN0MB4ImLfRo8eTVhYGIsXL77lNrNnz6ZPnz64ubnleL8ffPABixcvZsWKFbf8unHjxhEXF2d9nD17Ntf5RURswYrd53j2L0W+X/OKKvJFkJPRAXLiyJEjmM1mpk2bRrVq1Thw4ADDhg0jKSmJDz/80LrdRx99xKRJk5g4cSIhISGkpqZy4sQJA5NLUeXm7MhH3RtQs4wX7685gsUCP4Rd4PzVFKb3a0opL1ejI4qI5LkxY8awcuVK/vjjD8qXL3/TbTZu3Eh4eDhLly7N8X4//PBD3n33XX799Vfq169/y+1cXV1xddXrq4jYt+W7zvHC8n1cv+izf4uKvNmpDiaTinxRY7MD4E2cOJEpU6ZYy/qVK1coV64cq1at4r777rvj/WqwHMlrvx66yNNL9pCcnn1Wvlxxd2YPDKbmtTnqRUTultHHLovFwpgxY/j2229Zv3491atXv+W2AwcO5MCBA+zcuTNH+544cSLvvPMOP//8M82bN89VLqN/LiIiee3rnWd58Zswa5Ef0KIib6jI2x27HwAvLi7uhgFw1q5di9ls5vz589SqVYvy5cvTvXv3f73EToPlSH67v3Zpvh7RggCfPwfGe3zKFtaFRxucTEQkb4waNYoFCxawaNEivLy8iIqKIioqipSUlBu2i4+P5+uvv2bo0KE33U///v0ZN26cdfmDDz5g/PjxzJ49m0qVKln3m5iYmK/fj4hIYbTsb0V+YMtKKvJFnE2W+ePHjzN58mRGjBhhXXfixAnMZjPvvvsun3zyCcuXLyc2NpYHHnjgpvfWX6fBcqQg1Cnrw/ej7qF++exp6hLTMhkydwdzN580OJmIyN2bMmUKcXFxtG3bloCAAOvj75fSL1myBIvFQq9evW66nzNnznDhwp8zgHz55Zekp6fzxBNP3LDfv95iJyJSFCzbcZaX/lbkX3+0top8EWfoZfZvvPEGb7755m232bFjB02bNrUuR0ZGcu+993Lvvfcyc+ZM6/p3332XV199lZ9//pn27dsD2fPclilThh9//JEOHTrcdP9paWmkpaVZl+Pj4wkMDNQleZIvUtKz+O/Svaw5GGVd179FRf7vkdo4OdrkZ2siUgjocvKb089FROzBsh1neWnFn0V+0D2V+L9HVOTtWU6PX4YOgDd69Gh69ux5220qVapk/e/IyEjatWtHixYtmD59+g3bBQQEAFC7dm3rulKlSuHn58eZM2duuX8NliMFyd3FkS/7NObDX8L5cv1xAL7aeprTl5OZ3LsR3m7OBicUERERkcJi2c4bi/zgeyrz2iO1VOQFMLjM+/n54efnl6Ntz58/T7t27WjSpAlz5szBweHGs5j33HMPAOHh4dYRdGNjY4mJiaFixYp5G1zkLjg4mHixYxCV/Tx55dv9ZGRZ2HD0Ek9M2cKsAcEE+noYHVFEREREDLZ817kbLq1XkZe/s4nreiMjI2nbti2BgYF8+OGHXLp0yToIznU1atSgc+fOjB07li1btnDgwAEGDBhAUFAQ7dq1MzC9yM11axrI/CEhFPfIPht/9GIiXb7YzK7TVwxOJiIiIiJG+nbPjdPPDbqnkoq8/INNlPlffvmFiIgIfv/9d8qXL3/DIDh/9dVXXxESEsLDDz/Mvffei7OzM2vWrMHZWZcuS+HUvEpJvh15D1X8PAG4nJROrxnb+H7veYOTiYiIiIgRvttznueW7bthsDvdIy83Y7PzzOcXDZYjRrianM5TC3az9cRl67pn7q/O2Puq64VbRP6Vjl03p5+LiNia7/ee579L92K+1tD6Na/IW501/VxRY/fzzIvYk+IeLswb3IweTf+cGvGTX4/x36V7ScvMMjCZiIiIiBSEVfsibyjyfZtXUJGX21KZFykkXJwceO/xerzyUBDXX7O/2xtJ/1mhxCVnGBtORERERPLND2GRPPOXIt87pAJvdaqrIi+3pTIvUoiYTCaGt6nKlD5NcHPO/vPcfjKWrlM2czY22eB0IiIiIpLXftx/gbFL9pJ1rcn3ahbIO53r4uCgIi+3pzIvUgh1rFuGxcOaU9LTBYDjl5J47MvN7Dt71dhgIiIiIpJn1hy4wNOL91iLfI+mgfyvSz0VeckRlXmRQqpRhRLZI92Xyh7pPiYxnZ7Tt7H20EWDk4mIiIjI3fr5YBSjF+0h81qR79akPBO6qshLzqnMixRiFUp6sOKpljSr5AtASkYWT87fybwtp4wNJiIiIiJ37JeDUYxauNta5B9vXJ73H6+vIi+5ojIvUsgV93DhqyHNeLRBWQDMFnh95UHe+eEQZrNmlhQRERGxJb8eusioRX8W+a6NyvHBEyryknsq8yI2wM3ZkU97NGRk26rWdTM3nWTUot2kZmjqOhERERFb8PuRizy1cBcZWdlFvkvDskzs1gBHFXm5AyrzIjbCwcHEix2DePexetYX/J8ORNF7xjYuJ6YZnE5EREREbmddeDQj5u+2FvnODcsyqXtDFXm5YyrzIjamd0gFZg1oiqeLIwC7z1yl65QtnIxJMjiZiIiIiNzM+vBonpy/i/QsMwCPNijLJJ2Rl7ukMi9ig9rW9GfZiBaU9nYF4PTlZLp+uZmdp2INTiYiIiIif/XH0UsMn7+L9MzsIv9w/QA+7t4AJ0dVMbk7+g0SsVF1yvrw7ch7CCrjBcCV5Ax6z9zOj/svGJxMRERERAA2HrvEsK92/lnk6wXwaY+GKvKSJ/RbJGLDyhZ3Z9mIFrSu7gdAeqaZUYt2M3fzSYOTiYiIiBRtmyNiGDpvJ2nXivyDdcvwSU8Veck7+k0SsXHebs7MHhhMtyblAbBY4I1Vh3jvpyNYLJq6TkRERKSgbYmIYci8HdYi36FOaT7r1QhnFXnJQ/ptErEDzo4OfPBEfUa3q2ZdN3XDcZ5bto+MawOtiIiIiEj+23r8MoPn7SA1I/s92AO1SzO5V2MVeclz+o0SsRMmk4nnO9Tk7S51MV0bGHXFnvMMnruDxLRMY8OJiIiIFAHbT1xm8Nw/i/z9tfz5ondjXJxUuyTv6bdKxM70a16RKX2aWA8aG4/F0HP6Vi4laC56ERERkfwSejKWQXN3kJKRBcB9Qf580UdFXvKPfrNE7FDHumVYODQEH3dnAA6cj+dxzUUvIiIiki92nIpl4JxQktOzi3y7mqX4sm9jXJ0cDU4m9kxlXsROBVfyZfmIFpT1cQPgTGwyj0/Zwt6zV40NJiIiImJHdp2OZeDsP4t825qlmNK3iYq85DuVeRE7Vr20Fyv+Mhd9bFI6vaZvY92RaIOTiYiIiNi+3WeuMGD2DpKuFfk2NUoxtW8T3JxV5CX/qcyL2LkyPm4sfbIFzav4ApCSkcXQr3aybMdZg5OJiIiI2K49Z64wYFaodaDh1tX9mN5PRV4Kjsq8SBHg4+7MvMHNeLh+AABZZgsvfhPG5N+OaS56ERERkVwKO3eV/rNDSbhW5FtV82NG/6Yq8lKgVOZFighXJ0cm92zEoHsqWddNWnuU8d8dIMusQi8iIiKSEwfOx9FvVigJqdlFvmXVkiryYgiVeZEixMHBxP89UptxDwZZ1y3cfoYxi3eTlpllYDIRERGRwu9IVDz9Zm0nLiUDgGaVfZk5oCnuLiryUvBU5kWKGJPJxJP3VuXjHg1wcjAB8OP+KIbM3Wm950tEREREbnTsYgJ9ZmznSnJ2kW9asQRzBgbj4eJkcDIpqlTmRYqoxxqVZ8aAprg5Z78MbIqIoc+MbcQmpRucTERERKRwiYhOpNeM7Vy+9j6pUYXizBkUjKerirwYR2VepAhrV9OfhUND8HbLPhDtOxfHE1O3cP5qisHJRERERAqHkzFJ9J6xjZjENADql/dh3uBmeLk5G5xMijqVeZEirklFX74e0RJ/L1cATlxK4okpW4iITjA4mYiIiIixzlxOpveMbUQnZBf52gHezB8cgreKvBQCKvMiQs0yXnzzVEsqlfQA4EJcKt2mbmXv2avGBhMRERExyLkryfSasY0LcakABJXxYuHQEHw8VOSlcFCZFxEAAn09WP5US+qU9QbgSnIGvWdsY+OxSwYnExERESlYkVdT6DVjm/XWw+r+xVgwNIQSni4GJxP5k8q8iFj5FXNlyfDmhFT2BSA5PYvBc3fwQ1ikwclERERECkZUXCq9Z2zjbGx2ka9SypOFw0LwK+ZqcDKRG6nMi8gNvNycmTe4Ge1rlwYgI8vCmMV7mL/ttMHJRERERPJXdEIqvWdu49TlZAAq+3myeFhz/L3cDE4m8k8q8yLyD27OjnzZpzE9mgYCYLHAa98d4NNfj2GxWAxOJyIiIpL3YhLT6DNjOycuJQFQwdeDRcNCKO2tIi+Fk8q8iNyUk6MD7z1ejxH3VrWu+/jXo7y56hBmswq9iIiI2I/YpHT6ztzOsehEAMoVd2fRsBACfNwNTiZyayrzInJLJpOJlx8M4pWHgqzr5m45xfNf7yMzy2xgMhEREZG8cTU5u8gficqeljfAx43Fw5pTvoSHwclEbk9lXkT+1fA2VZn4RH0cHUwArNhznlGLdpOWmWVwMhEREZE7F5eSQb9ZoRy6EA9AaW9XFg9rToWSKvJS+KnMi0iOdGsayBe9G+HsmF3ofz54kaHzdpKcnmlwMhEREZHcS0jNYMDsUPafjwOglJcri4Y1p5Kfp8HJRHJGZV5Ecqxj3QBmDgjGzTn7pWPjsRj6zwolPjXD4GQiIiIiOZeUlsnAOTvYe/YqACU9XVg0NISqpYoZG0wkF1TmRSRX7q1RivlDQvBydQJg5+kr9J6xjcuJaQYnExEREfl3yemZDJq7g12nrwBQwsOZhcNCqF7ay+BkIrmjMi8iuRZcyZfFw5vj6+kCwIHz8fSYvo2ouFSDk4mIiIjcWkp6FkPn7ST0ZCwAPu7OLBgaQlAZb4OTieSeyryI3JG65XxY9mRzSnu7AhARnUi3aVs4cznZ4GQiIiIi/5SakcXw+TvZcvwyAF5uTiwYEkKdsj4GJxO5MyrzInLHqvl7sXxESwJ9s+dgPRubQrdpWzh2McHgZCIiIiJ/SsvMYsSCXWw8FgNAMVcnvhrcjHrlVeTFdqnMi8hdCfT14OsnW1LNP3vAmIvxaXSftpX95+IMTiYiIiIC6ZlmRi3czfrwSwB4uDgyb3AwjSqUMDiZyN1RmReRu1bGx41lT7agbrns+82uJGfQe8Y26/1oIiIiIkbIyDLz9OI9/Ho4GgB3Z0fmDAymSUVfg5OJ3D2VeRHJE76eLiwa1pzgStmfciekZdJ/9nY2HL1kcDIREREpijKzzDyzdC9rDkYB4OrkwKwBTQmpUtLgZCJ5Q2VeRPKMt5szXw0OoXV1PwBSM8wMnbeDX64dREWkcAgPD+eNN97gvvvuo2rVqgQEBFC/fn0GDBjAokWLSEvTVJMiYtuyzBae+3ofq8MuAODi5MDMAU1pWc3P4GQieUdlXkTylLuLIzMHNKVjnTIAZGRZGLlwt/VgKiLG2bNnDw888AANGjTgjz/+IDg4mGeeeYa3336bvn37YrFYePXVVylbtizvv/++Sr2I2CSz2cKLy8P4fm8kAC6ODkzr14TW1UsZnEwkbzkZHUBE7I+rkyOf927EC8vD+HbPeTLNFsYs3k16VgMea1Te6HgiRVaXLl144YUXWLp0Kb6+t75fdOvWrXz88cdMmjSJV155pQATiojcHbPZwivf7ueb3ecAcHY08WWfxrSr6W9wMpG8pzIvIvnCydGBD7s1wNnRxLKd5zBb4Nll+8jItNA9ONDoeCJF0rFjx3BxcfnX7Vq0aEGLFi1IT08vgFQiInnDYrHw+sqDLNlxFgBHBxOTezXm/tqlDU4mkj90mb2I5BtHBxPvda1P3+YVALBY4MVvwpi/7bTByUSKppwU+bvZXkTEKBaLhXdWH7a+x3Awwac9G9KxbhmDk4nkH52ZF5F85eBg4u3OdXFxdGT25pMAvPbdAdIzzQxpVdngdCJF0+XLlwkLC6NBgwb4+voSExPDrFmzSEtLo1u3btSqVcvoiCIiOWaxWJj4czizNmW/zzCZYFL3BjxSv6zByUTyl8q8iOQ7k8nEa4/UwtXZgSnrjwPw9g+HSM8081TbqganEylaQkNDad++PfHx8RQvXpy1a9fSrVs3nJycsFgsvPfee2zatInGjRsbHVVEJEcm/x7Bl9feXwC817WexuiRIkGX2YtIgTCZTLzYoSbP3F/duu79NUf45NejWCwWA5OJFC2vvvoq3bp1Iy4ujldeeYUuXbpw3333cfToUY4dO0bv3r15++23jY4pIpIj0zYc56O1R63Lb3euQ4/gCgYmEik4KvMiUmBMJhPP3F+DFzvWtK775NdjTPw5XIVepIDs2rWLZ599Fi8vL8aOHUtkZCTDhg2zPj9q1Ch27NhhYEIRkZyZu/kkE346Yl0e/3At+rWoZFwgkQKmMi8iBW5k22q89kht6/KX64/zzurDKvQiBSA9PR13d3cAnJ2d8fDwwM/Pz/p8yZIluXz5slHxRERyZNH2M7yx6pB1+YUONRnauoqBiUQKnsq8iBhiSKvKvN2lrnV51qaTvPb9AcxmFXqR/BQYGMiJEyesy0uWLCEgIMC6fOHChRvKvYhIYbN81zle/W6/dfnp+6ozql01AxOJGEMD4ImIYfo1r4irowMvrQjDYoEF286QkWlhQtd6ODiYjI4nYpd69uxJdHS0dfnhhx++4fmVK1fSrFmzgo4lIpIjK/dF8uLyfVy/mO/JNlX471/G4xEpSnRmXkQM1T04kI+7N+R6d1+68ywvrwjTGXqRfPL666/Ts2fPWz7/6quvsmjRolztc8KECQQHB+Pl5YW/vz9dunQhPDz8hm1MJtNNHxMnTrztvr/55htq166Nq6srtWvX5ttvv81VNhGxH2sORPHfpXu5/hZhYMtKvPxgECaTTgBI0aQyLyKG69KoHJN7NcbxWqNftvMcL32jQi9iBA8PD1xdXXP1NRs2bGDUqFFs27aNtWvXkpmZSfv27UlKSrJuc+HChRses2fPxmQy8fjjj99yv1u3bqVHjx7069ePffv20a9fP7p378727dvv+PsTEdu07kg0YxbvJuvae4NezSrw+qO1VeSlSDNZNOLUDeLj4/Hx8SEuLg5vb2+j44gUKT/uv8CYxXusB+puTcrz/uP1dcm9yL+4k2OXxWJh+fLlrFu3jujoaMxm8w3Pr1ix4o7zXLp0CX9/fzZs2ECbNm1uuk2XLl1ISEjgt99+u+V+evToQXx8PD/99JN1XceOHSlRogSLFy/+1xw6povYh03HYhg8bwfpmdmvU483Ls/EJ/T+QOxXTo9fOjMvIoXGQ/UCmNyrkfUM/de7zvHiN2HWci8ieWfs2LH069ePkydPUqxYMXx8fG543I24uDgAfH19b/r8xYsXWb16NUOGDLntfrZu3Ur79u1vWNehQwe2bNly0+3T0tKIj4+/4SEitm37icsM/erPIv9og7J8oCIvAmgAPBEpZB6qF4AJGLN4D5lmC8t3nQPg/cfrW0u+iNy9BQsWsGLFCh566KE83a/FYuHZZ5+lVatW1K1b96bbzJs3Dy8vL7p27XrbfUVFRVG6dOkb1pUuXZqoqKibbj9hwgTefPPNOwsuIoXOrtNXGDx3B6kZ2UW+Q53SfNS9gd4PiFyjM/MiUug8eO0MvdO1g/XyXed4cbnO0IvkJR8fH6pUyfs5mUePHk1YWNhtL4OfPXs2ffr0wc3N7V/39/f7YS0Wyy3vkR03bhxxcXHWx9mzZ3MXXkQKjbBzVxk4O5Sk9CwA/hPkz+RejXF2VH0RuU5/DSJSKD1YL4DPe/9Z6L/ZfY4Xlu9ToRfJI2+88QZvvvkmKSkpebbPMWPGsHLlStatW0f58uVvus3GjRsJDw9n6NCh/7q/MmXK/OMsfHR09D/O1l/n6uqKt7f3DQ8RsT2HIuPpNyuUhLRMAFpV8+PLPo1xcVJ1Efkr/UWISKHVsW4An/dubC30K3afV6EXySPdunXjypUr+Pv7U69ePRo3bnzDIzcsFgujR49mxYoV/P7771SuXPmW286aNYsmTZrQoEGDf91vixYtWLt27Q3rfvnlF1q2bJmrfCJiO45dTKDfrO3EpWQA0KyyL9P7N8HN2dHgZCKFj+6ZF5FCrWPdMnzeuzGjF+0m02xhxe7zYIGJ3XTPnMjdGDhwILt27aJv376ULl36rqZ3GjVqFIsWLeL777/Hy8vLejbdx8cHd3d363bx8fF8/fXXTJo06ab76d+/P+XKlWPChAlA9iB9bdq04f3336dz5858//33/Prrr2zatOmOs4pI4XUyJoneM7dzOSkdgEYVijN7YDAeLqosIjejvwwRKfT+Uej3nAdU6EXuxurVq/n5559p1arVXe9rypQpALRt2/aG9XPmzGHgwIHW5SVLlmCxWOjVq9dN93PmzBkcHP68aLBly5YsWbKE8ePH89prr1G1alWWLl1KSEjIXWcWkcLlbGwyvWds41JCGgD1yvkwd1AzirmqrojciuaZ/xvNSStSeP18MIpRC7MLPcBjjcrxoQq9yB0du4KCgli2bBn169fP53TG0TFdxDZEXk2hx/StnI3NHsMjqIwXi4c1p4Sni8HJRIyheeZFxO50qFOGL/r8eQ/9t3uy76E36x56kVybNGkSL774IqdOnTI6iogUYdHxqfSZud1a5Kv5F2PB0BAVeZEc0HUrImJTOtQpw5d9GjNq0W4ysrLvoXdxdODdx+rhoDP0IjnWt29fkpOTqVq1Kh4eHjg7O9/wfGxsrEHJRKSoiElMo/fM7ZyMSQKgUkkPFg0Nwa+Yq8HJRGzDXZf5tLQ0XF31ByciBad9nTJM7pVd6LPMFpbsOIuTo4m3O9e9q0G8RIqSjz/+WH8vImKYq8np9J25nYjoRADKl3Bn0bDm+Hu7GZxMxHbkusz//PPPLF68mI0bN3LmzBnMZjMeHh40btyY9u3bM2jQIMqWLZsfWUVErDrWLcMnPRoydskezBZYsO0Mzo4O/N8jtVVQRG7jl19+oV27djcMTCciUpDiUzPoNyuUI1EJAAT4uLF4WHPKFnf/l68Ukb/K8T3z3333HTVr1mTAgAE4ODjwwgsvsGLFCn7++WdmzZrFvffey6+//kqVKlUYMWIEly5dys/cIiI82qAsk7o34Hp3n7P5FO/9dASN6ylyayNGjKBUqVL06NGDRYsWcfXqVaMjiUgRkpiWycDZoew/HwdAKS9XFg4NIdDXw+BkIrYnx6PZN2vWjNdee42HH374hmlj/u78+fN8+umnlC5dmueeey7PghYUjXwrYnuW7TzLi8vDrMtj/lON59rXNDCRSMHK7bErLCyMlStXsnLlSsLCwrjnnnvo3LkznTp1olKlSvkfuIDomC5SuKSkZzFgTiihJ7PH5PD1dGHp8OZUL+1lcDKRwiWnxy9NTfc3OvCL2KaF20/z6rcHrMvPPlCDp++rbmAikYJzN8euyMhIa7Fft24dNWrUsBb7pk2b5lPigqFjukjhkZqRxbCvdrLxWAwAPu7OLB7WnNpl9bcp8neamk5EipQ+IRV549Ha1uWP1h5lyvrjBiYSsQ1ly5ZlxIgR/Pjjj8TExPDaa69x6tQpOnbsyLvvvmt0PBGxA+mZZkYu3G0t8l6uTswf0kxFXuQu3dFo9haLheXLl7Nu3Tqio6Mxm803PL9ixYo8CScikhsD76lMptnCO6sPA/D+miM4O5oY2rqKwclEbIOnpydPPPEETzzxBGazmcuXLxsdSURsXGaWmbFL9vD7kWgAPFwcmTu4GfXLFzc2mIgduKMyP3bsWKZPn067du0oXbq0Ro4WkUJjaOsqpGWamfhzOADvrD6Mi5MD/VtUMjaYSCHz2Wef3XS9yWTCzc2N6tWr07p16wJOJSL2xGy28MLyMH46EAWAm7MDswcG06RiCYOTidiHOyrzCxYsYMWKFTz00EN5nUdE5K6NaleNjCwzn/x6DID/+/4gzo4O9GpWweBkIoXHxx9/zKVLl0hOTqZEiRJYLBauXr2Kh4cHxYoVIzo6mipVqrBu3ToCAwONjisiNsZisfDqdwf4ds95AFwcHZjerynNq5Q0OJmI/bije+Z9fHyoUkWXrYpI4TX2vuqMalfVuvzKt/tZvuucgYlECpd3332X4OBgjh07xuXLl4mNjeXo0aOEhITw6aefcubMGcqUKcN///tfo6OKiI2xWCy8/cNhFoeeAcDRwcTnvRvRpkYpg5OJ2Jc7KvNvvPEGb775JikpKXmdR0QkT5hMJp5vX5PhbbI/eLRY4MXl+/hx/wWDk4kUDuPHj+fjjz+matU/P/SqVq0aH374IePGjaN8+fJ88MEHbN682cCUImKLJv1ylNmbTwJgMsHHPRrSvk4Zg1OJ2J87usy+W7duLF68GH9/fypVqoSzs/MNz+/evTtPwomI3A2TycS4B4NIy8hi3tbTmC0wdske3F0caVfT3+h4Ioa6cOECmZmZ/1ifmZlJVFT2/a1ly5YlISGhoKOJiA37Yl0En6+LsC6//3h9OjUoa2AiEft1R2V+4MCB7Nq1i759+2oAPBEp1EwmE68/Woek9CyW7zpHRpaFEfN3MW9wM923J0Vau3btePLJJ5k5cyaNGjUCYM+ePTz11FP85z//AWD//v1UrlzZyJgiYkNmbzppHYAW4K3OdejeVGNuiOSXOyrzq1ev5ueff6ZVq1Z5nUdEJM85OJh4//H6pKRnsXr/BdIyzQyZu4OFw5rTMLC40fFEDDFr1iz69etHkyZNrFfYZWZmct999zFr1iwAihUrxqRJk4yMKSI2YknoGd764ZB1+eUHgzSTjEg+u6MyHxgYiLe3d15nERHJN44OJj7u0ZDk9EzWhV8iKT2LAbNDWTK8ObUC9HomRU+ZMmVYu3Yt4eHhhIeHY7FYCAoKombNmtZt2rVrZ2BCEbEV3+05z7hv91uXx95XnRH3Vr3NV4hIXrijAfAmTZrEiy++yKlTp/I4johI/nFxcmBK3yY0r+ILQFxKBv1mbefEpUSDk4kUvLCwMABq1qxJp06d6Ny5s7XIf/fddwYmExFbsubABZ77eh8WS/by8DZVeOb+6saGEiki7qjM9+3bl3Xr1lG1alW8vLzw9fW94SEiUli5OTsyc0Cw9fL6mMR0+s7czrkrycYGEylgHTp04MSJE/9Y/80339CnTx8DEomIrVl3JJoxi/eQZc5u8v2aV2Tcg0EaT0ukgNzRZfaffPJJHscQESk4xVydmDeoGT2mb+VIVAKRcan0mbmdr59sgb+3m9HxRArEU089xX333ceWLVsICAgAYOnSpQwePJi5c+caG05ECr0tx2MYsWAXGVnZRf6JJuV5s1MdFXmRAmSyWK5fFCMA8fHx+Pj4EBcXp3EBROzcpYQ0ekzbyomYJABqlC7GkuEt8PV0MTiZSO7c6bFr7Nix/Prrr2zcuJE1a9YwdOhQ5s+fz+OPP56PaQuOjuki+WPX6Vj6zQolOT0LgIfrB/BZz0Y4OqjIi+SFnB6/cnyZfVJSUq4C5Hb72zl16hRDhgyhcuXKuLu7U7VqVV5//XXS09Ot28ydOxeTyXTTR3R0dJ5lERH7UcrLlQVDQyhX3B2AoxcTGTA7lITUDIOTiRSMTz/9lMaNG9O8eXOGDRvG4sWL7abIi0j+OHA+joGzd1iL/P21/PmkR0MVeRED5Pgy+2rVqjFmzBgGDhxI2bJlb7qNxWLh119/5aOPPqJNmzaMGzcuT0IeOXIEs9nMtGnTqFatGgcOHGDYsGEkJSXx4YcfAtCjRw86dux4w9cNHDiQ1NRU/P398ySHiNifssXdWTg0hG7TtnIpIY395+MYPHcH8wY3w8Plju5EEim0Vq5c+Y91Xbp0YcOGDfTq1QuTyWTdplOnTgUdT0QKufCoBPrN2k5CWiYArav78Xnvxjg73tEwXCJyl3J8mX14eDjjx49n5cqVNGzYkKZNm1K2bFnc3Ny4cuUKhw4dYuvWrTg7OzNu3DiGDx+Oo6NjvgWfOHEiU6ZMuengPQCXLl2iXLly1nl0byUtLY20tDTrcnx8PIGBgbokT6SIOXoxgR7TtnIlOfus/L01SjGjf1NcnPQGRQq/HF+O55Cz32eTyURWVlZexTOMLrMXyTsnY5LoNnUrMYnZ75uDK5XQB98i+STPL7OvWbMmX3/9NcePH6dnz55ERkayfPlyZsyYwfr16ylXrhwzZszg1KlTPPXUU/la5AHi4uJuO3L+V199hYeHB0888cRt9zNhwgR8fHysj8DAwLyOKiI2oEZpL74aHIKXa/abkg1HL/H81/swmzWsiNgPs9mco4c9FHkRyTvnriTTZ8Y2a5FvUN6H2QODVeRFDGaTA+AdP36cxo0bM2nSJIYOHXrTberUqcO9997Ll19+edt96cy8iPzV9hOX6T87lLRMMwADWlTkDY3OK4WczkDfnH4uInfvYnwq3aZu5Uxs9hSuQWW8WDK8OcU9NFisSH7J8zPz+eGNN9645aB11x87d+684WsiIyPp2LEj3bp1u2WR37p1K4cOHWLIkCH/msHV1RVvb+8bHiJSdIVUKcnnvRtbB/KZt/U0n/x6zOBUInljyZIlOd727NmzbN68OR/TiEhhdzkxjT4zt1uLfJVSniwYGqIiL1JIGFrmR48ezeHDh2/7qFu3rnX7yMhI2rVrR4sWLZg+ffot9ztz5kwaNmxIkyZNCuLbEBE780Dt0nzweH3r8qe/HWPu5pMGJhLJG1OmTCEoKIj333+fw4cP/+P5uLg4fvzxR3r37k2TJk2IjY01IKWIFAZxyRn0nRVKRHQiAIG+7iwa2hy/Yq4GJxOR6wy90cXPzw8/P78cbXv+/HnatWtHkyZNmDNnzi0H8UlMTGTZsmVMmDAhL6OKSBHzeJPyXE3J4O0fDgHwxqpDFPdwoUujcgYnE7lzGzZs4IcffmDy5Mm88soreHp6Urp0aetgtlFRUZQqVYpBgwZx4MABzQYjUkQlpmUyYE4ohy/EAxDg48aioc0p4+NmcDIR+SubGLUiMjKStm3bUqFCBT788EMuXbpkfa5MmTI3bLt06VIyMzPp06dPQccUETszpFVlriSl8/m6CACe/3ofPu7OtAtSwRHb9cgjj/DII49w+fJlNm3axKlTp0hJScHPz49GjRrRqFGjHI96LyL2JyU9i8Fzd7D37FUA/Iq5sGBoCIG+HsYGE5F/sIky/8svvxAREUFERATly5e/4bm/j983a9YsunbtSokSJQoyoojYqefa1yA2OZ1F28+Qabbw1MJdLBgSQtNKt55NQ8QWlCxZks6dOxsdQ0QKkbTMLIbP30noyexbbIp7OLNgaAhVSxUzOJmI3MwdjWZ/7Ngxvv/+e06dOoXJZKJy5cp06dKFKlWq5EfGAqWRb0Xk77LMFp5esofVYRcA8HZzYumTLagVoNcIKRx07Lo5/VxEci4jy8zIhbtZe+giAF6uTiwcFkL98sWNDSZSBOXbaPYTJkygdu3avPTSS3zzzTd8/fXXvPDCCwQFBfHhhx/eVWgRkcLI0cHEx90b0rp69hgf8amZ9J8dyunLSQYnExERuXtZZgvPLttnLfLuzo7MGRSsIi9SyOWqzK9bt47x48fz6quvEhMTw4ULF4iKiuLSpUu8/PLLvPzyy/zxxx/5lVVExDAuTg5M7duERhWKA3ApIY1+s0KJjk81NpiIiMhdMJstjFsRxqp9kUD28W7mgKa6nUzEBuTqMvsePXpQvHhxpk2bdtPnhw8fTkJCAosXL86zgAVNl+SJyO1cTU6n+7StHL2YPVVPUBkvlj7ZAh93Z4OTSVGmY9fN6ecicnsWi4U3Vh5k3tbTADg5mJjevwn/CSptcDKRoi1fLrMPDQ2lX79+t3y+X79+bNu2LTe7FBGxKcU9XPhqcAjlirsDcCQqgWFf7SQ1I8vgZCK5s379eqMjiIiBLBYL768JtxZ5BxN82rORiryIDclVmb948SKVKlW65fOVK1cmKirqbjOJiBRqZXzcWDA0hJKeLgCEnozlv0v3kmXO9XiiIobp2LEjVatW5Z133uHs2bNGxxGRAjb59wimbjhuXZ74RAMerh9gYCIRya1clfnU1FRcXFxu+byzszPp6el3HUpEpLCr7OfJnEHBeLg4AvDTgSjeXHXwH9NlihRWkZGRjB07lhUrVlC5cmU6dOjAsmXLdBwXKQJmbjzBR2uPWpf/91hdHm9S/jZfISKFUa7umXdwcOCdd96hWLGbzzWZkJDA//3f/5GVZbuXm+r+OhHJjfXh0Qydt5PMa2flX+hQk1HtqhmcSoqauz127d27l9mzZ7N48WLMZjN9+vRhyJAhNGjQIB/SFhwd00X+aXHoGcat2G9dHv9wLYa2tv3ppUXsSU6PX7kq85UqVcJkMv3rdidPnszpLgsdHfhFJLe+2XWO577eZ12e+ER9ujUNNDCRFDV5ceyKjIxk+vTpvPfeezg5OZGamkqLFi2YOnUqderUyePEBUPHdJEbfb/3PM8s3cv1d///vb8GY++vbmwoEfmHnB6/nHKz01OnTt1tLhERu/N4k/JEJ6Tx/pojALy8Yj9+Xq60q+lvcDKR28vIyOD7779n9uzZrF27lqZNm/L555/Tq1cvYmNjeemll+jWrRuHDh0yOqqI3KW1hy7y7LJ91iI/vE0Vnr5PV5KJ2LJc3TP/b86ePcvgwYPzcpciIjZhxL1VGNiyEgBZZgsjF+xm79mrhmYSuZ0xY8YQEBDAiBEjqFGjBnv27GHr1q0MHToUT09PAgMDee+99zhy5IjRUUXkLm2OiGHUot3WgVp7NavAuAeDcnTFrYgUXnla5mNjY5k3b15e7lJExCaYTCZee6Q2D9fLHgk4JSOLwXN3cDImyeBkIjd36NAhJk+eTGRkJJ988gl169b9xzZly5Zl3bp1BqQTkbyy63QsQ+ftJD3TDEDnhmV5p0tdFXkRO5CnZV5EpChzdDAxqXsDQir7AhCblE7/2duJTkg1OJnIP73++ut069btH7PUZGZm8scffwDg5OTEvffea0Q8EckDByPjGDhnBykZ2YNTP1C7NB92a4Cjg4q8iD1QmRcRyUNuzo5M79+UoDJeAJyNTWHQnB0kpmUanEzkRu3atSM2NvYf6+Pi4mjXrp0BiUQkL0VEJ9J/VigJqdnHn1bV/JjcqxHOjnr7L2Iv9NcsIpLHfNydmTuoGeWKuwNwMDKeEfN3WS9xFCkMLBbLTS+zvXz5Mp6engYkEpG8cjY2mb4zt3M5KR2AJhVLML1/E9ycHQ1OJiJ5KVej2Xft2vW2z1+9evVusoiI2I0yPm7MGxzM41O2EpeSwaaIGF5cvo+PujfEQZc3ioGuH8tNJhMDBw7E1dXV+lxWVhZhYWG0bNnSqHgicpcuxqfSZ+Z2ouKzb/GqU9ab2QOD8XDJ1dt+EbEBufqr9vHx+dfn+/fvf1eBRETsRTV/L2YPbErvGdtJyzTz3d5IypVw54UOQUZHkyLs+rHcYrHg5eWFu7u79TkXFxeaN2/OsGHDjIonInchNimdPjO3cyY2GYBq/sX4anAzfNydDU4mIvkhV2V+zpw5+ZVDRMQuNanoy+RejXhywS4sFvhi3XHKFfegd0gFo6NJEXX9WF6pUiWef/55XVIvYifiUzPoP3s7EdGJAAT6urNgSAgli7n+y1eKiK3SPfMiIvmsfZ0yvP5Ibevya98fYF14tIGJRLJHs1eRF7EPyemZDJ6zgwPn4wEo7e3KwiHNKePjZnAyEclPuTozP3jw4BxtN3v27DsKIyJirwbeU5lzV1KYuekkWWYLoxfuZtmIFtQpe/vbl0TyUuPGjfntt98oUaIEjRo1uu0807t37y7AZCJyp9Iys3hy/i52nr4CgK+nCwuHhlChpIfByUQkv+WqzM+dO5eKFSvSqFEjLBZLfmUSEbFLrzxUi/NXU/jpQBRJ6VkMnruDb0feQ9ni7v/+xSJ5oHPnztYB77p06WJsGBG5axlZZsYs2sPGYzEAeLk58dXgZlTz9zI4mYgUBJMlF6185MiRLFmyhAoVKjB48GD69u2Lr69vfuYrcPHx8fj4+BAXF4e3t7fRcUTEzqRmZNF7xjZ2n7kKQM3SXnz9VAu83TQ4kdw5HbtuTj8XsWdms4Vnl+3lu72RALg7OzJ/SDOaVrKv9+YiRVFOj1+5umf+yy+/5MKFC7z00kusWrWKwMBAunfvzs8//6wz9SIiOeDm7MiM/k2peO3yx/CLCYxcsJuMLM1BLwXr7NmznDt3zrocGhrKM888w/Tp0w1MJSI5YbFYeO37A9Yi7+LowIz+TVXkRYqYXA+A5+rqSq9evVi7di2HDh2iTp06jBw5kooVK5KYmJgfGUVE7ErJYq7MHdSMEh7ZZ+M3RcQwbsV+fSgqBap3796sW7cOgKioKO6//35CQ0N55ZVXeOuttwxOJyK3YrFYeO+nIyzcfgYARwcTn/duRKvqfgYnE5GCdlej2ZtMJkwmExaLBbNZZ5VERHKqsp8nMwc0xcUp+2V4+a5zfPZbhMGppCg5cOAAzZo1A2DZsmXUq1ePLVu2sGjRIubOnWtsOBG5pc9/j2DaHycAMJngo+4NaF+njMGpRMQIuS7zaWlpLF68mAceeICaNWuyf/9+Pv/8c86cOUOxYsXyI6OIiF1qUtGXT3o05PqA4h//epTlu87d/otE8khGRoZ1MLxff/2VTp06ARAUFMSFCxeMjCYitzBr00kmrT1qXf5fl3p0bljOwEQiYqRclfmRI0cSEBDA+++/zyOPPMK5c+f4+uuveeihh3Bw0JT1IiK59VC9AF55sJZ1+eVvwtgcEWNgIikq6tSpw9SpU9m4cSNr166lY8eOAERGRlKyZEmD04nI3y3dcYa3fzhkXX71oVr0DqlgYCIRMVquRrN3cHCgQoUK/zo37YoVK/IknBE08q2IFDSLxcLrKw/y1dbTAHi5OvHNyJbUKK2phSRn7uTYtX79eh577DHi4+MZMGAAs2fPBuCVV17hyJEjNn0sv07HdLEXq/ZF8vSSPVx/1/70fdV59oEaxoYSkXyT0+NXruaZ79+//21LvIiI5J7JZOL1R+sQeTWFXw9Hk5CWyeC5O/hu1D34FXM1Op7YqbZt2xITE0N8fDwlSpSwrh8+fDgeHh4GJhORv/rt8EX+u3SvtcgPaVWZ/95f3dhQIlIo5OrMfFGgT/FFxCjJ6Zn0mLaN/efjAGhcoTiLhjXHzdnR4GRS2OnYdXP6uYit23I8hoFzdpCemT3QdM/gQCZ0raeTayJ2Ll/mmRcRkfzj4eLEzAFNKePtBsDuM1d5cXmYpqyTfHHx4kX69etH2bJlcXJywtHR8YaHiBhr95krDJ2301rkH21Qlv89piIvIn9SmRcRKURKe7sxc0BT3K+djV+5L5JPfztmcCqxRwMHDmT37t289tprLF++nBUrVtzwyI0JEyYQHByMl5cX/v7+dOnShfDw8H9sd/jwYTp16oSPjw9eXl40b96cM2fO3Hbfn3zyCTVr1sTd3Z3AwED++9//kpqamqt8IrbmUGQ8A2eHkpyeBcD9tfz5qHsDHB1U5EXkT7m6Z15ERPJf3XI+fNarEcPn78RigU9+PUZlP09NPyR5atOmTWzcuJGGDRve9b42bNjAqFGjCA4OJjMzk1dffZX27dtz6NAhPD09ATh+/DitWrViyJAhvPnmm/j4+HD48GHc3Nxuud+FCxfy8ssvM3v2bFq2bMnRo0cZOHAgAB9//PFd5xYpjI5fSqT/7O3Ep2YC0LJqST7v3RhnR52DE5EbqcyLiBRCD9QuzSsP1uJ/Px4G4IXlYZQv4U6Tir4GJxN7ERgYmGe3cKxZs+aG5Tlz5uDv78+uXbto06YNAK+++ioPPfQQH3zwgXW7KlWq3Ha/W7du5Z577qF3794AVKpUiV69ehEaGponuUUKm7OxyfSduZ2YxHQAGlUozoz+TTV2iojclD7iExEppIa2rkzP4EAA0jPNDP9qF2djkw1OJfbik08+4eWXX+bUqVN5vu+4uOxBHH19sz98MpvNrF69mho1atChQwf8/f0JCQnhu+++u+1+WrVqxa5du6zl/cSJE/z44488/PDDN90+LS2N+Pj4Gx4itiI6PpW+s7ZzIS77NpJaAd7MHdgMT1edexORm1OZFxEppEwmE293qUvLqiUBuJyUzuC5O4hPzTA4mdiDHj16sH79eqpWrYqXlxe+vr43PO6UxWLh2WefpVWrVtStWxeA6OhoEhMTee+99+jYsSO//PILjz32GF27dmXDhg233FfPnj15++23adWqFc7OzlStWpV27drx8ssv33T7CRMm4OPjY30EBgbe8fchUpCuJKXTd9Z2Tl/O/sC2SilP5g9pho+Hs8HJRKQw09R0f6NpbESksIlLzuCxKZs5cSkJgNbV/ZgzMBgn3T8p19zJsWvevHm3fX7AgAF3lGXUqFGsXr2aTZs2Ub58eQAiIyMpV64cvXr1YtGiRdZtO3XqhKenJ4sXL77pvtavX0/Pnj155513CAkJISIigrFjxzJs2DBee+21f2yflpZGWlqadTk+Pp7AwEAd06VQS0jNoM/M7YSdy76ipVxxd5Y/1YIAH3eDk4mIUXJ6XNd1OyIihZyPhzOzBwTT5cvNXE3OYOOxGN5cdYi3OtfRFEVyx+60rN/OmDFjWLlyJX/88Ye1yAP4+fnh5ORE7dq1b9i+Vq1abNq06Zb7e+211+jXrx9Dhw4FoF69eiQlJTF8+HBeffVVHBxu/EDL1dUVV1fXPPyORPJXakYWQ+fttBb5Ul6uLBoWoiIvIjmi0zoiIjagkp8n0/o2wdkxu7zP33aauVtOGRtKbN7x48cZP348vXr1Ijo6GsgezO7gwYO52o/FYmH06NGsWLGC33//ncqVK9/wvIuLC8HBwf+Yru7o0aNUrFjxlvtNTk7+R2F3dHTEYrHk2eB9IkbJyDIzauFutp+MBaC4hzMLhoRQsaSnwclExFaozIuI2IiQKiWZ0LW+dfntHw6xLjzawERiyzZs2EC9evXYvn07K1asIDExEYCwsDBef/31XO1r1KhRLFiwgEWLFuHl5UVUVBRRUVGkpKRYt3nhhRdYunQpM2bMICIigs8//5xVq1YxcuRI6zb9+/dn3Lhx1uVHH32UKVOmsGTJEk6ePMnatWt57bXX6NSpE46OGt1bbFeW2cJzy/bx25Hs13BPF0fmDmpGzTJeBicTEVuiMi8iYkOeaFKekW2rAmC2wNOL9hARnWBwKrFFL7/8Mu+88w5r167FxcXFur5du3Zs3bo1V/uaMmUKcXFxtG3bloCAAOtj6dKl1m0ee+wxpk6dygcffEC9evWYOXMm33zzDa1atbJuc+bMGS5cuGBdHj9+PM899xzjx4+ndu3aDBkyhA4dOjBt2rS7+M5FjGWxWPi/7w+wcl8kAC5ODswY0JSGgcWNDSYiNkcD4P2NBsATkcLObLYwcuFu1hyMAqBSSQ++H9VKox4XYXdy7CpWrBj79++ncuXKeHl5sW/fPqpUqcKpU6cICgoiNTU1n1PnPx3TpTD6YM0Rvlx/HABHBxNT+zbhgdqlDU4lIoVJTo9fOjMvImJjHBxMTOregKBrl2OeupzMqEW7ycwyG5xMbEnx4sVvOAt+3Z49eyhXrpwBiUTs35T1x61F3mSCSd0aqMiLyB1TmRcRsUGerk7MHNCUkp7Zl0dviojhndWHDU4ltqR379689NJLREVFYTKZMJvNbN68meeff57+/fsbHU/E7izcfpr31xyxLr/VqQ5dGumDMxG5cyrzIiI2qnwJD6b8ZYT7uVtOsST0jMGpxFb873//o0KFCpQrV47ExERq165NmzZtaNmyJePHjzc6nohd+X7vecZ/d8C6/EKHmvRrUcm4QCJiF3TP/N/o/joRsTVLQs/w8or9ADg7mlg0rDnBlXwNTiUF6W6OXcePH2fPnj2YzWYaNWpE9erV8yllwdMxXQqD349cZPhXu8g0Z7/lfvLeKrzcMQiTyWRwMhEprHJ6/HIqwEwiIpIPejarwJGoBOZuOUVGloUR83fx/eh7KF/Cw+hoYgOqVq1K1apVjY4hYpe2nbjMUwt2W4t8r2YVVORFJM+ozIuI2IHxD9ciIjqRTRExXE5KZ+i8nXzzVEs8XfUyL3969tlnc7ztRx99lI9JROxf2LmrDJ23k7TM7MFJH6kfwDtd6qrIi0ie0bs8ERE74OTowOe9G9Hli82cupzMkagEnlu2jy/7NMbBQW8cJduePXtuWN61axdZWVnUrFkTgKNHj+Lo6EiTJk2MiCdiN45dTGDA7FAS0zIBaFezFB91b4ijXo9FJA9pADwRETtR3MOFmQOC8bp2Nn7NwSg++e2YwamkMFm3bp318eijj9K2bVvOnTvH7t272b17N2fPnqVdu3Y8/PDDRkcVsVlnY5PpO2s7V5IzAGhW2Zcv+zTBxUlvu0Ukb+lVRUTEjlTzL8ZnvRtx/eTPZ78dY3XYP+cSF5k0aRITJkygRIkS1nUlSpTgnXfeYdKkSQYmE7Fd0fGp9J21nYvxaQDUK+fDrAFNcXdxNDiZiNgjlXkRETvTrqY/4x6sZV1+7uu9HDgfZ2AiKYzi4+O5ePHiP9ZHR0eTkJBgQCIR23Y1OZ1+s0I5fTkZyP5wdd7gZni5ORucTETslcq8iIgdGtq6Ml0blwMgNcPMk/N3EZuUbnAqKUwee+wxBg0axPLlyzl37hznzp1j+fLlDBkyhK5duxodT8SmJKZlMmDODsIvZn8QVr6EOwuGhODr6WJwMhGxZyrzIiJ2yGQy8e5j9WgYWByA81dTGL1oN5lZZmODSaExdepUHn74Yfr27UvFihWpWLEiffr04cEHH+TLL780Op6IzUjNyGL4VzvZd/YqAH7FXFkwJIQyPm7GBhMRu6cyLyJip9ycHZnatwmlvFwB2HL8Mu/9dMTgVFJYeHh48OWXX3L58mX27NnD7t27iY2N5csvv8TT09PoeCI2ITPLzJjFe9hy/DIAPu7OLBjajEp++hsSkfynMi8iYsfK+LgxpU9jnK6NiDdz00m+33ve4FRSmHh6elK/fn0aNGigEi+SC2azhReXh7H2UPbYEx4ujswZFExQGW+Dk4lIUaEyLyJi55pW8uX1TnWsyy99E6YB8URE7oLFYuHNVQdZsSf7w1EXRwdm9G9K4wol/uUrRUTyjsq8iEgR0DekAj2aBgIaEE9E5G59tPYo87aeBsDRwcTk3o24p5qfwalEpKhRmRcRKQJMJhNvdamjAfFERO7SjD9OMPn3COvyxCfq06FOGQMTiUhRpTIvIlJEuDplD4jnV+zPAfHeX6MB8UREcmpJ6Bn+9+Nh6/KbnerQtXF5AxOJSFGmMi8iUoSU8XFjSt8/B8SbsVED4omI5MQPYZGM+3a/dfm5B2owoGUl4wKJSJGnMi8iUsQEV/Ll9UdrW5df+iaMg5EaEE9E5FbWhUfz36V7sViyl4e1rszo/1QzNpSIFHkq8yIiRVDf5hXp3jT70tDUDDPDv9KAeCIiNxN6MpanFuwiIyu7yfdoGsgrD9XCZDIZnExEijqVeRGRIshkMvFW57o0+MuAeGMWa0A8EZG/OnA+jiFzd5Cakf3a+HC9AN7tWk9FXkQKBZV5EZEiys3Zkal9G+NXzAWAzRGX+fCXowanEhEpHCKiE+k/O5SEtEwA7q1Rio97NMTRQUVeRAoHlXkRkSIswMedL/s0sQ6IN3XDcX4+GGVwKhERY527kky/Wduttx8FVyrB1L5NcHHSW2cRKTz0iiQiUsQ1q+zLKw/Vsi4/v2wfJ2OSDEwkImKcSwlp9J25nQtxqQDUKevNrIHBuLs4GpxMRORGKvMiIsKgeyrxSP0AABLSMhkxfxfJ6ZkGpxIRKVhxKRn0nx3KqcvJAFQp5cm8wc3wdnM2OJmIyD+pzIuICCaTifcfr081/2IAhF9M4NVvD2C5Pg+TiIidS0nPYsjcHRy+EA9AueLuLBgSgl8xV4OTiYjcnMq8iIgA4OnqxNS+TfC8dinpt3vOs2DbaYNTiYjkv/RMM08t3MXO01cAKOnpwvwhzShb3N3gZCIit6YyLyIiVtX8i/HBEw2sy2/9cIjdZ64YmEhEJH+ZzRae/3of68MvAeDl6sS8wc2oUqqYwclERG5PZV5ERG7wcP0AhrSqDEBGloVRC3dzOTHN4FQiInnPYrHw+sqDrNwXCYCrkwMzBzSlbjkfg5OJiPw7lXkREfmHlx8MIrhSCQAuxKUydslessy6f15E7MtHa48y/9rtRI4OJr7s05iQKiUNTiUikjMq8yIi8g/Ojg580buxdeCnTRExfLQ23OBUIiJ5Z9amk0z+PcK6PKlbA+6rVdrARCIiuaMyLyIiN+Xv7cYXvRvh6GAC4It1x/n10EWDU4mI3L3lu87x9g+HrMtvPFqbLo3KGZhIRCT3VOZFROSWQqqU5OWOQdbl/y7by+nLSQYmEhG5O78cjOKlb8Ksy8/cX52B91Q2MJGIyJ1RmRcRkdsa2royD9UrA0BCaiYjFuwmNSPL4FQiIrm39fhlRi/eYx0DZGDLSoy9r7rBqURE7ozKvIiI3JbJZOKDJxpQpZQnAIcvxPPmqkP/8lUiIoVL2LmrDPtqJ+mZZgAea1SO/3ukNiaTyeBkIiJ3RmVeRET+VTFXJ6b0aYKbc/ZhY3HoGb7fe97gVCIiORMRncjAOTtITMsE4P5a/nzwRH0cHFTkRcR2qcyLiEiO1CzjxVud61qXx63YT0R0ooGJRET+3fmrKfSbtZ3YpHQAmlX25fPejXF21NtgEbFtehUTEZEc6940kMcblwcgOT2LUQt3k5Ku++dFpHCKSUyj38ztXIhLBaBOWW9mDmiKm7OjwclERO6eyryIiOTK213qUN2/GADhFxN4Y+VBgxOJiPxTQmoGA+eEciImewaOKn6ezBvcDG83Z4OTiYjkDZV5ERHJFQ8XJ77s0xj3a2e2lu48y4rd5wxOJSLyp9SMLIbO28mB8/EABPi48dWQZvgVczU4mYhI3lGZFxGRXKte2ov/Pfbn/fOvfnuAYxcTDEwkIpItI8vM6EW72X4yFoASHs7MH9KM8iU8DE4mIpK3VOZFROSOdG1cnh5NAwFIychi5MLdJKdnGpxKRIoys9nCS8vD+PVwNACeLo7MHdSMav5eBicTEcl7KvMiInLH3uxch6Ay2W+Sj0Un8tp3un9eRIxhsVh4e/UhVuzJnjbTxdGBGf2b0iCwuLHBRETyicq8iIjcMTdnR77o0xhPl+z757/ZfY5lO88anEpEiqLJv0cwZ/MpABxMMLl3I1pW8zM2lIhIPlKZFxGRu1K1VDHe7VrPuvx/3x8gPEr3z4tIwZm35RQfrT1qXX7v8fp0qFPGwEQiIvlPZV5ERO5a54bl6B1SAYDUDDNPLdxFUprunxeR/Pf93vO8/pcpMsc/XIvu18bzEBGxZzZR5k+dOsWQIUOoXLky7u7uVK1alddff5309PQbttuxYwf33XcfxYsXp0SJErRv3569e/caE1pEpIj5v0dqUzvAG4ATl5JueHMtIpIffj9ykeeW7bMuj2pXlaGtqxiYSESk4NhEmT9y5Ahms5lp06Zx8OBBPv74Y6ZOncorr7xi3SYhIYEOHTpQoUIFtm/fzqZNm/D29qZDhw5kZGQYmF5EpGj4+/3zy3ed4/u95w1OJSL2KvRkLE8t2E2m2QJAn5AKPN++psGpREQKjslisViMDnEnJk6cyJQpUzhx4gQAO3fuJDg4mDNnzhAYmH1p1f79+6lfvz4RERFUrVo1R/uNj4/Hx8eHuLg4vL298y2/iIi9WrH7HM9eO1NWzNWJH8a0opKfp8Gp7JuOXTenn4v9OhgZR89p20i4djvPI/UD+LRnIxwdTAYnExG5ezk9ftnEmfmbiYuLw9fX17pcs2ZN/Pz8mDVrFunp6aSkpDBr1izq1KlDxYoVb7mftLQ04uPjb3iIiMid69q4PF0blwMgMS2Tp5fsIT3TbHAqEbEXJ2OSGDA71Frk761Rio+6N1SRF5EixybL/PHjx5k8eTIjRoywrvPy8mL9+vUsWLAAd3d3ihUrxs8//8yPP/6Ik5PTLfc1YcIEfHx8rI/rZ/VFROTOvd25LpWvnY0POxfHxJ+PGJxIROzBhbgU+s7cTkxi9rhJTSqWYErfxrg42eRbWhGRu2LoK98bb7yByWS67WPnzp03fE1kZCQdO3akW7duDB061Lo+JSWFwYMHc88997Bt2zY2b95MnTp1eOihh0hJSbllhnHjxhEXF2d9nD2r+ZFFRO6Wp6sTk3s1wtkx+0zZjI0nWRcebXAqEbFlV5LS6T8rlPNXs9/XBZXxYvaAYDxcbn3SRkTEnhl6z3xMTAwxMTG33aZSpUq4ubkB2UW+Xbt2hISEMHfuXBwc/vwsYtasWbzyyitcuHDBuj49PZ0SJUowa9YsevbsmaNMur9ORCTvzN50krd+OARASU8XfhrbGn9vN4NT2R8du25OPxf7kZiWSZ+Z29l39ioAFXw9WD6ihV5PRMQu5fT4ZehHmX5+fvj5+eVo2/Pnz9OuXTuaNGnCnDlzbijyAMnJyTg4OGAy/Xm/1PVls1n3aoqIGGHQPZXYHBHDb0eiuZyUzn+X7WX+4BAcdG+riORQWmYWT87faS3y/l6uLBgSoiIvIkWeTdxgFBkZSdu2bQkMDOTDDz/k0qVLREVFERUVZd3mgQce4MqVK4waNYrDhw9z8OBBBg0ahJOTE+3atTMwvYhI0WUymZjYrQGlvV0B2BxxmSkbjhucSkRsRWaWmbGL97I54jIAPu7OzB8SQoWSHgYnExExnk2U+V9++YWIiAh+//13ypcvT0BAgPVxXVBQEKtWrSIsLIwWLVrQunVrIiMjWbNmzQ3biYhIwfL1dOGTHo24fuHUR2uPsut0rLGhRKTQs1gsvPLtftYczD554+7syOyBwdQs42VwMhGRwsFm55nPL7q/TkQkf3z0Szif/R4BQLni7vw4tjU+7s4Gp7IPOnbdnH4utm3Cj4eZ9scJAJwdTcwaEEybGqUMTiUikv/sfp55ERGxLU/fV53gSiUAOH81hXErwtDnySJyM9M2HLcWeZMJPunRSEVeRORvVOZFRKRAODk68EnPRtaz8T/uj2LJDk0HKiI3WrbzLBN+OmJdfqdLXR6ur1smRUT+TmVeREQKTLni7rz/eH3r8lurDnHiUqKBiUSkMPnlYBQvfxNmXX6+fQ36hFQ0MJGISOGlMi8iIgWqY90y9A6pAEBKRhZjl+wlPVNTiIoUddtOXGb04j2Yr919M7BlJUa1q2ZsKBGRQkxlXkRECtz4h2tRxc8TgP3n4/jk16MGJ5K7MWHCBIKDg/Hy8sLf358uXboQHh7+j+0OHz5Mp06d8PHxwcvLi+bNm3PmzJnb7vvq1auMGjWKgIAA3NzcqFWrFj/++GN+fStikAPn4xg2b6f1g70uDcvyf4/UxnR9GgwREfkHlXkRESlwHi5OfNqzEU4O2W/Up2w4zrYTlw1OJXdqw4YNjBo1im3btrF27VoyMzNp3749SUlJ1m2OHz9Oq1atCAoKYv369ezbt4/XXnsNNze3W+43PT2dBx54gFOnTrF8+XLCw8OZMWMG5cqVK4hvSwrIqZgkBs4JJSEtE4C2NUsxsVsDHBxU5EVEbkdT0/2NprERESk4U9Yf5/012QNdlfVx46exbfDx0HR1uVXYjl2XLl3C39+fDRs20KZNGwB69uyJs7Mz8+fPz/F+pk6dysSJEzly5AjOzrn/vShsPxf5p+j4VB6fuoWzsSkANK5QnAVDQ/BwcTI4mYiIcTQ1nYiIFHrD21SheRVfACLjUnn1u/2ars4OxMXFAeDrm/3/1mw2s3r1amrUqEGHDh3w9/cnJCSE77777rb7WblyJS1atGDUqFGULl2aunXr8u6775KVlXXT7dPS0oiPj7/hIYVXXEoG/WeHWot8jdLFmD0wWEVeRCSHVOZFRMQwjg4mPureEG+37DfvP4Rd4Ns95w1OJXfDYrHw7LPP0qpVK+rWrQtAdHQ0iYmJvPfee3Ts2JFffvmFxx57jK5du7Jhw4Zb7uvEiRMsX76crKwsfvzxR8aPH8+kSZP43//+d9PtJ0yYgI+Pj/URGBiYL9+j3L2U9CyGztvBkagEIHumi68Gh1Dcw8XgZCIitkOX2f+NLskTESl4P4RFMnrRHgCKuTrx49OtqVDSw+BUtqMwHbtGjRrF6tWr2bRpE+XLlwcgMjKScuXK0atXLxYtWmTdtlOnTnh6erJ48eKb7qtGjRqkpqZy8uRJHB0dAfjoo4+YOHEiFy5c+Mf2aWlppKWlWZfj4+MJDAwsFD8X+VNGlpkR83fx25FoAEp6urD8qZZUvjYopohIUafL7EVExGY8Ur8sjzfOLn6JaZk8s3QPmVmars7WjBkzhpUrV7Ju3TprkQfw8/PDycmJ2rVr37B9rVq1bjuafUBAADVq1LAW+etfExUVRXp6+j+2d3V1xdvb+4aHFC5ms4WXvgmzFvlirk7MG9xMRV5E5A6ozIuISKHwRqfaBPq6A7D7zFW+WHfc4ESSUxaLhdGjR7NixQp+//13KleufMPzLi4uBAcH/2O6uqNHj1KxYsVb7veee+4hIiICs9l8w9cEBATg4qLLsW2NxWLh3R8Ps2J39q00Lo4OTO/fhLrlfAxOJiJim1TmRUSkUPByc+aTHo1wvDYd1We/H2PX6SsGp5KcGDVqFAsWLGDRokV4eXkRFRVFVFQUKSkp1m1eeOEFli5dyowZM4iIiODzzz9n1apVjBw50rpN//79GTdunHX5qaee4vLly4wdO5ajR4+yevVq3n33XUaNGlWg35/kjSkbjjNz00kAHEzwWa9GtKzqZ3AqERHbpTIvIiKFRpOKJRjzn2oAZJktPLN0DwmpGQankn8zZcoU4uLiaNu2LQEBAdbH0qVLrds89thjTJ06lQ8++IB69eoxc+ZMvvnmG1q1amXd5syZMzfcCx8YGMgvv/zCjh07qF+/Pk8//TRjx47l5ZdfLtDvT+7ektAzfLDmzysz3n2sHh3rljEwkYiI7dMAeH9TmAYREhEpijKzzHSftpXdZ64C0KNpIO8/Ud/YUIWcjl03p59L4bDmwAVGLtyN+do7zhc71mRk22rGhhIRKcQ0AJ6IiNgkJ0cHPunRCE+X7EHPlu48y9pDFw1OJSJ3YsvxGJ5evNda5Ie2qsxT91Y1NpSIiJ1QmRcRkUKnQkkP/u/RP0c+H7cijMuJabf5ChEpbA6cj2P4V7tIvzYzRdfG5XjloVqYTCaDk4mI2AeVeRERKZS6Nw3k/lr+AMQkpvPKt/vRnWEituHEpUQGzA4lMS0TgPuC/Hn/8fo4OKjIi4jkFZV5EREplEwmExO61sfXM3sKsp8PXrROaSUihdfF+FT6zQrlclI6AMGVSvBFn8Y4O+ptp4hIXtKrqoiIFFqlvFz5X5e61uU3Vh7k/NWU23yFiBgpLjmD/rNCrX+nQWW8mDkgGDdnR4OTiYjYH5V5EREp1B6sF0DXRuUASEjL5IWv92E263J7kcImJT2LwfN2EH4xAYBAX3e+GtwMH3dng5OJiNgnlXkRESn0Xu9Uh7I+bgBsOX6ZuVtOGRtIRG6QkWXmqYW72HX6CgB+xVyZPzgEf283g5OJiNgvlXkRESn0fNydmditgXX5/TVHiIhOMDCRiFxnNlt44et9rA+/BICXqxPzBgdTyc/T4GQiIvZNZV5ERGzCPdX8GHRPJQDSMs08u2wfGdemvBIRY1gsFt5efYjv9kYC4OLkwIwBTalT1sfgZCIi9k9lXkREbMZLHYOoWir7bF/YuTg+/z3C4EQiRdsX6yKYs/kUAA4m+LxXI5pXKWlsKBGRIkJlXkREbIabsyMf92iI47W5qj9fF8G+s1eNDSVSRC3cfpoPfzlqXX7v8fq0r1PGwEQiIkWLyryIiNiU+uWLM+Y/1QDIMlv477K9pGZkGZxKpGj5cf8Fxn93wLo87sEgujcNNDCRiEjRozIvIiI2Z1S7atQvn31P7olLSUz6JdzgRCJFx+aIGJ5ZshfLtRkin2xThSfvrWpsKBGRIkhlXkREbI6zowMfdW+Ii1P2YWzmppPsOh1rcCoR+xd27irDv9pJ+rXBJ7s1Kc/LDwYZnEpEpGhSmRcREZtUzb8Yzz1QAwCLBV74OkyX24vko4joRAbO2UFSevbf2QO1SzOhaz1MJpPByUREiiaVeRERsVlDW1ehUYXiAJyISeLDn3W5vUh+uBCXwoDZocQmpQPQrLIvk3s1wslRbyVFRIyiV2AREbFZjg4mJj7RwHq5/azNJ9lxSpfbi+SlK0np9J8VyvmrKQDUDvBm5oCmuDk7GpxMRKRoU5kXERGbVs2/GM+3/+vl9vtISdfl9iJ5ITk9k8HzdnAsOhGAiiU9mDe4Gd5uzgYnExERlXkREbF5Q1pVofG1y+1PXU5moi63F7lrGVlmRi7czZ4zVwEo5eXK/MEhlPJyNTaYiIgAKvMiImIHHB1MTOzWANdrl9vP2XKS0JO63F7kTpnNFl5aHsb68EsAeLk68dXgZlQo6WFwMhERuU5lXkRE7ELVUsV4vn1NIPty+xeX63J7kTv13pojrNhzHgAXJwdmDmhKrQBvg1OJiMhfqcyLiIjdGNyqMk0qlgCyL7f/4OcjBicSsT3T/zjO9D9OAOBggsm9GhFSpaTBqURE5O9U5kVExG5kj25f33q5/dwtp3S5vUgufLPrHO/++OeHYP97rB4d6pQxMJGIiNyKyryIiNiVKqWK8UKHPy+3f2H5PpLTMw1OJVL4rTsSzYvfhFmXn3ugBr2aVTAwkYiI3I7KvIiI2J1B91Sm6bXL7U9fTuaDNRrdXuR2dp2+wlMLd5FltgAwoEVFRv+nmsGpRETkdlTmRUTE7jg6mPjgL5fbz9t6il2ndbm9yM0cu5jA4Lk7SM0wA/Bw/QD+79E6mEwmg5OJiMjtqMyLiIhd+vvl9i8uDyM1Q6Pbi/xV5NUU+s8OJS4lA4B7qpXko+4NcHRQkRcRKexU5kVExG4NuqcyDQKLA3D8UhJfrIswNpBIIXIlKZ3+s0O5EJcKQL1yPkzr1xRXJ0eDk4mISE6ozIuIiN1ydDDx/uP1cLp2lnHK+uMciow3OJWI8ZLTMxk8bwcR0YkAVCrpwZxBwRRzdTI4mYiI5JTKvIiI2LWgMt6MbJc9kFem2cJL34SRmWU2OJWIcTKyzIxauJs9Z64CUMrLlflDQvAr5mpsMBERyRWVeRERsXuj2lWlun8xAPafj2PWppMGJxIxhsWS/YHWuvBLAHi5OjFvUDMCfT0MTiYiIrmlMi8iInbP1cmR95+oz/XBuT9ae5STMUnGhhIxwHs/HWHF7vMAuDg5MGNAU2qX9TY4lYiI3AmVeRERKRIaVyjBoJaVAUjLNPPyN2GYr82pLVIUzPjjBNP+OAGAgwk+69mQ5lVKGpxKRETulMq8iIgUGc93qEH5Eu4AbD8Zy+IdZwxOJFIwVuw+x/9+PGxdfqdLPTrWDTAwkYiI3C2VeRERKTI8XJx4r2t96/KEH49wIS7FwEQi+W9deDQvLg+zLj/7QA16h1QwMJGIiOQFlXkRESlSWlX3o1uT8gAkpmUy/tsDWCy63F7s0+4zVxi5YDeZ124p6d+iImP+U83gVCIikhdU5kVEpMgZ/3BtSnllT8P125FoVoVdMDiRSN6LiE5g8NwdpGRkAfBwvQBef7QOpusjQYqIiE1TmRcRkSLHx8OZtzvXsS6/sfIgsUnpBiYSyVuRV1PoPyuUq8kZALSsWpKPejTA0UFFXkTEXqjMi4hIkdSxbgAP1i0DQGxSOu+sPmRwIpG8cTU5nQGzQ4mMSwWgbjlvpvVrgquTo8HJREQkL6nMi4hIkfVm5zp4uTkBsGL3ebZExBicSOTupKRnMWTeTo5FJwJQsaQHcwY2w8vN2eBkIiKS11TmRUSkyPL3cuPlB4Osy69+d4DUa/cXi9iajCwzoxbtZtfpKwD4FXNl/uAQ6/gQIiJiX1TmRUSkSOsVXIEmFUsAcDImiS/WRRicSCT3LBYLL3+zn9+PRAPg5erEvMHBVCjpYXAyERHJLyrzIiJSpDk4mHj3sXo4XRsYbOqG4xy7mGBwKpHceW/NEb7ZfQ4AFycHpvdvSp2yPganEhGR/KQyLyIiRV7NMl48eW8VADKyLLzy7X7MZs09L7Zh5sYTTNtwAgCTCT7r2ZAWVUsanEpERPKbyryIiAgw5j/VqXTtkuQdp66wdOdZgxOJ/Ltv95zjndWHrcvvdKlLx7oBBiYSEZGCojIvIiICuDk78k6XetblCT8e5lJCmoGJRG5vXXg0L3wdZl3+7/016BNS0cBEIiJSkFTmRURErmlV3Y/HGpUDID41k7d/0NzzUjjtOXOFkQt2k3ntdpB+zSvy9H3VDE4lIiIFSWVeRETkL8Y/XIviHtlzcq/cF8n68GiDE4ncKCI6gUFzd5BybRrFh+qV4Y1OdTCZTAYnExGRgqQyLyIi8hcli7nyykO1rMuvfX+AlHTNPS+Fw4W4FPrPCuVqcgYALauW5OMeDXF0UJEXESlqVOZFRET+pluT8oRU9gXgbGwKn/52zOBEInA1OZ3+s0KJjEsFoE5Zb6b1a4Krk6PByURExAgq8yIiIn9jMpl4t2s9XByzD5MzNp7g8IV4g1NJUZaSnsWQeTs5Fp0IQMWSHswd1AwvN2eDk4mIiFFU5kVERG6iaqlijGxXFYAss+aeF+NkZpkZvWg3u05fAcCvmCtfDW5GKS9Xg5OJiIiRVOZFRERu4am2ValSyhOAPWeuskxzz0sBs1gsvPrtAX47kj0QYzFXJ+YOCqZiSU+Dk4mIiNFU5kVERG7B1cmRdzrXtS6/t+YIsUnpBiaSouajtUdZeu1DJBdHB6b3b0Ldcj4GpxIRkcJAZV5EROQ2Wlbzo3PDsgBcTc7ggzVHDE4kRcX8raeY/HsEACYTfNSjAS2r+hmcSkRECguVeRERkX/x6kO18HJ1AmDJjrPWe5dF8stP+y/wfysPWpdff6Q2j9Qva2AiEREpbFTmRURE/oW/txvPtq9hXR7/3QEys8wGJhJ7tu3EZcYu2Yvl2niLT7WtysB7KhsbSkRECh2VeRERkRzo17witQO8ATh8IZ6vtp42OJHYoyNR8Qz7aifp1z4serxxeV7sUNPgVCIiUhipzIuIiOSAk6MD7zz252B4H609SnR8qoGJxN6cu5LMgNmhJKRmAtCuZinee7weJpPJ4GQiIlIYqcyLiIjkUOMKJejVLBCAxLRM3ll92OBEYi+uJKUzYHYoF+PTAGgQWJwv+jTG2VFv1URE5OZ0hBAREcmFFzsEUcLDGYCV+yLZHBFjcCKxdSnpWQyZt4Pjl5IAqOLnyZyBwXi4OBmcTERECjOVeRERkVwo4enCyw8GWZdf+/4AaZlZBiYSW5aZZWb0ot3sPnMVgFJerswb3AxfTxdjg4mISKGnMi8iIpJL3ZoE0rhCcQBOXEpi5saTxgYy2IQJEwgODsbLywt/f3+6dOlCeHj4P7Y7fPgwnTp1wsfHBy8vL5o3b86ZM2dy9G8sWbIEk8lEly5d8ji9cSwWC69+e4DfjkQD4OXqxLxBzQj09TA4mYiI2AKVeRERkVxycDDxTpd6OFwbl2zy78c4G5tsbCgDbdiwgVGjRrFt2zbWrl1LZmYm7du3JykpybrN8ePHadWqFUFBQaxfv559+/bx2muv4ebm9q/7P336NM8//zytW7fOz2+jwH289ihLd54FwMXRgWn9m1C7rLfBqURExFaYLJbrs5gKQHx8PD4+PsTFxeHtrQOqiIjc2purDjJn8ykA7q9VmpkDmhqSo7Aduy5duoS/vz8bNmygTZs2APTs2RNnZ2fmz5+fq31lZWVx7733MmjQIDZu3MjVq1f57rvvcvS1he3n8lfzt53mte8OAGAyweRejXikflmDU4mISGGQ0+OXzsyLiIjcoWcfqIG/lysAvx6+yLrwaIMTFQ5xcXEA+Pr6AmA2m1m9ejU1atSgQ4cO+Pv7ExISkqNS/tZbb1GqVCmGDBnyr9v+f3t3Hh5VYb59/J4kZBIgBJMgyZCFsJRdlM0CAUO1iCC4FVFke0EtLTu+KQht4ZVCIlGwgLJVQKUsV1tW3x+WyKZIkZgQIVSI7EigkRaSECHr+f3BxdQIQpZJzpzJ93Nd88eczEzuJ5OcJ885Z87Jz89XTk5OqZs7+ij9gn6/Od15f8bjrRnkAQDlxjAPAEAFBfjV0vR+rZz3X9v6zxp/MjzDMDR58mTFxMSobdu2kqSsrCxdvXpVCQkJ6tOnj7Zv366nnnpKTz/9tPbs2fOjr/XZZ5/p3Xff1fLly8v0vePj4xUYGOi8RUREuKQmV9p/8t8avy5NN4+L/FVsU43oHm1uKACAJTHMAwBQCQPaO9S58T2SpFOX8pyH3ddUY8eO1aFDh7R27VrnspKSEknSE088oUmTJun+++/X1KlT9fjjj2vJkiW3fZ3c3FwNGTJEy5cvV0hISJm+96uvvqrs7Gzn7dy5c5UvyIWOXszRS+9/oYKiGz+PZzqE6zePtjA5FQDAqriAKQAAlWCz2TRzQBv1X7hXJYa0cMfXeuqBRmpY7+4ndvM048aN05YtW/TJJ58oPDzcuTwkJEQ+Pj5q3bp1qce3atVKe/fuve1rnThxQqdPn1b//v2dy25uFPDx8dGxY8fUtGnTUs+x2+2y2+2uKselzl+5puErDij3epEkKbZFAyU80042m83kZAAAq2LPPAAAldTGEajBD0ZKkvIKihX/P1+ZnKh6GYahsWPHasOGDdq5c6eio0sfNu7r66vOnTvfcrm6jIwMRUVF3fY1W7ZsqcOHDystLc15GzBggHr16qW0tDS3PIT+x1zOK9Cwdz/Xv3LyJUntI+rrnRc6qJY3/4YBACqOPfMAALjAKz9voQ8PXdCV7wq1KS1TL/w0Sp0bB5kdq1qMGTNGa9as0ebNmxUQEKCLFy9KkgIDA+Xv7y9JiouL06BBg9SzZ0/16tVLH330kbZu3ardu3c7X2fYsGFq1KiR4uPj5efn5/zM/U3169eXpFuWu7NrBcUa9V6yTnx74zJ9TULqaOWIzqrty79gAIDKYZMwAAAucE8dX/3f3v/9/POMzUdUXFIzrv66ePFiZWdnKzY2VmFhYc7b+vXrnY956qmntGTJEs2dO1ft2rXTn/70J/3tb39TTEyM8zFnz57VhQsXzCihShQVl2jc2lSlnr0iSWoQYNd7I7soqI6vucEAAB6B68z/gDtfkxYA4N6KSwwNWLRXRzJvXBJt1pNtNfSntz+M3JXoXbdn5s/FMAxN/dthrf/ixkn4Auw+Wv/Lrmrt4P0BANyZR11n/vTp0xo1apSio6Pl7++vpk2basaMGSooKCj1uB07dqhbt24KCAhQWFiYpkyZoqKiIpNSAwBqGm8vm/7fgDbO+29uP6bLeQV3eAY81fykDOcg7+vtpaXDOjLIAwBcyhLD/NGjR1VSUqKlS5fqyJEjmj9/vpYsWaJp06Y5H3Po0CH17dtXffr00cGDB7Vu3Tpt2bJFU6dONTE5AKCm6dQ4SE/e75AkXfmuUG8mHbvLM+BpPth/Rgt2Hpck2WzSvEHt1a1p2S6vBwBAWVn2MPvExEQtXrxYJ0+elCRNmzZNSUlJSk5Odj5m06ZNev7555WVlaWAgIAyvS6HKgIAKutfOdf1szd2K6+gWF42aeu4GLVxBFbZ96N33Z4ZP5eP0i/oV39O1c3/rmb2b60R3aPv/CQAAL7How6zv53s7GwFBf33LMH5+fny8yt9TV9/f39dv35dKSkpP/o6+fn5ysnJKXUDAKAyGtbz07iHm0uSSgxp5pYjsui2c5TD5yf/rfHr0pyD/K9imzLIAwCqjCWH+RMnTmjhwoUaPXq0c9mjjz6qffv2ae3atSouLtb58+f1hz/8QZLueGbc+Ph4BQYGOm9Wum4tAMB9jewerSYhdSRJyacva3NapsmJUJWOXszRi+9/oYKiEknSMx3C9ZtHW9zlWQAAVJypw/zMmTNls9nuePviiy9KPSczM1N9+vTRwIED9eKLLzqX9+7dW4mJiRo9erTsdrt+8pOfqF+/fpIkb2/vH83w6quvKjs723k7d+5c1RQLAKhRfH289Pv+rZ3347d9pbx8Tsrqic5fuabhKw4o9/qN9ze2RQMlPNNONpvN5GQAAE9m6mfmL126pEuXLt3xMY0bN3YePp+ZmalevXrpwQcf1KpVq+Tldeu2CMMwdOHCBd1zzz06ffq0WrdurQMHDqhz585lysTnDgEArvTie1/o46/+JUka97NmeqW36/fW0rturzp+LpfzCvSLJft04ts8SVL7iPpa+9KDqu3rUyXfDwDg+crav0ztNCEhIQoJKdvZXc+fP69evXqpY8eOWrly5W0HeUmy2WxyOG6cRXjt2rWKiIhQhw4dXJYZAIDymN6vlfZkZKmw2NCyT05qUOcIhd9T2+xYcIFrBcUa9V6yc5BvElJHK0d0ZpAHAFQLS3xmPjMzU7GxsYqIiNAbb7yhb7/9VhcvXtTFixdLPS4xMVGHDx/WkSNHNGvWLCUkJGjBggV3PMweAICqFB1SRyO6NZYk5ReVKH7bUXMDwSWKiks0bm2qUs9ekSQ1CLDrvZFdFFTH19xgAIAawxKbjrdv367jx4/r+PHjCg8PL/W1739KYNu2bZo9e7by8/PVvn17bd68WY899lh1xwUAoJRxDzfXhtTz+ndegf7/oQsa3vU/6hIddPcnwi0ZhqHfbkrXx19lSZLq2n206v90VkQQR1wAAKqPJfbMjxgxQoZh3Pb2fTt37tSVK1d07do17d+/n0EeAOAW6vnVKvVZ+dc+PKKSEi5VZ1XzkzK0LvnGCXN9vb20bFhHtXEEmpwKAFDTWGKYBwDA6gZ1jlCrsBsnsUk/n6O/pnxjciJUxAf7z2jBzuOSJJtNmjeovbo1Ldv5fwAAcCWGeQAAqoG3l02/f/y/l6qb+/djyr1eaGIilNdH6Rf0+83pzvszHm+tx+9zmJgIAFCTMcwDAFBNujYNVp82oZKkS1fz9fauEyYnQll9fvLfGr8uTTc/4fer2KYa0T3a3FAAgBqNYR4AgGo0rW8r+XrfaL8r9p5S5pVrJifC3RiGodn/85UKikokSc90CNdvHm1xl2cBAFC1GOYBAKhGkcG1NapHtCKDamvB8w8oLNDP7Ei4C5vNpneHd1a7RoGKbdFACc+0k81mMzsWAKCGs8Sl6QAA8CTjf9ZcEx5uLr9a3mZHQRk1CLBr7cs/lZdNquXNvhAAgPkY5gEAqGb+vgzxVlTXzr9NAAD3waZlAAAAAAAshmEeAAAAAACLYZgHAAAAAMBiGOYBAAAAALAYhnkAAAAAACyGYR4AAAAAAIthmAcAAAAAwGIY5gEAAAAAsBiGeQAAAAAALIZhHgAAAAAAi2GYBwAAAADAYhjmAQAAAACwGIZ5AAAAAAAshmEeAAAAAACLYZgHAAAAAMBiGOYBAAAAALAYH7MDuBvDMCRJOTk5JicBAKBsbvasmz0MN9DTAQBWVNa+zjD/A7m5uZKkiIgIk5MAAFA+ubm5CgwMNDuG26CnAwCs7G593WawGb+UkpISZWZmKiAgQDabrVKvlZOTo4iICJ07d0716tVzUUJzeVpN1OP+PK0m6nF/VqzJMAzl5ubK4XDIy4tP0N3kyp4uWfN34048rR7J82qiHvfnaTVRj3soa19nz/wPeHl5KTw83KWvWa9ePUv98pSFp9VEPe7P02qiHvdntZrYI3+rqujpkvV+N+7G0+qRPK8m6nF/nlYT9ZivLH2dzfcAAAAAAFgMwzwAAAAAABbDMF+F7Ha7ZsyYIbvdbnYUl/G0mqjH/XlaTdTj/jyxJriGp/1ueFo9kufVRD3uz9Nqoh5r4QR4AAAAAABYDHvmAQAAAACwGIZ5AAAAAAAshmEeAAAAAACLYZgHAAAAAMBiGOar0DvvvKPo6Gj5+fmpY8eO+vTTT82OVCHx8fHq3LmzAgICdO+99+rJJ5/UsWPHzI7lMvHx8bLZbJo4caLZUSrl/PnzGjJkiIKDg1W7dm3df//9SklJMTtWhRQVFem3v/2toqOj5e/vryZNmui1115TSUmJ2dHK7JNPPlH//v3lcDhks9m0adOmUl83DEMzZ86Uw+GQv7+/YmNjdeTIEXPClsGd6iksLNSUKVPUrl071alTRw6HQ8OGDVNmZqZ5ge/ibu/P9/3yl7+UzWbTW2+9VW354H48padL9HUr8KSeLlm/r3taT5fo657S1xnmq8j69es1ceJETZ8+XQcPHlSPHj302GOP6ezZs2ZHK7c9e/ZozJgx2r9/v5KSklRUVKTevXsrLy/P7GiVlpycrGXLlum+++4zO0qlXL58Wd27d1etWrW0bds2/fOf/9Sbb76p+vXrmx2tQl5//XUtWbJEixYt0ldffaW5c+cqMTFRCxcuNDtameXl5al9+/ZatGjRbb8+d+5czZs3T4sWLVJycrJCQ0P185//XLm5udWctGzuVM93332n1NRU/e53v1Nqaqo2bNigjIwMDRgwwISkZXO39+emTZs26fPPP5fD4aimZHBHntTTJfq6u/O0ni5Zv697Wk+X6Ose09cNVIkuXboYo0ePLrWsZcuWxtSpU01K5DpZWVmGJGPPnj1mR6mU3Nxco3nz5kZSUpLx0EMPGRMmTDA7UoVNmTLFiImJMTuGy/Tr188YOXJkqWVPP/20MWTIEJMSVY4kY+PGjc77JSUlRmhoqJGQkOBcdv36dSMwMNBYsmSJCQnL54f13M6BAwcMScaZM2eqJ1Ql/Fg933zzjdGoUSMjPT3diIqKMubPn1/t2eAePLmnGwZ93d14Wk83DM/q657W0w2Dvm5l7JmvAgUFBUpJSVHv3r1LLe/du7f27dtnUirXyc7OliQFBQWZnKRyxowZo379+umRRx4xO0qlbdmyRZ06ddLAgQN177336oEHHtDy5cvNjlVhMTEx2rFjhzIyMiRJX375pfbu3au+ffuanMw1Tp06pYsXL5ZaR9jtdj300EMesY6QbqwnbDabZfcklZSUaOjQoYqLi1ObNm3MjgMTeXpPl+jr7sbTerrk2X29JvR0ib7urnzMDuCJLl26pOLiYjVs2LDU8oYNG+rixYsmpXINwzA0efJkxcTEqG3btmbHqbB169YpNTVVycnJZkdxiZMnT2rx4sWaPHmypk2bpgMHDmj8+PGy2+0aNmyY2fHKbcqUKcrOzlbLli3l7e2t4uJizZ49W88//7zZ0Vzi5nrgduuIM2fOmBHJpa5fv66pU6dq8ODBqlevntlxKuT111+Xj4+Pxo8fb3YUmMyTe7pEX3dHntbTJc/u657e0yX6ujtjmK9CNput1H3DMG5ZZjVjx47VoUOHtHfvXrOjVNi5c+c0YcIEbd++XX5+fmbHcYmSkhJ16tRJc+bMkSQ98MADOnLkiBYvXmzJxr9+/XqtXr1aa9asUZs2bZSWlqaJEyfK4XBo+PDhZsdzGU9cRxQWFuq5555TSUmJ3nnnHbPjVEhKSor++Mc/KjU11fLvB1zHE/9eJfq6O/K0ni7VjL7uqesI+rp74zD7KhASEiJvb+9btthnZWXdstXOSsaNG6ctW7Zo165dCg8PNztOhaWkpCgrK0sdO3aUj4+PfHx8tGfPHi1YsEA+Pj4qLi42O2K5hYWFqXXr1qWWtWrVyrInZ4qLi9PUqVP13HPPqV27dho6dKgmTZqk+Ph4s6O5RGhoqCR53DqisLBQzz77rE6dOqWkpCTLbr3/9NNPlZWVpcjISOc64syZM3rllVfUuHFjs+OhmnlqT5fo6+7K03q65Nl93VN7ukRftwKG+Srg6+urjh07KikpqdTypKQkdevWzaRUFWcYhsaOHasNGzZo586dio6ONjtSpTz88MM6fPiw0tLSnLdOnTrphRdeUFpamry9vc2OWG7du3e/5bJCGRkZioqKMilR5Xz33Xfy8iq9evL29rbMJWzuJjo6WqGhoaXWEQUFBdqzZ48l1xHSfxv+119/rY8//ljBwcFmR6qwoUOH6tChQ6XWEQ6HQ3Fxcfr73/9udjxUM0/r6RJ93d15Wk+XPLuve2JPl+jrVsFh9lVk8uTJGjp0qDp16qSuXbtq2bJlOnv2rEaPHm12tHIbM2aM1qxZo82bNysgIMC55TEwMFD+/v4mpyu/gICAWz4XWKdOHQUHB1v284KTJk1St27dNGfOHD377LM6cOCAli1bpmXLlpkdrUL69++v2bNnKzIyUm3atNHBgwc1b948jRw50uxoZXb16lUdP37cef/UqVNKS0tTUFCQIiMjNXHiRM2ZM0fNmzdX8+bNNWfOHNWuXVuDBw82MfWPu1M9DodDv/jFL5SamqoPP/xQxcXFzvVEUFCQfH19zYr9o+72/vzwn5ZatWopNDRULVq0qO6ocAOe1NMl+rq787SeLlm/r3taT5fo6x7T1807kb7ne/vtt42oqCjD19fX6NChg2Uv+SLptreVK1eaHc1lrHwJm5u2bt1qtG3b1rDb7UbLli2NZcuWmR2pwnJycowJEyYYkZGRhp+fn9GkSRNj+vTpRn5+vtnRymzXrl23/bsZPny4YRg3LmUzY8YMIzQ01LDb7UbPnj2Nw4cPmxv6Du5Uz6lTp350PbFr1y6zo9/W3d6fH/KUS9ig4jylpxsGfd0KPKmnG4b1+7qn9XTDoK97Sl+3GYZhuHLjAAAAAAAAqFp8Zh4AAAAAAIthmAcAAAAAwGIY5gEAAAAAsBiGeQAAAAAALIZhHgAAAAAAi2GYBwAAAADAYhjmAQAAAACwGIZ5AAAAAAAshmEegEvExsZq4sSJZscAAACVRE8HrIFhHgAAAAAAi2GYBwAAAADAYhjmAZRbXl6ehg0bprp16yosLExvvvlmqa+vXr1anTp1UkBAgEJDQzV48GBlZWVJkgzDULNmzfTGG2+Uek56erq8vLx04sQJSdLMmTMVGRkpu90uh8Oh8ePHV09xAADUIPR0wLoY5gGUW1xcnHbt2qWNGzdq+/bt2r17t1JSUpxfLygo0KxZs/Tll19q06ZNOnXqlEaMGCFJstlsGjlypFauXFnqNVesWKEePXqoadOm+utf/6r58+dr6dKl+vrrr7Vp0ya1a9euOksEAKBGoKcD1mUzDMMwOwQA67h69aqCg4P1/vvva9CgQZKk//znPwoPD9fLL7+st95665bnJCcnq0uXLsrNzVXdunV14cIFRUREaN++ferSpYsKCwvVqFEjJSYmavjw4Zo3b56WLl2q9PR01apVq5orBACgZqCnA9bGnnkA5XLixAkVFBSoa9euzmVBQUFq0aKF8/7Bgwf1xBNPKCoqSgEBAYqNjZUknT17VpIUFhamfv36acWKFZKkDz/8UNevX9fAgQMlSQMHDtS1a9fUpEkTvfTSS9q4caOKioqqqUIAAGoGejpgbQzzAMrlbgfz5OXlqXfv3qpbt65Wr16t5ORkbdy4UdKNQ/VuevHFF7Vu3Tpdu3ZNK1eu1KBBg1S7dm1JUkREhI4dO6a3335b/v7++vWvf62ePXuqsLCw6goDAKCGoacD1sYwD6BcmjVrplq1amn//v3OZZcvX1ZGRoYk6ejRo7p06ZISEhLUo0cPtWzZ0nminO/r27ev6tSpo8WLF2vbtm0aOXJkqa/7+/trwIABWrBggXbv3q1//OMfOnz4cNUWBwBADUJPB6zNx+wAAKylbt26GjVqlOLi4hQcHKyGDRtq+vTp8vK6sW0wMjJSvr6+WrhwoUaPHq309HTNmjXrltfx9vbWiBEj9Oqrr6pZs2alDvFbtWqViouL9eCDD6p27dr64IMP5O/vr6ioqGqrEwAAT0dPB6yNPfMAyi0xMVE9e/bUgAED9MgjjygmJkYdO3aUJDVo0ECrVq3SX/7yF7Vu3VoJCQm3XLLmplGjRqmgoOCWLfj169fX8uXL1b17d913333asWOHtm7dquDg4CqvDQCAmoSeDlgXZ7MHYJrPPvtMsbGx+uabb9SwYUOz4wAAgAqipwPVj2EeQLXLz8/XuXPn9PLLLyssLEx//vOfzY4EAAAqgJ4OmIfD7AFUu7Vr16pFixbKzs7W3LlzzY4DAAAqiJ4OmIc98wAAAAAAWAx75gEAAAAAsBiGeQAAAAAALIZhHgAAAAAAi2GYBwAAAADAYhjmAQAAAACwGIZ5AAAAAAAshmEeAAAAAACLYZgHAAAAAMBi/hcyCNBXdlMysgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1000x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=plt.subplots(1,2,figsize=[10,5],constrained_layout=True)\n",
    "ax[0].plot(kt.time,kt.zML*-1,lw=2,c='tab:blue')\n",
    "\n",
    "ax[0].set_xlabel('days')\n",
    "ax[0].set_ylabel('MLD (m)')\n",
    "# ax[0].legend()\n",
    "\n",
    "ax[1].plot(kt.time,kt.d-1000,lw=2,c='tab:blue')\n",
    "\n",
    "ax[1].set_xlabel('days')\n",
    "ax[1].set_ylabel('density (kg/m3)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}